At any point $P$, we define the divergence of a vector field $\FF$, which we write for now as $\Div\FF$, to be the flux of $\FF$ per unit volume “at” $P$, which is the (limiting value of the) flux per unit volume out of a small box around $P$.
Furthermore, we can chop up an arbitrary (closed) box into a suitable combination of (small) rectangular boxes. The outward fluxes through all interior sides will cancel, leaving just the flux out of the original box. Thus, \begin{equation} \Int_{\rm box} \FF \cdot d\SS = \Int_{\rm inside} \hskip-5pt \Div\FF \>\> dV \label{DivThm} \end{equation} This is the Divergence Theorem! But what is the divergence?
In order to find the divergence, we must compute the flux. Consider a (possibly large) rectangular box whose sides are parallel to the coordinate planes. What is the flux of $\FF$ out of this box? Consider first the vertical contribution, namely the flux up through the top plus the flux down through the bottom. These two sides each have area $dA=dx\,dy$, so we get \begin{eqnarray*} \Int_{\rm top+bottom} \hskip-7pt\FF \cdot d\SS &=& \DInt{R} \FF \bigg|_{\rm top}\kern-3pt\cdot\kk\> dx\,dy - \DInt{R} \FF \bigg|_{\rm bottom}\kern-10pt\cdot\kk\> dx\,dy \cr &=& \DInt{R} \FF \bigg|_{\rm bottom}^{\rm top}\kern-10pt\cdot\kk\> dx\,dy \end{eqnarray*} where $R$ denotes the common domain of the top and bottom faces in terms of $x$ and $y$. Note the minus sign, which is really due to taking the dot product of $\FF$ with $-\kk$, which is the outward unit normal vector on the bottom face. But along any vertical line $C$ from the bottom face to the top face, we have $$ \FF \bigg|_{\rm bottom}^{\rm top}\kern-10pt\cdot\kk = F_z \bigg|_{\rm bottom}^{\rm top} = \Lint dF_z = \Lint \grad F_z \cdot d\rr = \Lint \Partial{F_z}{z} \,dz $$ where the second equality is the Fundamental Theorem of Calculus (“adding small changes yields the total change”), the next equality is the master formula (1) of Section 3, and the final equality follows since $d\rr=dz\,\kk$ for a vertical path.
Inserting this last expression into the one before results in $$ \Int_{\rm top+bottom} \hskip-7pt\FF \cdot d\SS = \TInt{V} \Partial{F_z}{z} \>dz\,dx\,dy $$ where $V$ denotes the volume inside the box. Repeating this argument using the remaining pairs of faces, it follows that the total flux out of the box is $$ \hbox{total flux} = \Int_{\rm box} \FF \cdot d\SS = \Int_V \left( \Partial{F_x}{x} + \Partial{F_y}{y} + \Partial{F_z}{z} \right) \> dV $$ We can use this to compute the flux per unit volume through a small box around $P$. Suppose the box is small enough that the expression in parentheses is nearly constant. We then get $$ \hbox{total flux} \approx \left( \Partial{F_x}{x} + \Partial{F_y}{y} + \Partial{F_z}{z} \right) \Bigg|_P \> (\hbox{volume of $V$}) $$ from which we see that the flux per unit volume, and hence the divergence, is given in rectangular coordinates by \begin{equation} \Div\FF = \Partial{F_x}{x} + \Partial{F_y}{y} + \Partial{F_z}{z} \label{Divergence} \end{equation}
The formula (\ref{Divergence}) is often taken as the definition of the divergence, in which case the preceding calculation is used to prove the Divergence Theorem (\ref{DivThm}); we prefer to view (\ref{DivThm}) as the definition, from which (\ref{Divergence}) can be derived. This emphasizes the important fact that $${\rm divergence} = {{\rm flux}\over{\rm unit volume}}$$ For this reason, points where $\Div\FF>0$ are referred to as sources, and points where $\Div\FF<0$ are referred to as sinks.
It is often easier to work in coordinates adapted to the given problem. To use the Divergence Theorem in, say, spherical coordinates one must have an appropriate expression for the divergence. % For example, a spherically symmetric, radial vector field of the form $$\HH = h(r) \,\rhat = h(r) \,{\rr\over r}$$ has divergence $$\Div\HH = {1\over r^2} {\partial\over\partial r} \left(r^2 h(r)\right)$$ as can be checked by looking up the full formula for the divergence in spherical coordinates (e.g. in the book Div, Grad, Curl, and All That). The electric field of a point charge has this form. But where did this formula come from? You can derive it yourself by computing the flux of $\HH$ out of a small box adapted to spherical coordinates, that is, whose sides are surfaces with one of the spherical coordinates held constant.