Divergence in terms of Levi-Civita connection

The divergence of a vector field $X$ on a manifold $M$ is defined usually as the function $\text{Div}(.)$ such that $(\text{Div} X) \;\mu =L_X \mu$ for $\mu$ a volume form.

I know that there is also an alternative expression for $\text{Div} \; X$ in terms of the covariant derivative:

$\text{div} \; X= \text{trace}(\nabla X)= \nabla^aX_a.$

Can anyone explain to me how to prove this equivalence?

To be precise, given an oriented metric $g$ on an $n$-manifold with volume form $\mu_{b_1 \cdots b_n}$, the divergence of a vector field w.r.t. $\mu$ coincides with $\nabla^a X_a$, where $\nabla$ is the Levi-Civita connection of $g$.

One way to see this is as follows: We can write the Lie derivative of a $k$-tensor $\phi$ on $M$ w.r.t. $X$ in terms of $\nabla$ (any torsion-free connection will do) as $$(\mathcal{L}_X \phi)_{b_1 \cdots b_k} = X^c \nabla_c \phi_{b_1 \cdots b_k} + \sum_{i = 1}^k (\nabla_{b_i} X^c) \phi_{b_1 \cdots c \cdots b_k}.$$

If we take $\phi$ to be $\mu$, then the first term is zero (the volume form of a Riemannian metric is parallel with respect to the Levi-Civita connection). Now, we may write the second term as $$-(n + 1) \nabla_{[c} X^c \mu_{b_1 \cdots b_n]} + \nabla_c X^c \mu_{b_1 \cdots b_n}.$$ Here the bracket $[ \cdots ]$ indicates antisymmetrization over the enclosed indices, so the first term in this expression is the trace of a tensor of type $(n + 1, 1)$ which is skew in its $n + 1$ lower indices. But any such tensor on an $n$-manifold is zero, so only the second term survives, giving (now in classical notation, for emphasis) $$\mathcal{L}_X \mu = (\mathrm{tr} \nabla X) \mu = (\operatorname{div} X) \mu.$$

Though I probably would not have found this without Travis proof above, there is actually a way to show this without having to recourse to coordinate calculations. As far as I can tell, this is the most elegant way to show this:

It is correct that the covariant derivative has to be torsion-free for the calculation to work out like this: $$[X,Y] = \nabla_XY - \nabla_Y X$$

Using an argument I found in O'Neill's book "Semi-Riemannian Geometry", the first step in the proof is to realize that $\mu$ is a top degree form and so is $\mathcal{L}_X \mu$. Therefore, there has to exist a smooth (possibly zero) function $f$ such that $$\mathcal{L}_X \mu = f \, \mu$$ This is true regardless of the torsion of $\nabla$. Introducing a local frame field $(Y_1, \dots, Y_n)$ with $$\mu (Y_1, \dots, Y_n) =1 \, ,$$ we locally find \begin{align} f &= \left( \mathcal{L}_X \mu \right) (Y_1, \dots, Y_n) \\ &= \mathcal{L}_X \left( \mu (Y_1, \dots, Y_n) \right) - \mu (\mathcal{L}_X Y_1, Y_2, \dots, Y_n) - \dots - \mu (Y_1, \dots, Y_{n-1}, \mathcal{L}_X Y_n) \\ &= \nabla_X \left( \mu (Y_1, \dots, Y_n) \right) - \mu ([X, Y_1], Y_2, \dots, Y_n) - \dots - \mu (Y_1, \dots, Y_{n-1}, [X,Y_n]) \\ &= \left(\nabla_X \mu \right) (Y_1, \dots, Y_n) + \mu (\nabla_{Y_1} X, Y_2, \dots, Y_n) + \dots + \mu (Y_1, \dots, Y_{n-1}, \nabla_{Y_n}X) \end{align} For the last equality, one uses the 'product rule' for $\nabla_X$ and its torsion-freeness. Since $\mu$ is parallel, the first summand vanishes. For the other summands, only the component with respect to the $i$th 'missing' $Y_i$ matters. Thus indeed $$f = \sum_{i=1}^n \left( \nabla_{Y_i} X \right)^i = \operatorname{tr} \nabla X \, ,$$ which clearly has to hold everywhere.