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 $$

This does not solve the problem of finding divergence-free fields, but at least avoid globalization problems. You still have to solve a PDE (with more equations but also more degrees of freedom)

**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.


THANKS to Jean Van Schaftingen for pointing out an imprecision in my previous answer.