  Multilevel Monte Carlo is a key tool for approximating integrals involving
expensive scientific models. The idea is to use approximations of the integrand
to construct an estimator with improved accuracy over classical Monte Carlo. We
propose to further enhance multilevel Monte Carlo through Bayesian surrogate
models of the integrand, focusing on Gaussian process models and the associated
Bayesian quadrature estimators. We show, using both theory and numerical
experiments, that our approach can lead to significant improvements in accuracy
when the integrand is expensive and smooth, and when the dimensionality is
small or moderate. We conclude the paper with a case study illustrating the
potential impact of our method in landslide-generated tsunami modelling, where
the cost of each integrand evaluation is typically too large for operational
settings.
