An inconsistency between flux through surface and the divergence theorem

Question

There was recently a question about a flux through a surface. The idea is to use the formula:

$$\iint_S \mathbf F \cdot d\mathbf S,$$

with $\mathbf{F} = x \boldsymbol{\hat{\imath}}+ y \boldsymbol{\hat{\jmath}} - 2z \boldsymbol{\hat k}$ and $S$ is the semisphere with radius $a$ and $z≥0$.

Since $S$ is closed, you can use the divergence theorem, whose result is:

$$\iiint_V \nabla \cdot \mathbf F \, dV = -\dfrac{2}{3} \pi a^3.$$

However, there are two answers with different results. The first one uses spherical coordinates and cartesian coordinates to find the flux, which is zero; the second one uses parametrization of the surface in spherical coordinates to find the flux, which is not zero, is $ -2\pi a^3/3$. And just to be sure, the flux through the bottom of the semisphere is zero (when $z=0$):

$$\iint_{S_1} \mathbf F \cdot \boldsymbol{\hat n} \, dS = \iint_{S_1} \mathbf F \cdot (-\boldsymbol{\hat k}) \, dS = \iint_{S_1} 2z \, dS = 0,$$

so we are comparing the flux through the semisphere.

If you want, you can reproduce the answers to check the inconsistency, so my question is: why this inconsistency is happening, with one method you get an answer, with the other you get another answer.

Answer 1

It appears there must be a mistake in your computation in spherical coordinates that found $\iint_S \mathbf F \cdot d \mathbf S = -\frac{2}{3}\pi a^2$. I'm not terribly interested in tracking down where it was, but if the correction to your divergence computation in the comments didn't convince you the answer should be zero, here's another pretty quick way to compute the surface integral:

$$\iint_S \mathbf F \cdot d \mathbf S = \iint_S \mathbf F \cdot \frac{x \hat{\mathbf{i}} + y \hat{\mathbf{j}} + z \hat{\mathbf{k}}}{a} dS = \frac{1}{a}\iint_S (x^2+y^2-2z^2) dS.$$

By symmetry under the reflection $\mathbf x \mapsto -\mathbf x$, this is the same as (denoting by $\tilde{S}$ the full sphere)

$$\iint_S \mathbf F \cdot d \mathbf S = \frac{1}{2a} \iint_\tilde{S} (x^2 +y^2-2z^2) dS.$$

Now, by applying a rotation interchanging the coordinate axes and noting $\tilde{S}$ is invariant under rotations, it's easy to see $$\iint_{\tilde{S}} x^2 dS = \iint_{\tilde{S}} y^2 dS = \iint_{\tilde{S}} z^2 dS,$$ giving the result $\iint_S \mathbf F \cdot d \mathbf S = 0.$ So, the resolution is that there's no inconsistency between the flux integral and the divergence theorem's result.