Since we are interested in a local result, we may
assume without loss of generality that \(M\)  is an open disk
in \(\mathbb {R}^2\) , centered at \(p\in \mathbb {R}^2\) , with a Riemannian metric \(g\) .
By Darboux' theorem [1]}, there exists a local
coordinate system \(y = (y^1, \, y^2)\)  such that
\( \mathrm {vol}_g= \mathrm {d}y^1\wedge \mathrm {d}y^2\) 
Let \(A = (A_k^j)_{j,k=1}^2\in \mathbb {R}^{2\times 2}\)  be the unique
(constant) matrix such that
\( \frac{\partial }{\partial y^j}_{|p} A_k^j = \hat{e}_k\qquad \textrm {for } k = 1, \, 2.\) 
The matrix \(A\)  is invertible, because \(\lbrace \hat{e}_1, \, \hat{e}_2\rbrace \) 
is a basis for \(T_pM\) . We define local coordinates \((x^1, \, x^2)\) 
in such a way that
\( y^j = A_k^j x^k \qquad \textrm {for } j = 1, \, 2.\) 
Then, the chain rule and (REF ) imply
\( \frac{\partial }{\partial x^k}_{|p}= \frac{\partial }{\partial y^j}_{|p} A_k^j = \hat{e}_k\qquad \textrm {for } k = 1, \, 2.\) 
It only remains to show that the coordinate
fields \(\frac{\partial }{\partial x^k}\)  are divergence-free,
for \(k=1, \, 2\) .
Thanks to (REF ),
(REF ) and (REF ), we obtain
\((\mathrm {div}\,X)\mathrm {vol}_g= \mathrm {d}\!\left(\mathbf {i}_{X}(\mathrm {d}y^1\wedge \mathrm {d}y^2)\right)= (\det A) \, \mathrm {d}\!\left(\mathbf {i}_{X}(\mathrm {d}x^1\wedge \mathrm {d}x^2)\right)\) 
for any smooth vector field \(X\) .
Taking \(X = \frac{\partial }{\partial x^1}\) , we obtain
\(\begin{split}\left(\mathrm {div}\,\frac{\partial }{\partial x^1}\right)\mathrm {vol}_g= (\det A) \, \mathrm {d}\!\left(\mathbf {i}_{\frac{\partial }{\partial x^1}}\mathrm {d}x^1\wedge \mathrm {d}x^2\right) = (\det A) \, \mathrm {d}(\mathrm {d}x^2) = 0\end{split}\) 
In a similar way, \(\mathrm {div}\,(\frac{\partial }{\partial x^2}) = 0\) .
