I am trying to read some basic constructions on Riemannian manifolds, and I came across the definition of the divergence of a vector field. It is given by the formula,
$$ d( i_{X} d V_g ) = ( \text{div} X ) dV_g. $$
Here $i_{X} d V_g$ denotes the interior multiplication of the Riemannian volume form with $X$, giving us a $(n-1)$ form. $d$ denotes taking the exterior derivative, giving us a $n$ form, and the divergence is specified to be the smooth function multiplying with the resulting $n$ form.
Can anyone she some shed on this formula? I am finding it very hard to motivate this coordinate independent definition.
This is a sketch rather than a complete answer. You did not specify what is for you the intuitive meaning of divergence (usually some version of "flux per unit of volume") so it is a bit hard to answer. But I can think of three possibly helpful things
write the expression in local coordinate and convince yourself that you get the same thing as the divergence you are familiar with - essentially you are going to get $\nabla_\mu X^\mu\ \mathrm{vol}$.
since the volume element is closed, the expression you wrote is equal to $L_X\ \mathrm{vol}$. So, what you have is how the volume element changes when Lie transported along $X$, which is kinda of the dual viewpoint of the flux ingoing/outgoing a fixed volume element perspective.
By Stokes' theorem \begin{equation} \int_V d(i_X\ \mathrm{vol})=\int _{\partial V}i_X \ \mathrm{vol}. \end{equation} Now you have to convince yourself that $i_X \ \mathrm{vol}= *X$ is the flux of $X$ across $\partial V$. Something along the following lines shoud work. Write $\mathrm{vol}=n \wedge W$ with $n$ a 1-form such that the associated vector field is unit normal to $\partial V$ and $W$ is a volume form on $\partial V$. This way $i_X\ \mathrm{vol}$ once restricted to $\partial V$ is equal to $\langle X,n \rangle W $.