This holds in complete generality. I sketch the case in which, instead of the ball $B(x,h)$, you consider the metric sphere $S(x,h)$. In particular (here $\dim M = n$ and the Laplacian is $\Delta = \mathrm{div}\circ \mathrm{grad}$):
$$ (\Delta u)(x) = \lim_{h\to 0} \frac{2n}{h^2}\frac{1}{|S(x,h)|}\int_{S(x,h)} [u(y)-u(x)] dy .$$
The case of the ball is similar with a different constant. Choose normal coordinates in a neighbhorhood of $x$, such that
$$u(y)-u(x) = \sum_{i=1}^k (\partial_i u)(x) y_i + \frac{1}{2} \sum_{i,j=1}^n (\partial_{ij}^2 u)(x) y_i y_j + O(|y|^3). $$
The Riemannian metric in these coordinates is $g_{ij}= \delta_{ij}+O(|y|^2)$, hence up to higher orders the measure is the Euclidean one. Also the metric sphere, in these coordinates, for small $h$, coincides with the Euclidean sphere $\mathbb{S}^{n-1}$.
When you take the average over the sphere of each term, the linear part averages to zero (as the integral of a linear function over the sphere). The integral of the quadratic part averages to the sum of second derivatives. To see this, use the fact that for any symmetric matrix $M$
$$\int_{\mathbb{S}^{n-1}} x^* Q x \,d\mu_{\mathbb{S}^{n-1}} = \frac{|\mathbb{S}^{n-1}|}{n} \mathrm{Tr}(Q).$$
Hence, at your point $x$, in normal coordinates, you have
$$ \lim_{h\to 0} \frac{2n}{h^2}\frac{1}{|S(x,h)|}\int_{S(x,h)} [u(y)-u(x)] dy = \sum_{i=1}^n (\partial_i^2 u)(x).$$
At the point $x$ this coincides with $(\Delta u)(x)$ written in normal coordinates centered in $x$. Since your averaging expression is coordinate invariant, the equality holds in general.
