Theorem (Divergence Theorem): For a simple solid region $E$ with outward-oriented boundary $S = \partial E$, $\oiint_S \mathbf{F} \cdot d\mathbf{S} = \iiint_E \nabla \cdot \mathbf{F}\,dV$

Write $\mathbf{F} = (P, Q, R)$. It suffices to prove: $\oiint_S R\,\hat{\mathbf{k}} \cdot d\mathbf{S} = \iiint_E \frac{\partial R}{\partial z}\,dV \quad (*)$ and the analogous statements for $P$ and $Q$, then add.

Step 1: Set up the region.

Assume $E$ is a type-$z$ region: $E = \{(x,y,z) : (x,y) \in D, \; g_1(x,y) \leq z \leq g_2(x,y)\}$ where $D$ is the projection of $E$ onto the $xy$-plane, and $g_1, g_2$ are continuous functions.

The boundary $S$ consists of three parts:

• $S_1$ (bottom): $z = g_1(x,y)$, outward normal points downward
• $S_2$ (top): $z = g_2(x,y)$, outward normal points upward
• $S_3$ (lateral): the vertical sides connecting $S_1$ and $S_2$

Step 2: Compute the volume integral.

$\iiint_E \frac{\partial R}{\partial z}\,dV = \iint_D \left[\int_{g_1(x,y)}^{g_2(x,y)} \frac{\partial R}{\partial z}\,dz\right] dA = \iint_D [R(x,y,g_2(x,y)) - R(x,y,g_1(x,y))]\,dA$

Step 3: Compute the flux integral.

On $S_2$ (top, $z = g_2(x,y)$, outward normal upward): parametrize by $(x,y) \in D$ with $\mathbf{r}(x,y) = (x, y, g_2(x,y))$: $d\mathbf{S} = (-\partial_x g_2, -\partial_y g_2, 1)\,dA$ The $R\hat{\mathbf{k}}$ contribution: $R(x,y,g_2) \cdot 1 \cdot dA$, so $\iint_{S_2} R\,\hat{\mathbf{k}} \cdot d\mathbf{S} = \iint_D R(x,y,g_2(x,y))\,dA$

On $S_1$ (bottom, $z = g_1(x,y)$, outward normal downward): $d\mathbf{S} = (\partial_x g_1, \partial_y g_1, -1)\,dA$ The $R\hat{\mathbf{k}}$ contribution: $R(x,y,g_1) \cdot (-1) \cdot dA$, so $\iint_{S_1} R\,\hat{\mathbf{k}} \cdot d\mathbf{S} = -\iint_D R(x,y,g_1(x,y))\,dA$

On $S_3$ (lateral surface): the outward normal is horizontal (no $\hat{\mathbf{k}}$ component), so $\iint_{S_3} R\,\hat{\mathbf{k}} \cdot d\mathbf{S} = 0$

Step 4: Combine.

$\oiint_S R\,\hat{\mathbf{k}} \cdot d\mathbf{S} = \iint_D R(x,y,g_2)\,dA - \iint_D R(x,y,g_1)\,dA = \iiint_E \frac{\partial R}{\partial z}\,dV$

This proves $(*)$. The analogous results for $P$ (using type-$x$ regions) and $Q$ (using type-$y$ regions) follow by identical arguments with permuted coordinates.

Adding all three: $\oiint_S \mathbf{F} \cdot d\mathbf{S} = \iiint_E \left(\frac{\partial P}{\partial x} + \frac{\partial Q}{\partial y} + \frac{\partial R}{\partial z}\right) dV = \iiint_E \nabla \cdot \mathbf{F}\,dV \quad \square$