  We analyze several Galerkin approximations of a Gaussian random field
$\mathcal{Z}\colon\mathcal{D}\times\Omega\to\mathbb{R}$ indexed by a Euclidean
domain $\mathcal{D}\subset\mathbb{R}^d$ whose covariance structure is
determined by a negative fractional power $L^{-2\beta}$ of a second-order
elliptic differential operator $L:= -\nabla\cdot(A\nabla) + \kappa^2$. Under
minimal assumptions on the domain $\mathcal{D}$, the coefficients
$A\colon\mathcal{D}\to\mathbb{R}^{d\times d}$,
$\kappa\colon\mathcal{D}\to\mathbb{R}$, and the fractional exponent $\beta>0$,
we prove convergence in $L_q(\Omega; H^\sigma(\mathcal{D}))$ and in
$L_q(\Omega; C^\delta(\overline{\mathcal{D}}))$ at (essentially) optimal rates
for (i) spectral Galerkin methods and (ii) finite element approximations.
Specifically, our analysis is solely based on
$H^{1+\alpha}(\mathcal{D})$-regularity of the differential operator $L$, where
$0<\alpha\leq 1$. For this setting, we furthermore provide rigorous estimates
for the error in the covariance function of these approximations in
$L_{\infty}(\mathcal{D}\times\mathcal{D})$ and in the mixed Sobolev space
$H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})$, showing convergence which is
more than twice as fast compared to the corresponding $L_q(\Omega;
H^\sigma(\mathcal{D}))$-rate. For the well-known example of such Gaussian
random fields, the original Whittle-Mat\'ern class, where $L=-\Delta +
\kappa^2$ and $\kappa \equiv \operatorname{const.}$, we perform several
numerical experiments which validate our theoretical results.
