Solution (source code)

= Solution

The <divergence> of $\mathbf F=\nabla z$ is the <Laplacian>
$$
\nabla\cdot\mathbf F=6xy(x^2+y^2-8).
$$
Its area <integral> over the square is
$$
\int_0^2\int_0^2 6xy(x^2+y^2-8)\,dx\,dy
=6(8+8-32)=\boxed{-96}.
$$
For the boundary <flux>, $\int_0^2A(s)\,ds=-4$. On $x=0$, $\mathbf F\cdot\mathbf n=4A(y)$, giving $-16$; on $x=2$ it is $8A(y)$, giving $-32$. The corresponding horizontal sides give $-16$ and $-32$. The total is again $-96$, verifying the <divergence theorem> with outward normals. With a counterclockwise boundary, this is the same orientation as $M\,dy-N\,dx$.