This question, metric determinant and its partial and covariant derivative, seems to indicate $$\nabla_a \sqrt{g}=0.$$ Why is this the case? I've always learned that $$\nabla_a f= \partial_a f,$$ hence surely $$\nabla_a \sqrt{g}= \partial_a\sqrt{g} \neq 0. $$
Where's the hole in my logic?
With $g:=\det(g_{\mu\nu})$ note that $\sqrt{|g|}$ transforms as a density $$\begin{align}\rho^{(x)}\quad\longrightarrow&\quad \rho^{(y)}~=~\frac{\rho^{(x)}}{|J|}, \cr J~:=~&\det\frac{\partial y^{\nu}}{\partial x^{\mu}}, \end{align}$$ rather than a scalar $$s^{(x)}\quad\longrightarrow\quad s^{(y)}~=~s^{(x)},$$ under general coordinate transformations $$x^{\mu}\quad\longrightarrow\quad y^{\nu}~=~f^{\nu}(x).$$ In particular, the covariant derivative of $\sqrt{|g|}$ does not necessarily coincide with the partial derivative of $\sqrt{|g|}$.
It is consistent to define that a connection $\nabla$ acts on a density $\rho>0$ as $$ \nabla_{\mu}\ln\rho~:=~\partial_{\mu}\ln\rho-\Gamma^{\nu}_{\mu\nu},$$ because it then transforms covariantly $$ \frac{\partial y^{\nu}}{\partial x^{\mu}}\nabla^{(y)}_{\nu}\ln\rho^{(y)}~=~\nabla^{(x)}_{\mu}\ln\rho^{(x)}.$$ A connection $\nabla$ and a density $\rho$ are by definition compatible if $$\nabla_{\mu}\ln\rho~=~0\quad\Leftrightarrow\quad\partial_{\mu}\ln\rho~=~\Gamma^{\nu}_{\mu\nu} .$$ Using linearity and Leibniz rule, we conclude that $$\nabla_{\mu}f(\ln\rho) ~=~0$$ for a sufficiently nice (e.g. real analytic) function $f$.
A connection $\nabla$ is by definition compatible with a metric $$\mathbb{\bf g}~=~g_{\mu\nu}\mathrm{d}x^{\mu}\odot \mathrm{d}x^{\nu}$$ if $$\nabla_{\lambda}\mathbb{\bf g} ~=~0 \quad\Leftrightarrow\quad (\nabla_{\lambda}\mathbb{\bf g})_{\mu\nu}~=~0.$$ E.g. the Levi-Civita connection is compatible with the metric $\mathbb{\bf g}$. Metricity implies that$^1$ $$\begin{align} 2\Gamma^{\kappa}_{\lambda\kappa} ~=~&g^{\mu\nu}\partial_{\lambda}g_{\nu\mu}\cr ~=~&g^{\mu\nu}\eta_{\nu\sigma}\partial_{\lambda}\eta^{\sigma\tau}g_{\tau\mu}\cr ~=~&\partial_{\lambda}{\rm tr}{\rm Ln}(\eta^{\sigma\tau}g_{\tau\mu})\cr ~=~&\partial_{\lambda}\ln\det(\eta^{\sigma\tau}g_{\tau\mu})\cr ~=~&\partial_{\lambda}\ln|g|,\end{align}$$ which in turn shows that the density $\rho=\sqrt{|g|}$ is compatible with a metric connection $\nabla$. This proves OP's title question. Note that it is unnecessary to assume that the connection $\nabla$ is torsionfree. $\Box$
Below is an alternative approach that doesn't use the transformation rule for densities: A metric connection satisfies $$\nabla_{\lambda}(\mathbb{\bf g}^{\otimes n}) ~=~0\quad\Leftrightarrow\quad (\nabla_{\lambda}(\mathbb{\bf g}^{\otimes n}))_{\mu_1\nu_1\ldots\mu_n\nu_n}~=~0.$$ Using linearity and Leibniz rule, the covariant derivative $\nabla_{\lambda}$ then annihilates any sufficiently nice (e.g. real analytic) function $f(g_{00},g_{01}, \ldots)$ of the metric. In particular, the square root of the determinant $\sqrt{|g|}$, so that $$\nabla_{\lambda}\sqrt{|g|}~=~0.$$ $\Box$
--
$^1$ The above metric $\eta_{\mu\nu}$ is merely a trick to deal with the absolute value $|g|$ of $g$. The metric $\eta_{\mu\nu}$ is assumed to
(i) have the same signature as $g_{\mu\nu}$,
(ii) have $\det(\eta_{\mu\nu})=\pm 1$, and
(iii) have constant components in the pertinent local coordinate neighborhood.
There is an analogous trick to deal with the absolute value $|J|$ of the Jacobian determinant $J$.
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$.
Actually, the very notation $\nabla_\alpha \sqrt{g}$ is wrong, because $g = \det g_{\mu\nu}$ is coordinate dependent, and thus does not form a scalar field.
To make it more concrete, consider a 2-dimensional Euclidean space with polar coordinates $(r, \theta)$. $g = \mathrm{diag} (1, r^2)$. So $\sqrt{g} = r$. If $\nabla_r \sqrt{g}$ were meaningful, we would have $\nabla_r \sqrt{g} = \partial_r r = 1$, contradicting the claim in your original question.
I can only guess that you asked this question when you were trying to calculate derivatives of tensor fields such as $\epsilon_{\mu\nu\rho\sigma}$, in which $\sqrt{g}$ is used to turn a tensor density into a proper tensor. My answer to another question might be helpful.
