Why is the covariant derivative of the determinant of the metric zero?
Comments to the post (v2):
Note that $\sqrt{|g|}$ transforms as a density rather than a scalar under general coordinate transformations. In particular, the covariant derivative of $\sqrt{|g|}$ does not necessarily coincide with the partial derivative of $\sqrt{|g|}$.
Here is a heuristic explanation using local coordinates. The Levi-Civita connection is compatible with the metric $g_{\mu\nu}\mathrm{d}x^{\mu}\odot \mathrm{d}x^{\nu}$. That a connection $\nabla$ is compatible with a metric means that $\nabla_{\lambda}g_{\mu\nu}=0$. Using linearity and Leibniz rule, the covariant derivative $\nabla_{\lambda}$ then annihilates any sufficiently nice function $f(g_{00},g_{01}, \ldots)$ of the metric. In particular, the square root of the determinant $\sqrt{|g|}$, so $\nabla_{\lambda}\sqrt{|g|}=0$.
OK. Let us take ordinary derivative of determinant of some covariant 2-tensor $A_{\mu\nu}$. Let call it $A$. But it is more convenient to allow us to think about $A_{\mu\nu}$ like a matrix with covariant indices. So $$\det{A_{\mu\nu}} = A$$ Next, let's do the following calculations: $$\delta\ln{\det{A_{\mu\nu}}} = \ln{\det{(A_{\mu\nu}}+\delta A_{\mu\nu})}-\ln{\det{A_{\mu\nu}}} = \ln{\det({A^{\mu\sigma}(A_{\sigma\nu}+\delta A_{\sigma\nu}))}},$$ where $\delta$ is like differential and $A^{\mu\sigma}$ denotes contravriant 2-tensor with the following property: $A^{\mu\sigma}A_{\sigma\nu} = \delta^{\mu}_{\,\nu}$, in other words, "inverse" tensor.
Let's continue $$\ln{\det({A^{\mu\sigma}(A_{\sigma\nu}+\delta A_{\sigma\nu}))}} = \ln{\det{(I+A^{\mu\sigma}\delta A_{\sigma\nu})}} = \ln{(1 + \mathrm{Tr}\,{A^{\mu\sigma}\delta A_{\sigma\nu}})} = \mathrm{Tr}\,{A^{\mu\sigma}\delta A_{\sigma\nu}}.$$ But $$\mathrm{Tr}\,{A^{\mu\sigma}\delta A_{\sigma\nu}} = A^{\mu\sigma}\delta A_{\sigma\mu}$$
Divided by $dx^{\lambda}$ it gives $$\partial_{\lambda}\ln{\det{A_{\mu\nu}}} = A^{\mu\sigma}\partial_{\lambda} A_{\sigma\mu}.$$
Therefore, $$\frac{\partial_{\lambda}A}{A} = A^{\mu\sigma}\partial_{\lambda} A_{\sigma\mu}$$
Or $$\partial_{\lambda}g = g g^{\mu\sigma}\partial_{\lambda} g_{\sigma\mu}$$ The following step is pretty fun. Let's replace all ordinary partials by absolute (covariant). So we have $$\nabla_{\lambda}g = g g^{\mu\sigma}\nabla_{\lambda} g_{\sigma\mu}.$$ But $$\nabla_{\lambda} g_{\sigma\mu} = 0.$$ QED. The last is not hard exercise. Indeed, in the geodesic coordinates it is always true because in these coordinates $\nabla_{\nu} = \partial_{\nu}$. But if some tensor is equal to zero in one reference frame then it is zero in every reference frame.
But I am not sure about the same trick with arbitrary matrix (although it may turn out the same). It would be better to use the following. Since $$\det{g^{\mu\nu}A_{\nu\sigma}}$$ is a scalar, we can use ordinary derivative for this. But on the other hand, we could use covariant derivative for it. For scalar it is the same. So $$\nabla_{\nu}(\det{g^{\mu\nu}A_{\mu\nu}}) = g^{-1}\nabla_{\nu}A + A\nabla_{\nu}g^{-1} = g^{-1}\partial_\nu A + A\partial_\nu g^{-1}$$ Let us continue calculations $$\nabla_{\nu}A = \partial_{\nu} A - A\frac{\partial_\nu g}{g}$$ Where we used $\nabla_\nu g = 0$. Partial derivatives we can find from the previous equations.
Here's a heuristic calculation: Let $\{E_i\}$ be an orthonormal frame ($g(E_i,E_j)=\epsilon_i\delta_{ij}, \epsilon_i=\pm 1$). Then $\mu$ is the canonical volume form $\sqrt{g}\,\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n$ iff $\mu(E_1,\dotsc, E_n)=1$. Then $$(\nabla_X\mu)(E_1,\dotsc,E_n)=\nabla_X(\mu(E_1,\dotsc,E_n))-\sum \mu(E_1,\dotsc,\nabla_X E_i,\dotsc,E_n)=-\sum \epsilon_ig(E_i,\nabla_X E_i)\mu(E_1,\dotsc,E_i,\dotsc,E_n)=-\sum \epsilon_ig(E_i,\nabla_X E_i)=-\frac{1}{2}\sum \epsilon_i\nabla_X g(E_i,E_i)=0$$ for all vector fields $X$. Then using the derivation property of the connection, and $\nabla_X(\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n)=0$ for all $X$, one has $$\nabla_X\mu=(\nabla_X\sqrt{g})\,\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n=0$$ whence $\nabla_X\sqrt{g}=0$ for all $X$ and in some chart.
To be absolutely pedantic, one should adapt the defintion of the connection in terms of parallel transport to tensor densities. This is done in e.g. Straumann, General Relativity (2013). For a scalar density $\rho$ one finds in local coordinates $\nabla_i\rho=(\partial_i-\Gamma^l{}_{il})\rho$. From the standard expression for $\Gamma^l{}_{li}$ it is easy to verify that $\nabla_i\sqrt{g}=0$.