We analyze several types of Galerkin approximations of a Gaussian random field Z: D× Ω→ R indexed by a Euclidean domain D⊂ Rd whose covariance structure is determined by a negative fractional power L-2β of a second-order elliptic differential operator L: = - ∇ · (A∇) + κ2. Under minimal assumptions on the domain D, the coefficients A: D→ Rd×d, κ: D→ R, and the fractional exponent β> 0 , we prove convergence in Lq(Ω; Hσ(D)) and in Lq(Ω; Cδ(D¯)) at (essentially) optimal rates for (1) spectral Galerkin methods and (2) finite element approximations. Specifically, our analysis is solely based on H1+α(D) -regularity of the differential operator L, where 0 < α≤ 1. For this setting, we furthermore provide rigorous estimates for the error in the covariance function of these approximations in L∞(D× D) and in the mixed Sobolev space Hσ,σ(D× D) , showing convergence which is more than twice as fast compared to the corresponding Lq(Ω; Hσ(D)) -rate. We perform several numerical experiments which validate our theoretical results for (a) the original Whittle–Matérn class, where A≡IdRd and κ≡ const. , and (b) an example of anisotropic, non-stationary Gaussian random fields in d= 2 dimensions, where A: D→ R2 × 2 and κ: D→ R are spatially varying.