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 $X_1,\dotsc, X_N$. 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 $g$ on $M$, 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 $(\psi_k)_{k\geq 0}$ of $L^2(M, dV_g)$
$$\Delta \psi_k=\lambda_k,\;\; \Vert\psi_k\Vert_{L^2}=1,\;\;\forall k\geq 0. $$
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.