  Spherical Whittle--Mat\'ern Gaussian random fields are considered as
solutions to fractional elliptic stochastic partial differential equations on
the sphere. Approximation is done with surface finite elements. While the
non-fractional part of the operator is solved by a recursive scheme, a
quadrature of the Dunford--Taylor integral representation is employed for the
fractional part. Strong error analysis is performed, and the computational
complexity is bounded in terms of the accuracy. Numerical experiments for
different choices of parameters confirm the theoretical findings.
