For the sake of argument, suppose we have some closed surface $S$ in a vector field $\vec{F}(x, y, z)$. Since $S$ is a closed surface, it forms the boundary of some 3D region $D$, where $D$ has some finite but non-zero volume. If we compute the surface integral of $\vec{F}$ across $S$ (dotting with the surface's normal vectors), we get the flux of $\vec{F}$ across the surface $S$.
Now, let's pick an arbitrary point, $P$, such that $P{\in}D$. Suppose that we perform a kind of uniform compression of $S$ such that all points on $S$ have their distance to $P$ reduced by a certain scaling factor $\lambda$. That is, if $P_S$ represents some general point on $S$, and $\tilde{P}_S$ represents that point after this transformation, then $\text{dist}(\tilde{P_S},P)=\frac{\text{dist}(P_S,P)}{\lambda}$. The question then is what happens to the surface integral as we let $\lambda\rightarrow\infty$. It seems reasonable to me that this limit should approach $(\nabla\cdot\vec{F})(P)$, but I cannot even figure out how one would rigorously write this limit of a scaling transformation down, notwithstanding how to go about solving it.
EDIT: I think I figured out a way to at least rigorously define the transformation I am talking about. The definition I came up with is that $T_\lambda(S)=\{\tilde{P_S}\in\mathbb{R}^3\mid{}\exists\space{}P_S\in{}S\mid{}\tilde{P_S}=\frac{1}{\lambda}P_S+\frac{\lambda-1}{\lambda}P\}$, where $P$ is a point in the region enclosed by $S$. The limit as $\lambda\rightarrow\infty$ is $P$, which makes sense. I am still working on the problem of applying the limit of this transformation to the surface integral though.