  We develop numerical methods for computing statistics of stochastic processes
on surfaces of general shape with drift-diffusion dynamics $d\mathbf{X}_t =
a(\mathbf{X}_t)dt + \mathbf{b}(\mathbf{X}_t)d\mathbf{W}_t$. We formulate
descriptions of Brownian motion and general drift-diffusion processes on
surfaces. We consider statistics of the form $u(\mathbf{x}) =
\mathbb{E}^{\mathbf{x}}\left[\int_0^\tau g(\mathbf{X}_t)dt \right] +
\mathbb{E}^{\mathbf{x}}\left[f(\mathbf{X}_\tau)\right]$ for a domain $\Omega$
and the exit stopping time $\tau = \inf_t \{t > 0 \; |\; \mathbf{X}_t \not\in
\Omega\}$, where $f,g$ are general smooth functions. For computing these
statistics, we develop high-order Generalized Moving Least Squares (GMLS)
solvers for associated surface PDE boundary-value problems based on
Backward-Kolmogorov equations. We focus particularly on the mean First Passage
Times (FPTs) given by the case $f = 0,\, g = 1$ where $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
statistics are influenced by the surface geometry, drift dynamics, and
spatially dependent diffusivities.
