This holds in complete generality. I sketch the case in which, instead of $B$ you consider the metric sphere. The case with the ball can be obtained similarly.

In fact, choose normal coordinates in a neighbhorhood of $x$, such that

$$u(x+y)-u(x) = \sum_{i=1}^k a_i y_i + \frac{1}{2} b_{ij} 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. 

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.

Hence, in your point $x$, in normal coordinates, you have

$$ \lim_{h\to 0} \frac{1}{h^2}(u(x) - \frac{1}{|S(x,h)|}\int_{S(x,h)} u(y) dy \propto \sum_{i=1}^n (\partial_i^2 u)(x)$$

At the point $x$ this coincides with $-\Delta_g u$ written in normal coordinates centered in $x$. Since your averaging expression is coordinate invariant, the equality holds in general.