I want to calculate the flux of the vector field f[x,y,z] = {x,y,z} across the sphere of radius “a” centered at the origin with outward-Orientated Surface normal

In[]:=
VectorPlot3D[{x,y,z},{x,-2,2},{y,-2,2},{z,-2,2}]
The answer appears to be
"4
3
a
π
” from the calculations below as well as in the last calculation using SurfaceIntegrate. However, I cannot seem to get the Integrate function to give me the correct answer and I do not understand why.
In[]:=
Div[{x,y,z},{x,y,z}]
Out[]=
3
In[]:=
3
4
3
π
3
(a)
Out[]=
4
3
a
π
In[]:=
Integrate[Div[{x,y,z},{x,y,z}],{x,y,z}∈Sphere[{0,0,0},a]]
Out[]=
12
2
a
π
In[]:=
Integrate[
∇
{x,y,z}
.{x,y,z},{x,y,z}∈Sphere[{0,0,0},a]]
Out[]=
12
2
a
π
In[]:=
SurfaceIntegrate[{x,y,z},{x,y,z}∈Sphere[{0,0,0},a]]
Out[]=
4
3
a
π