You can construct   examples a dime a dozen. $\newcommand{\bR}{\mathbb{R}}$

 Here  is a first  simple way. Fix $N$ smooth functions $f_1,f_2,\dotsc, f_N:M\to\bR$ and $N$ independent Gaussian random variables. Then

$$f(x)=\sum_{k=1}^N X_k f_k(x) $$

is a Gaussian random field and the sample functions are a.s. smooth.

You can   allow  infinitely many functions  in the above example, i.e., $N=\infty$, but then  you need to make some assumptions on the functions and on the variances of  the variables $X_k$ to guarantee  the convergence of the resulting series. 

We know from Kolmogorov that the convergence is a $0-1$ phenomenon. The  three-series theorem  tells us what these conditions should be. (In the Gaussian case we can  say a bit more.)

Here is an example of this kind when $M$ is compact. Fix a  Riemann metric, denote by $\Delta$ the resulting Laplacian. Its eigenvalues (multiplicities included) are

$$  0=\lambda_0<\lambda_1\leq \lambda_2\leq \cdots. $$

Fix an orthonormal eigenbasis of $L^2(M, dV_g)$ $(\psi_k)_{k\geq 0}$

$$\Delta \psi_k=\lambda_k,\;\; \Vert\psi_k\Vert_{L^2}=1. $$

Next choose  independent  Gaussian random variables $(X_k)_{k\geq 0}$. We denote by $v_k$ the variance of $X_k$.   If $v_k$ goes to $0$ sufficiently fast, then the random series 

$$f(x)=\sum_{k\geq 0} X_k \psi_k(x) $$

defines a  Gaussian random field on $M$.  The regularity of the sample functions  of this random field   depends on the decay rate of  $v_k$. The faster $v_k$ decays as $k\to\infty$, the more regular is the random function. For example, if 

$$\lim_{k\to \infty}k^\alpha v_k=0,\;\;\forall \alpha>0, $$

then the random function $f(x)$ is a.s. smooth.

Ultimately, the most general construction  of a Gaussian random function on a manifold  is via Gaussian measures on the space of   distributions (i.e. generalized functions) on $M$.  I refer to   [*this paper*][1]  for  more details and additional references. I would start with reference [6] in this paper.


  [1]: http://www3.nd.edu/~lnicolae/GB2.pdf