Let ${\mathbf {E} }({\mathbf {x} })=[E_{x}(x,y,z),E_{y}(x,y,z),E_{z}(x,y,z)]$ be a smooth (differentiable) three-component vector field on the three dimensional space and $\nabla \cdot \mathbf {E} ={\frac {\partial E_{x}}{\partial x}}+{\frac {\partial E_{y}}{\partial y}}+{\frac {\partial E_{z}}{\partial z}}$ is its divergence then the field divergence integral over the arbitrary three dimensional volume $V$ equals to the integral over the surface $S$ of the field itself projected onto the unite length vector field $\mathbf {n} (\mathbf {x} )$ always perpendicular to the surface and pointing outside the surface which contains this volume or otherwise the inner values of the field divergence make virtually no contributions to the integral over the volume i.e,

$\int _{V}\nabla \cdot \mathbf {E} \,\,dV=\int _{S}\mathbf {E} \cdot d\mathbf {S}$

where $d\mathbf {S} =\mathbf {n} dS$ and the $S$ wraps the $V$.

We can approximate the integral of the divergence over the volume by the finite sum by dividing densely the space inside the volume $V$
into small cubes with the edges $dx=dy=dz$ and the corners $[x_{i},y_{j},z_{k}]$ as well as approximating three of the coordinate derivatives
by their difference quotients. We will keep the edges coordinate names for the convenience even if they are equal.
We get

$\int _{V}\nabla \cdot \mathbf {E} \,\,dV=\sum _{i,j,k}\left[{\frac {E_{x}(x_{i+1},y_{j},z_{k})-E_{x}(x_{i},y_{j},z_{k})}{dx}}+{\frac {E_{y}({x_{i},y_{j+1},z_{k}})-E_{y}(x_{i},y_{j},z_{k})}{dy}}+{\frac {E_{z}(x_{i,j,k+1})-E_{z}(x_{i,j,k})}{dz}}\right]dxdydz+\Theta (dxdydz)$

Let us focus on a single contribution to this sum related to the derivative with respect to a chosen coordinate for example $x$ i.e.
for example $\sum _{i,j,k}{\frac {E_{x}(x_{i+1},y_{j},z_{k})-E_{x}(x_{i},y_{j},z_{k})}{dx}}dxdydz$. For a fixed $j,k$ we have

$\sum _{i}{\frac {E_{x}(x_{i+1},y_{j},z_{k})-E_{x}(x_{i},y_{j},z_{k})}{dx}}dxdydz=\sum _{i}[E_{x}(x_{i+1},y_{j},z_{k})-E_{x}(x_{i},y_{j},z_{k})]dydz$

Note now that because of the alternating signs the vast majority of terms in the right sum vanish and we have

$\sum _{i}[E_{x}(x_{i+1},y_{j},z_{k})-E_{x}(x_{i},y_{j},z_{k})]dydz=[E_{x}(x_{1},y_{j},z_{k})-E_{x}(x_{0},y_{j},z_{k})+E_{x}(x_{2},y_{j},z_{k})-E_{x}(x_{1},y_{j},z_{k})+...+E_{x}(x_{n},y_{j},z_{k})-E_{x}(x_{n-1},y_{j},z_{k})]dydz$

which reduces only to two terms or

$\sum _{i}[E_{x}(x_{i+1},y_{j},z_{k})-E_{x}(x_{i},y_{j},z_{k})]dydz=[E_{x}(x_{n},y_{j},z_{k})-E_{x}(x_{1},y_{j},z_{k})]dydz$

where the bordering $[x_{n},y_{j},z_{k}]$ and $[x_{1},y_{j},z_{k}]$ with the first coordinate obviously depending on the choice of $j$ and $k$ are such that those points are the closed to the
surface $S$ containing the volume $V$.

Also note that while $dydz$ is an infinitesimal (small) element of the surface parallel to the $y-z$ plane and for the unite vector $\mathbf {n} _{x}=[1,0,0]$ perpendicular to it $\mathbf {E} (x_{n},y_{j},z_{k})\cdot \mathbf {n} _{x}=E_{x}(x_{n},y_{j},z_{k})$ and so for the second point the right side is an approximate to the growth $\mathbf {E} \cdot d\mathbf {S}$ of the surface integral $\int _{S}\mathbf {E} \cdot d\mathbf {S}$ i.e.

$E_{x}(x_{n},y_{j},z_{k})dydz=\mathbf {E} \cdot d\mathbf {S} +\Theta (dydz)$,

$-E_{x}(x_{1},y_{j},z_{k})]dydz=\mathbf {E} \cdot d\mathbf {S} +\Theta (dydz)$,

Repeating the estimate for the two other dimensions and coming back to the original sum we get

$\int _{V}\nabla \cdot \mathbf {E} \,\,dV=\sum _{j,k}[E_{x}(x_{n_{j,k}},y_{j},z_{k})-E_{x}(x_{1_{j,k}},y_{j},z_{k})]dydz+\sum _{i,k}[E_{y}(x_{i},y_{n_{i,k}},z_{k})-E_{y}(x_{i},y_{1_{i,k}},z_{k})]dxdz+\sum _{i,j}[E_{z}(x_{i},y_{j},z_{n_{i,j}})-E_{z}(x_{i},y_{j},z_{1_{i,j}})]dxdy+\Theta (dxdydz)$

so the right side is the approximate surface integral (sum over the surfaces of the cubes closest to the surface $S$) of the field itself projected on the outward unit vector field which proves the therem i.e.

$\int _{V}\nabla \cdot \mathbf {E} \,\,dV=\int _{S}\mathbf {E} \cdot d\mathbf {S}$.