  We develop numerical methods for computing statistics of stochastic processes
on surfaces of general shape with drift-diffusion dynamics $d{X}_t = a({X}_t)dt
+ {b}({X}_t)d{W}_t$. We consider on a surface domain $\Omega$ the statistics
$u(\mathbf{x}) = \mathbb{E}^{\mathbf{x}}\left[\int_0^\tau g(X_t)dt \right] +
\mathbb{E}^{\mathbf{x}}\left[f(X_\tau)\right]$ with the exit stopping time
$\tau = \inf_t \{t > 0 \; |\; X_t \not\in \Omega\}$. Using Dynkin's formula, we
compute statistics by developing high-order Generalized Moving Least Squares
(GMLS) solvers for the associated surface PDE boundary-value problems. We focus
particularly on the mean First Passage Times (FPTs) given by the special case
$f = 0,\, g = 1$ with $u(\mathbf{x}) =
\mathbb{E}^{\mathbf{x}}\left[\tau\right]$. We perform studies for a variety of
shapes showing our methods converge with high-order accuracy both in capturing
the geometry and the surface PDE solutions. We then perform studies showing how
FPTs are influenced by the surface geometry, drift dynamics, and spatially
dependent diffusivities.
