Does the covariant derivative still give divergence in skew coordinates? I stumbled upon the formula $\,\,div \, \vec{F}=F^{\mu}_{\,\,\,;\mu}$. Does this still hold true in skew coordinates? I can picture it working geometrically in orthogonal coordinates, but in skew coordinates I'm having trouble. I'm imagining a small parallelepiped with edges parallel to the coordinate axes and considering the changes in flux through the sides. With orthogonal coordinates, we only have to consider the change in the $x^{\mu}$ component of the vector field in the $x^{\mu}$ direction--any "transverse" changes to the vector don't alter the flux through that face. But in skew coordinates that would seem not to be so.
1) Does the formula still hold in skew coordinates?
2) If so, is there a good geometric picture as to why?
 A: Since the "covariant" derivative has been constructed precisely to be coordinate-independent, yes, the formula $\text{div}X=X^\mu_{\ ;\mu}$ is of course valid.
Note that if the covariant derivative arises from the Levi-Civita connection (eg. it is torsionless and metric-compatible), then it can also be expressed as $$ \text{div}X=\frac{1}{\sqrt{|g|}}\partial_\mu(X^\mu\sqrt{|g|}), $$ where $g$ denotes the determinant of the metric tensor's matrix representation.

On the other hand, you seem to be more concerned with whether the usual "visual" interpretation of the divergence is still valid. The answer is yes, but to satisfactorially see this, one needs to adopt a different point of view. If you do not know what a Lie derivative or a differential form is, you might want to look them up, though I'll try to avoid being overly technical.
Assume that the manifold is equipped with a preferred volume form $\Omega$, which, in a local chart is expressed as $$ \Omega=\rho(x)dx^1\wedge...\wedge dx^n. $$ If you do not know what this is, but know what a scalar density (of weight 1) is, you can imagine that the "component" $\rho$ is a nowhere-vanishing scalar density.
If $X$ is a vector field, we can take the Lie derivative $$ \mathcal L_X\Omega, $$ which, by the properties of the Lie derivative, is the same type of object as $\Omega$ is (alternatively, from the coordinate-based point of view, you can take the Lie derivative of the scalar density). Because $n$-forms/scalar densities have only one component, any $n$-form/scalar density can be expressed as a scalar field times $\Omega$ (or $\rho$), so we have $$ \mathcal L_X\Omega=\text{div}(X)\cdot\Omega, $$ where the factor of proportionality is called the divergence of $X$.
It can be shown that the divergence can be calculated as $$ \text{div}X=\frac{1}{\rho}\partial_\mu(X^\mu\rho) $$ in any coordinate system.
It is clear that the divergence operator is associated to a volume measure/form/density. However if you are given a metric, then you also have a volume form $$ đ\mu=\sqrt{|g|}dx^1\wedge...\wedge dx^n, $$ and then the formula for divergence is the same as the formula for the divergence when it is calculated from the Levi-Civita connection.

This aside is useful for interpretation, because the Lie derivative $$ \mathcal L_X\Omega $$ expresses how the volume measure changes under the flow of the vector field $X$. So if this is zero, $X$ is volume-preserving, and if it is nonzero, then $X$ is not volume-preserving.
So, the divergence measures how much a vector field fails to be volume-preserving, which is very much like the usual "flux" interpretation in vector calculus.
Note: This interpretation works when you take the divergence with respect to a volume form, or with respect to the Levi-Civita connection (and so, with respect to the Riemannian volume-form). This interpretation however fails if $D_\mu$ is a completely arbitrary connection. I do not know of any "pretty" interpretation of the divergence then.
As far as I can remember, if the connection $D_\mu$ is $\mathfrak{sl}(n,\mathbb R)$-valued (happens if and only if the curvature tensor of $D$ is traceless, as in $R^\sigma_{\ \sigma\mu\nu}=0$), then there exists a preferred volume form, which is preserved by $D_\mu$. I think (but am not totally sure) that in that case, the divergence with respect to $D_\mu$ agrees with the divergence with respect to this preferred volume form.
