Below is a very concise proof of Gauss's theorem (from the book Vector and Tensor Analysis with Applications by Borisenko, Tarapov). Unfortunately, I'm having trouble understanding it, despite staring at it for a few days. I'd really appreciate it if someone would explicitly show a couple of the intermediate steps. I've marked these with $(*)$:
Gauss's Theorem:
Given a volume $V$ bounded by a closed surface $S$, suppose the functions $$P(x_1,x_2,x_3),\ \ Q(x_1,x_2,x_3),\ \ R(x_1,x_2,x_3)$$ and their derivatives $$\frac{\partial P}{\partial x_1},\ \frac{\partial Q}{\partial x_2},\ \frac{\partial R}{\partial x_3}$$ are continuous on $V \cup S$. Then $$\iiint \limits_V \left( \frac{\partial P}{\partial x_1} + \frac{\partial Q}{\partial x_2} + \frac{\partial R}{\partial x_3} \right)dV = \iint \limits_S \left[ P\cos(\mathbf{n},x_1) + Q\cos(\mathbf{n},x_2)+ R\cos(\mathbf{n},x_3)\right]dS,$$ where $\mathbf{n}$ is the unit exterior normal to $S$.
Proof: Suppose no line parallel to the $x_1$-axis intersects $S$ in more than two points $M^\prime$ and $M^{\prime\prime}.$ Then, if $S_{23}$ is the projection of $S$ onto the $x_2x_3$-plane, we have $$(*)\iiint\limits_V \frac{\partial P}{\partial x_1}dV = \iint\limits_{S_{23}} \left( \int\frac{\partial P}{\partial x_1}dx_1\right)dS_{23} = \iint \limits_{S_{23}}\left[ P(M^\prime) - P(M^{\prime\prime})\right]dS_{23}.$$ But the element $dS_{23}$ of the projection $S_{23}$ can be expressed in terms of the elements of the surface $S$ at the points $M^\prime$ and $M^{\prime\prime}$: $$dS_{23} = dS(M^\prime)\cos[\mathbf{n}(M^\prime),x_1] = -dS(M^{\prime\prime}) \cos[\mathbf{n}(M^{\prime\prime}),x_1].$$ Therefore $$ (**) \iiint\limits_V \frac{\partial P}{\partial x_1}dV = \iint \limits_S P(M)\cos[\mathbf{n}(M),x_1]dS(M),$$ where $M$ is a variable point of the surface $S$. The formulas $$\iiint\limits_V \frac{\partial Q}{\partial x_2}dV = \iint \limits_S Q\cos[\mathbf{n},x_2]dS,$$ $$\iiint\limits_V \frac{\partial R}{\partial x_3}dV = \iint \limits_S R\cos[\mathbf{n},x_3]dS$$ are proved in the same way, provided no line parallel to the $x_2$ or $x_3$-axis intersects $S$ in more than two points. Adding these formulas, we obtain $$\iiint \limits_V \left( \frac{\partial P}{\partial x_1} + \frac{\partial Q}{\partial x_2} + \frac{\partial R}{\partial x_3} \right)dV = \iint \limits_S \left[ P\cos(\mathbf{n},x_1) + Q\cos(\mathbf{n},x_2)+ R\cos(\mathbf{n},x_3)\right]dS.$$
The first $(*)$ is just Fubini's theorem: we decompose the volume $V$ into its projection and the fibres of this projection. Write $\pi :V \to S_{23}$ for the projection, and take a point $s\in S_{23}$. The integral $\int \partial P /\partial x_1 dx_1$ has bounds the two endpoints of the closed interval $\pi^{-1}(s)$, hence the formula $\int \partial P /\partial x_1 dx_1 = P(M') - P(M'')$.
The second $(**)$ is obtained by the change of variables formula for surface integrals, to express the integral $\iint_{S_{23}}$ as an integral $\iint_S$. Under assumption, locally the projection $S \to S_{23}$ has two different sections, and we can lift the integral $\iint_{S_{23}}$ in two different ways to $S$. We have:
$$\iint_{S_{23}}P(M') dS = \iint_{S'} P(M') \cos(n(M'), x) dS'$$ $$\iint_{S_{23}}P(M'') dS = - \iint_{S''} P(M'') \cos(n(M''), x) dS''$$ where I used $S'$ and $S''$ to denote respectively the "upper" half of $S$ (i.e. the part that you would see by looking at $S$ from high above on the $x_1$-axis) and the bottom half. Here we are viewing $M'$ and $M''$ as functions of $s \in S_{23}$. Taking the difference of these two equations and using $(*)$, the two minus signs cancel and we have:
$$\iiint(...) = \iint_{S_{23}}(P(M')-P(M'')) = \left\{\iint_{S'}(\dots) + \iint_{S''}(\dots)\right\} = \iint_S(\dots).$$