  We interpret uncertainty in a model for seismic wave propagation by treating
the model parameters as random variables, and apply the Multilevel Monte Carlo
(MLMC) method to reduce the cost of approximating expected values of selected,
physically relevant, quantities of interest (QoI) with respect to the random
variables. Targeting source inversion problems, where the source of an
earthquake is inferred from ground motion recordings on the Earth's surface, we
consider two QoI that measure the discrepancies between computed seismic
signals and given reference signals: one QoI, $\hbox{QoI}_E$, is defined in
terms of the $L^2$-misfit, which is directly related to maximum likelihood
estimates of the source parameters; the other, $\hbox{QoI}_W$, is based on the
quadratic Wasserstein distance between probability distributions, and
represents one possible choice in a class of such misfit functions that have
become increasingly popular to solve seismic inversion in recent years. We
simulate seismic wave propagation, including seismic attenuation, using a
publicly available code in widespread use, based on the spectral element
method. Using random coefficients and deterministic initial and boundary data,
we present benchmark numerical experiments with synthetic data in a
two-dimensional physical domain and a one-dimensional velocity model where the
assumed parameter uncertainty is motivated by realistic Earth models. Here, the
computational cost of the standard Monte Carlo method was reduced by up to 97%
for $\hbox{QoI}_E$, and up to 78% for $\hbox{QoI}_W$, using a relevant range of
tolerances. Shifting to three-dimensional domains is straight-forward and will
further increase the relative computational work reduction.
