Can the Laplace operator on $n-$ manifolds be represented as a sum of $n$ second order derivational operators
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.
You can use the general local formula for the Laplace-Beltrami operator in terms of any local orthonormal frame:
$$\Delta = \sum_{i=1}^n W_i^2 +\mathrm{div}(W_i)W_i$$
where the $W_i$'s are seen as derivations on functions.
You can always find a local frame of vector fields $W_1,\ldots,W_n$ that are divergence-free at a given point $q$. In terms of this frame, the Laplacian at the point $q$ is just a "sum of squares".
Locally, the construction of a local divergence-free, orthonormal frame leads to a system of first order PDEs. The integrability conditions then give a local obstruction.
UPDATE WITH $N \geq \dim M$ FIELDS
Unless your manifold is parallelizable you can't find a global orthonormal frame. Still, the above formula works also if the number of vector fields $W_1,\ldots,W_N$ is greater than the dimension of the manifold $n=\dim M \leq N$. To see this practically, pick a orthonormal frame $X_1,\ldots,X_n$ (local on $U \subset M$). We have
$$ W_I = \sum_{j=1}^n A_{Ij} X_j, \qquad I=1,\ldots,N $$
for some smooth family of $N \times n$ matrix $A : U \to M_{N\times n}$. Assume that
$$A^T A = \mathbb{I}_n$$
on $U$. Then you can check that for any function $f \in C^\infty(M)$
$$ \sum_{I=1}^N W_I(f) W_I = \sum_{I=1}^N \sum_{i,j=1}^n A_{Ij}A_{Ii} X_j(f) X_i(f) = \sum_{i=1}^n X_i(f) X_i = \nabla f$$
where $\nabla f $ is the Riemannian gradient of $f$. Then
$$\Delta f = \mathrm{div}(\nabla(f)) = \sum_{I=1}^N W_I^2(f) + \mathrm{div}(W_I)W_I(f) \tag{1}$$
That is the formula at the beginning of this answer. Observe that all of this (starting from the definition of the matrix $A$) is local, since the $X_i$'s are local, but clearly formula (1) holds wherever the $W_i$'s are defined (i.e. globally).
More abstractly, the initial formula holds true for any set of vector fields $W_1,\ldots,W_N$ (local or global) such that the symbol (as a function on $T^*M$) is written
$$ \lambda \mapsto \sum_{I=1}^N \langle\lambda, W_I\rangle^2, \qquad \lambda \in T^*M,$$
where $\langle \lambda, \cdot\rangle$ denotes the action of covectors on vectors. Equivalently, any set of vector fields $W_1,\ldots,W_N$ (local or global) such that
$$ \|Z\|^2 = \sum_{I=1}^N g(Z,W_I)^2, \qquad Z \in \Gamma(TM) $$
This indeed puts constraints on your $W_I$ as, for example, $\|W_I\| \leq 1$ and at least one (actually 2) of them will have $\|W_I\| < 1$ as soon as $N > n$.
EXPLICIT EXAMPLE on $\mathbb{S}^2$
As an example, on the $2$-sphere $\mathbb{S}^2 \subset \mathbb{R}^3$, take three global vector fields $W_1,W_2,W_3$ obtained by taking the orthogonal projection of the fields $\partial_x,\partial_y,\partial_z$ of $\mathbb{R}^3$ on the sphere. If you work out the details you obtain (in spherical coordinates)
\begin{eqnarray*} W_1 &=& \cos\theta\cos\phi \partial_\theta - \frac{\sin\phi}{\sin\theta}\partial_\phi \\ W_2 &=& \cos\theta\sin\phi \partial_\theta + \frac{\cos\phi}{\sin\theta}\partial_\phi \\ W_3 &=& -\sin\theta\partial_\theta \end{eqnarray*}
and you can check that the standard spherical Laplacian on $\mathbb{S}^2$ is
$$\Delta_{\mathbb{S}^2} = W_1^2+W_2^2+W_3^2$$
In particular the ''divergence part'' is zero with this particular construction.
This construction indeed works for any manifold by taking an isometric embedding on an $R^N$ of sufficiently large dimension.
BACK TO THE ORIGINAL QUESTION
This does not solve the problem of finding divergence-free fields, but at least is a way to possibly avoid globalization problems. You still have to solve a PDE, to kill the first order part $\sum_{i=1}^N \mathrm{W}_I X_I$. A naive parameter counting shows that you have $n$ equations with $\frac{n(2N-n-1)}{2}$ degrees of freedom, so, unless there is some hidden extra constraint the problem seems easier when $N > n$.
THANKS to Jean Van Schaftingen and Robert Bryant for pointing out an imprecision in my previous answer.
The point of this answer is to flesh out my comment above: Any such $X_i$ must be everywhere linearly independent, so they only exist if the manifold $M$ is parallelizable. Suppose, for the sake of contradiction, that the $X_i$ become linearly dependent at some point $p$. Pass to local coordinates $(x_1, \ldots, x_n)$ with $p=(0,0,\ldots, 0)$ and write $$X_i = \sum_j a_{ij}(x) \frac{\partial}{\partial x_j}.$$ Let $A(x)$ denote the matrix $(A_{ij}(x))$. Our hypothesis is that $A(0)$ is not of full rank.
We compute the principal symbol of $\sum (\partial/\partial X_i)^2$. We have $$\left( \sum_j a_{ij}(x) \frac{\partial}{\partial x_j} \right)^2 = \sum_{j_1, j_2} a_{i j_1}(x) a_{i j_2}(x) \frac{\partial^2}{(\partial x_{j_1}) (\partial x_{j_2})} + \mbox{first order operators}.$$ So the principal symbol is $$\sum_i \sum_{j_1, j_2} a_{i j_1}(x) a_{i j_2}(x) y_{j_1} y_{j_2} = y^T A^T(x) A(x) y$$ where $y$ is the vector $(y_1, y_2, \ldots, y_n)^T$. If $A$ is not of full rank, neither is this quadratic form.
But the symbol of $\Delta$ is positive definite, a contradiction.