  We propose an alternative method to generate samples of a spatially
correlated random field with applications to large-scale problems for forward
propagation of uncertainty. A classical approach for generating these samples
is the Karhunen-Lo\`{e}ve (KL) decomposition. However, the KL expansion
requires solving a dense eigenvalue problem and is therefore computationally
infeasible for large-scale problems. Sampling methods based on stochastic
partial differential equations provide a highly scalable way to sample Gaussian
fields, but the resulting parametrization is mesh dependent. We propose a
multilevel decomposition of the stochastic field to allow for scalable,
hierarchical sampling based on solving a mixed finite element formulation of a
stochastic reaction-diffusion equation with a random, white noise source
function. Numerical experiments are presented to demonstrate the scalability of
the sampling method as well as numerical results of multilevel Monte Carlo
simulations for a subsurface porous media flow application using the proposed
sampling method.
