We consider a class of density-driven flow problems. We are particularly interested in the problem of the salinization of coastal aquifers. We consider the Henry saltwater intrusion problem with uncertain porosity, permeability, and recharge parameters as a test case. The reason for the presence of uncertainties is the lack of knowledge, inaccurate measurements,and inability to measure parameters at each spatial or time location. This problem is nonlinear and time-dependent. The solution is the salt mass fraction, which is uncertain and changes in time. Uncertainties in porosity, permeability, recharge, and mass fraction are modeled using random fields. This work investigates the applicability of the well-known multilevel Monte Carlo (MLMC) method for such problems. The MLMC method can reduce the total computational and storage costs. Moreover, the MLMC method runs multiple scenarios on different spatial and time meshes and then estimates the mean value of the mass fraction. 'The parallelization is performed in both the physical space and stochastic space. To solve every deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion. We use the solution obtained from the quasi-Monte Carlo method as a reference solution.

We present our work [1], Adaptive Multilevel Monte Carlo (AMLMC). It is based upon the work [2], which developed convergence rates for an adaptive algorithm based on isoparametric d-linear quadrilateral finite element approximations and the dual weighted residual error representation in the deterministic setting. The present work can also be seen as an extension of [3]. Our AMLMC algorithm uses as a building block an auxiliary sequence of deterministic, non-uniform meshes, which are generated by the above-mentioned deterministic adaptive algorithm, and they correspond to a geometrically decreasing sequence of tolerances. Specifically, for a given realization of the diffusivity coefficient and a given accuracy level, the AMLMC constructs its approximate sample as the one using the first mesh in the hierarchy that satisfies a corresponding bias accuracy constraint.

Nested integration arises when a nonlinear function is applied to an integrand, and the result is integrated again, which is common in engineering problems, such as optimal experimental design, where typically neither integral has a closed-form expression. Using the Monte Carlo method to approximate both integrals leads to a double-loop Monte Carlo estimator, which is often prohibitively expensive, as the estimation of the outer integral has bias relative to the variance of the inner integrand. For the case where the inner integrand is only approximately given, additional bias is added to the estimation of the outer integral. Variance reduction methods, such as importance sampling, have been used successfully to make computations more affordable. Furthermore, random samples can be replaced with deterministic low-discrepancy sequences, leading to quasi-Monte Carlo techniques. Randomizing the low-discrepancy sequences simplifies the error analysis of the proposed double-loop quasi-Monte Carlo estimator. To our knowledge, no comprehensive error analysis exists yet for truly nested randomized quasi-Monte Carlo estimation (i.e., for estimators with low-discrepancy sequences for both estimations). We derive asymptotic error bounds and a method to obtain the optimal number of samples for both integral approximations. Then, we demonstrate the computational savings of this approach compared to standard nested (i.e., double-loop) Monte Carlo integration when estimating the expected information gain via an example from Bayesian optimal experimental design involving an experiment from solid mechanics.