  The numerical approximation of solutions to stochastic partial differential
equations with additive spatial white noise on bounded domains in
$\mathbb{R}^d$ is considered. The differential operator is given by the
fractional power $L^\beta$, $\beta\in(0,1)$, of an integer order elliptic
differential operator $L$ and is therefore non-local. Its inverse $L^{-\beta}$
is represented by a Bochner integral from the Dunford-Taylor functional
calculus. By applying a quadrature formula to this integral representation, the
inverse fractional power operator $L^{-\beta}$ is approximated by a weighted
sum of non-fractional resolvents $(I + t_j^2 L)^{-1}$ at certain quadrature
nodes $t_j>0$. The resolvents are then discretized in space by a standard
finite element method.
  This approach is combined with an approximation of the white noise, which is
based only on the mass matrix of the finite element discretization. In this
way, an efficient numerical algorithm for computing samples of the approximate
solution is obtained. For the resulting approximation, the strong mean-square
error is analyzed and an explicit rate of convergence is derived. Numerical
experiments for $L=\kappa^2-\Delta$, $\kappa > 0$, with homogeneous Dirichlet
boundary conditions on the unit cube $(0,1)^d$ in $d=1,2,3$ spatial dimensions
for varying $\beta\in(0,1)$ attest the theoretical results.
