You can use the general local formula for the Laplace-Beltrami operator in terms of any local orthonormal frame:
$$\Delta = \sum_{i=1}^n X_i^2 +\mathrm{div}(X_i)X_i$$
where the $X_i$'s are seen as derivations on functions.
You can always find a local frame of vector fields $X_1,\ldots,X_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
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_i $$
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$
$$ \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$ is the Riemannian gradient. 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 my 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) is true 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 \lambda(W_I)^2, \qquad \lambda \in T^*M$$
Equivalently, any set of vector fields $W_1,\ldots,W_N$ (local or global) such that
$$ \|Z\|^2 = \sum_{I=1}^N g(X,W_I)^2 $$
EXAMPLE
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. 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 avoid globalization problems. You still have to solve a PDE, that now has $N$ equations, and $\frac{n(2N-n-1)}{2}$ degrees of freedom.
THANKS to Jean Van Schaftingen for pointing out an imprecision in my previous answer.