Divergence formula on manifolds using a physical approach

Question

In rectangular coordinates, one can define the divergence of a vector field $\vec{V}$ as the following limit $$ \text{Div } V = \lim_{\epsilon \rightarrow 0} \frac{1}{\epsilon^3} \oint_{\partial C_{\epsilon}(p)} \vec{V} \cdot d \vec{A} \: \: \: \: (1), $$ where $C_{\epsilon}(p)$ is the cube of length $\epsilon$ and with corners at $p$ and $p+\epsilon(1,1,\ldots,1)$.

In the setting of Riemannian manifolds we have the formula $$ \text{Div } V = \frac{1}{\sqrt{|g|}} \partial_i \left( \sqrt{|g|} V^i \right). \: \: \: \: (2)$$ By the divergence theorem (for manifolds) one also obtains a definition for $\text{Div } V$ as in $(1)$. My question is how to deduce $(2)$ from such a definition.

I'd like to have an answer which doesn't use differential forms or Lie derivatives.

Answer 1

The proof of (2) from (1) is entirely local, that is, it can be carried out in an arbitrarily small neighbourhood of $p$.

EDIT THE FOLLOWING IS WRONG (See comments)

Since all Riemannian manifolds are locally isometric to Euclidean space, if you can prove (2) from (1) in an arbitrary curvilinear coordinate system in $\mathbb R^n$ then you are done. And it seems to me that you have already done this last step. (If not, then you can find it in books in vector calculus or fluid mechanics).