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 examples, but then you need to make some assumptions on the functions and the variables $X$ to guarantee the convergence of the resulting series.
We know from Kolmogorov that the convergence is a $0-1$ phenomenon. The the three-series theorem tells us what these conditions should be.
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 for more details and additional references. I would start with reference [6] in this paper.