  This work describes a domain embedding technique between two non-matching
meshes used for generating realizations of spatially correlated random fields
with applications to large-scale sampling-based uncertainty quantification. The
goal is to apply the multilevel Monte Carlo (MLMC) method for the
quantification of output uncertainties of PDEs with random input coefficients
on general, unstructured computational domains. We propose a highly scalable,
hierarchical sampling method to generate realizations of a Gaussian random
field on a given unstructured mesh by solving a reaction-diffusion PDE with a
stochastic right-hand side. The stochastic PDE is discretized using the mixed
finite element method on an embedded domain with a structured mesh, and then
the solution is projected onto the unstructured mesh. This work describes
implementation details on how to efficiently transfer data from the structured
and unstructured meshes at coarse levels, assuming this can be done efficiently
on the finest level. We investigate the efficiency and parallel scalability of
the technique for the scalable generation of Gaussian random fields in three
dimensions. An application of the MLMC method is presented for quantifying
uncertainties of subsurface flow problems. We demonstrate the scalability of
the sampling method with non-matching mesh embedding, coupled with a parallel
forward model problem solver, for large-scale 3D MLMC simulations with up to
$1.9\cdot 10^9$ unknowns.
