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
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$.
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 is a way to possibly avoid globalization problems. You still have to solve a PDE. 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 for pointing out an imprecision in my previous answer.