As Raziel wrote, the local question is whether one can find a local basis of orthonormal vector fields that are divergence-free.

It's true that, in dimension $2$, this can only be done if the metric is locally flat, which is the local obstruction. This is because this is an overdetermined problem; one has two equations for a single unknown.

However, in higher dimensions, it is not clear that there is a local obstruction because this is a system of $n$ first-order PDE for $\tfrac12n(n{-}1)$ unknowns, so, for $n=3$, this is a determined system while, when $n>3$ it is underdetermined.

One can prove that the system for $n=3$ is always locally solvable in the real-analytic case (even though the determined system cannot be written in Cauchy-Kowalevski form, even locally), so there cannot be any local obstruction that is computed on the basis of some kind of curvature condition or identity. (Presumably, it is also always locally solvable in the smooth case as well, but that would require further study.)

*Remark:* It is interesting to note that, if $(M^3,g)$ is real-analytic and possesses a real-analytic orthonormal frame field $X = (X_1,X_2,X_3)$ where each of the $X_i$ are divergence free, then $(M,g)$ can be isometrically embedded as a hypersurface in a Calabi-Yau surface $S$ in such a way that $X_i = I_i(N)$ where $N$ is the oriented unit normal and $I_1$, $I_2$, and $I_3$ are the orthogonal parallel complex structures that define the Calabi-Yau structure.

I expect that, in dimension $n>3$, the problem is so underdetermined that it is always locally solvable, though I have not yet carried out the analysis. *However, see below, where I do complete the analysis in the real-analytic case.*

Global solvability is, of course, much harder, and it's conceivable that there are counterexamples, even for metrics on $S^3$, though I don't know of one.

**The analysis via exterior differential systems:** At the OP's request, I will sketch the EDS analysis of this system. I don't have time to put in all the details, and, in any case, they won't make sense to anyone who doesn't already know Cartan-Kähler theory, but for those who do know this theory, the following explains the proof of the following results:

When $n=2$, local solutions exist if and only if the metric is flat.

When $n=3$, local solutions always exist if the metric is real-analytic. Moreover, the local solutions depend on 2 arbitrary functions of two variables. If the metric is real-analytic and the scalar curvature is positive, then every local solution is also real analytic.

When $n>3$, local solutions always exist if the metric is real-analytic, and the local solutions depend on ${n\choose2}{-}n$ functions of $n$ variables.

In particular, (2) and (3) imply that there are no curvature-type obstructions to local solvability when $n>2$. Whether we have solvability in the smooth case when $n>2$ will require further study.

Here is the argument: Let $(M^n,g)$ be a Riemannian manifold, and let $\pi:F\to M$ be the orthonormal frame bundle, so that an element $f\in F$ is an $n$-tuple $f = (e_1,\ldots,e_n)$ where $e_1,\ldots,e_n$ is an orthonormal basis of $T_xM$ where $x = \pi(f)$. We define the *canonical $1$-forms* $\omega_1,\ldots,\omega_n$ on $F$ so that the equation
$$
\pi'(v) = \omega_1(v)\,e_1 + \cdots \omega_n(v)\,e_n
$$
holds for all $v\in T_fF$ where $f = (e_1,\ldots,e_n)$. A standard result (cf. Kobayashi & Nomizu) then says that there exist unique $1$-forms $\phi_{ij}=-\phi_{ji}$ (where the indices run from $1$ to $n$) satisfying the *first structure equations*
$$
\mathrm{d}\omega_i = -\phi_{ij}\wedge\omega_j\,,
$$
that the forms $\omega_i$ and $\phi_{ij}$ ($i<j$) give a basis for the $1$-forms on $F$,
and that the $\phi_{ij}$ satisfy the *second structure equations*
$$
\mathrm{d}\phi_{ij}
= -\phi_{ik}\wedge\phi_{kj} + \tfrac12\,R_{ijkl}\,\omega_k\wedge\omega_l
$$
for some unique functions $R_{ijkl}=-R_{ijlk}$ on $F$.

A local orthonormal coframing $X = (X_1,\ldots,X_n)$ defined on an open set $U\subset M$
is simply a section of $F$ over $U$, and it satisfies $X^*\omega_i = \xi_i$ where the $\xi_i$ are the $1$-forms on $U$ dual to the $X_i$. The volume form of the metric is, up to a sign, the wedge product of the $\xi_i$, and so the condition that the $X_i$ be divergence free is that
$$
\mathrm{d}\left(\xi_1\wedge\cdots\wedge\widehat{\xi_i}\wedge\cdots\wedge\xi_n\right) = 0
$$
for all $i = 1,\ldots,n$. In other words, defining the $(n{-}1)$-forms
$$
\Omega_i
= (-1)^{i-1}\,\omega_1\wedge\cdots\wedge\widehat{\omega_i}\wedge\cdots\wedge\omega_n\,,
$$
we are requiring that $\mathrm{d}\left(X^*\Omega_i\right)=X^*\left(\mathrm{d}\Omega_i\right) = 0$, i.e., that the image of the section $X$ in $F$ should be an integral manifold of the differential ideal $\mathcal{I}$ on $F$ generated by the $n$ $n$-forms $\mathrm{d}\Omega_i$.

Unfortunately, $\mathcal{I}$ is not involutive. However, it turns out that $\mathcal{I}$ can be enlarged to an ideal $\mathcal{I}_+$ as follows: For $i<j$, define the $(n{-}2)$-forms
$$
\Omega_{ij}
= (-1)^{i+j-1}\,\omega_1\wedge\cdots\wedge\widehat{\omega_i}\wedge
\cdots\wedge\widehat{\omega_j}\wedge\cdots\wedge\omega_n\,,
$$
and set $\Omega_{ii}=0$ while $\Omega_{ji}=-\Omega_{ij}$. Now define the $(n{-}1)$-form
$$
\Upsilon = \phi_{ij}\wedge\Omega_{ij}\,.
$$
It is not hard to show that $\mathrm{d}\Omega_i = \pm\omega_i\wedge\Upsilon$, and, from this, one concludes that $X_i$ is divergence free for all $i$ if and only iff $X^*(\Upsilon)=0$. Thus, we can let $\mathcal{I}_+$ be the differential ideal generated by $\Upsilon$ (i.e., the exterior ideal generated by $\Upsilon$ and $\mathrm{d}\Upsilon$) and look for integral manifolds of this ideal instead.

When $n=2$, $\Upsilon = 2\phi_{12}$, so $\mathrm{d}\Upsilon
= 2R_{1212}\,\omega_1\wedge\omega_2 = 2K\,\mathrm{d}A$, so there are no sections $X$ that are integral manifolds unless $K=0$. (When $K=0$, of course, the sections that are integral manifolds of $\phi_{12}$ are exactly the parallel sections.)

When $n>2$, the structure equations show that $\mathcal{I}_+$, which encodes a system of $n{+}1$ first-order equations on the orthonormal coframing $X$, is involutive, with the Cartan characters of a regular flag being $s_i = 0$ for $i<n{-}2$, $s_{n-2}=1$, $s_{n-1} = n{-}1$, and $s_n = {n\choose2}-n$. Now apply Cartan-Kähler.