Solution 1:

You can always write $ z^2 = z^2 \hat{n}.\hat{n} $ where $\hat{n} $ is the unit normal then using divergence theorem your expression becomes $$ \iiint_{|z|\leq 1} \nabla.(z^2\hat{n})dV $$

Solution 2:

Indeed we can use divergence theorem if we could find an $F$ such that: $$ \iint_S g\,dS = \iint_S F\cdot n\,dS.\tag{1} $$

Though the example you gave is not very illustrative it seems. Say the unit cube is $D = [0,1]\times[0,1]\times [0,1]$. The faces perpendicular to the $x$-axis have unit outward normals $(\pm 1,0,0)$, if we want $F\cdot (1,0,0) = z^2$ on the face $\{x = 1\}\cap D$, while $F\cdot (-1,0,0) = z^2$ on the face $\{x = 0\}\cap D$ as well?

Here we can't find a well-defined $F$ such that $F\cdot n = z^2$ on all six faces of the unit cube, and we have to evaluate directly: $$ \mathrm{Top} + \mathrm{Bottom} = \iint_{[0,1]\times[0,1]\times \{1\}} 1 \,dxdy + \iint_{[0,1]\times[0,1]\times \{-1\}} 0 \,dxdy, \\ \mathrm{Right} + \mathrm{Left} = \iint_{[0,1]\times\{1\}\times [0,1]} z^2 \,dzdx + \iint_{[0,1]\times \{-1\}\times[0,1]} z^2 \,dzdx, \\ \mathrm{Front} + \mathrm{Back} = \iint_{\{1\}\times[0,1]\times [0,1]} z^2 \,dydz + \iint_{\{-1\}\times[0,1]\times [0,1]} z^2 \,dydz. $$


Based on my teaching experience, the problem to illustrate (1) is actually taking $S$ to be the unit 2-sphere: $$ S = S^2 = \{(x,y,z): x^2 + y^2 + z^2 = 1\}, $$ so that the unit outward normal is $n=(x,y,z)$: $$ \iint_S z^2 \,dS = \iint_S (0,0,z)\cdot (x,y,z) \,dS = \iiint_{B(0,1)} \nabla\cdot (0,0,z)\,dV = |B(0,1)|. $$