divergence of a vector field on a manifold

Let's work backwards.

The product rule gives $$\begin{align*} \frac{1}{\sqrt{\det g}}\partial_i\left(\sqrt{\det g}\,V^i \right) & = \frac{1}{\sqrt{\det g}}\left[ \sqrt{\det g}\ \partial_iV^i + \partial_i\left(\sqrt{\det g}\right)V^i \right] \\ & = \partial_iV^i + \frac{1}{\sqrt{\det g}}\partial_i\left(\sqrt{\det g}\right)V^i \\ & = \partial_iV^i + \frac{1}{\sqrt{\det g}} \frac{1}{2\sqrt{\det g}}\partial_i\left(\det g\right)V^i \end{align*}$$

We can now apply the formula for the derivative of a determinant to get: $$\begin{align*} \frac{1}{\sqrt{\det g}}\partial_i\left(\sqrt{\det g}\,V^i \right) & = \partial_iV^i + \frac{1}{2} \text{tr}(g^{-1}\partial_ig)\,V^i \\ & = \partial_iV^i + \frac{1}{2} \text{tr}(g^{jk}\partial_ig_{kl})\,V^i \\ & = \partial_iV^i + \frac{1}{2} g^{jk} \partial_i g_{kj}\,V^i \end{align*}$$

Now, it is a general fact that for the Levi-Civita connection, the Christoffel symbols satisfy $$\partial_ig_{kj} = \Gamma^l_{ik}g_{lj} + \Gamma^l_{ij}g_{kl}.$$ Therefore, $$\begin{align*} \frac{1}{\sqrt{\det g}}\partial_i\left(\sqrt{\det g}\,V^i \right) & = \partial_iV^i + \frac{1}{2} g^{jk} \left(\Gamma^l_{ik}g_{lj} + \Gamma^l_{ij}g_{kl} \right)\,V^i \\ & = \partial_iV^i + \frac{1}{2}\left(\Gamma^l_{ik}\delta^k_l + \Gamma^l_{ij}\delta^j_l \right)\,V^i \\ & = \partial_iV^i + \Gamma^l_{il}\,V^i \\ & = Div \cdot V \end{align*}$$ as was to be shown.

$\def\div[0]{\operatorname{div}}$ My favourite reason this formula is true: let $\phi$ be an arbitrary smooth function with compact support contained in a single chart $(U,x)$. Then integrating by parts we get

$$ \int_U \phi \div V d\mu_g= -\int_U V^i \partial_i \phi\, d\mu_g = -\int_U V^i \partial_i \phi \sqrt {\det g}\, dx.$$ Now integrate by parts again, but this time in the chart geometry (i.e. use the product rule for partial derivatives instead of the one for Riemannian divergence, and the measure $dx$ instead of $d \mu_g$):

$$ -\int_U \partial_i \phi\, V^i \sqrt {\det g}\, dx= \int_U \phi\, \partial_i (\sqrt{\det g}\, V^i) dx = \int_U \phi\, \frac 1{\sqrt{\det g}} \partial_i(\sqrt{\det g}\, V^i)\ d\mu_g.$$

Since this is true for every $\phi$ we get the desired formula for $\div V$.

In some sense this is just a weak version of Giuseppe's answer, but seeing it was the first time this formula really clicked for me, so I figured I'd write it down.

The metric volume form is $\omega=\sqrt{|\det(g)|}dx^1\wedge\ldots\wedge dx^n.$
The divergence of $V=V^i\partial_i$ is determined by $(\text{div}V)\omega=d(V\lrcorner\omega)\equiv V(\omega),$ hence we get: $$(\text{div}V)\omega=\left[V^i\partial_i(\sqrt{|\det(g)|})+\sqrt{|\det(g)|}\partial_iV^i\right]dx^1\wedge\ldots\wedge dx^n,$$ Where we used the obvious formula $V(dx^1\wedge\ldots\wedge dx^n)=(\partial_iV^i)dx^1\wedge\ldots\wedge dx^n.$
Therefore $$\text{div}V=\frac{1}{\sqrt{|\det(g)|}}\frac{\partial}{\partial x_i}\left[\sqrt{|\det(g)|}V^i\right]$$

Edit added in reply to Asaf Shachar's comment

Probably, in the original version of my answer, I should have spent a few words justifying the following formula $$V(dx^1\wedge\ldots\wedge dx^n)=(\partial_iV^i)dx^1\wedge\ldots\wedge dx^n.\tag{$\star$}$$ Here, on the left-hand side, as elsewhere in my answer, I am using the expression $V(\eta)$ to denote the Lie derivative, along the vector field $V$, of a differential form $\eta$.

Recall that the vector field $V$ is locally given by $V=V^i\partial_i$ in terms of the holonomic frame $\partial_1,\ldots,\partial_n$ corresponding to the local coordinate system $x^1,\ldots,x^n$, so that $$V(x^i)=V^i\tag{1}.$$ Moreover, since the Lie derivative $V(-)$ is a (degree $0$) derivation of the differential graded algebra $(\Omega^\bullet(X),\wedge,d_{dR})$, we get immediately that $$V(dx^1\wedge\ldots\wedge dx^n)=\sum_{i=1}^n dx^1\wedge\ldots\wedge d(V(x^i))\wedge\ldots\wedge dx^n.\tag{2}$$ Now, from $(1)$ and $(2)$, it follows $(\star)$, as we needed.