Uncertainty quantification in the Henry problem using the multilevel Monte Carlo method

We investigate the applicability of the well-known multilevel Monte Carlo (MLMC) method to the class of density-driven flow problems, in particular the problem of salinisation of coastal aquifers. As a test case, we solve the uncertain Henry saltwater intrusion problem. Unknown porosity, permeability and recharge parameters are modelled by using random fields. The classical deterministic Henry problem is non-linear and time-dependent, and can easily take several hours of computing time. Uncertain settings require the solution of multiple realisations of the deterministic problem, and the total computational cost increases drastically. Instead of computing of hundreds random realisations, typically the mean value and the variance are computed. The standard methods such as the Monte Carlo or surrogate-based methods is a good choice, but they compute all stochastic realisations on the same, often, very fine mesh. They also do not balance the stochastic and discretisation errors. These facts motivated us to apply the MLMC method. We demonstrate that by solving the Henry problem on multi-level spatial and temporal meshes, the MLMC method reduces the overall computational and storage costs. To reduce the computing cost further, parallelization is performed in both physical and stochastic spaces. To solve each deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion.


Introduction
We investigate the applicability of the well-known multilevel Monte Carlo (MLMC) method to density-driven flow problems, in particular the problem of salinisation of coastal aquifers.As a numerical test case, we consider the well-known Henry problem in a stochastic setting.The novelty here is that parameters such as porosity, permeability and recharge are not deterministic but stochastic (uncertain).The reason for the presence of uncertainties is the lack of knowledge, inaccurate measurements, and the inability to measure parameters at each spatial or temporal location.Although the deterministic Henry problem is well-known, the difficulties here are that it is not clear how input uncertainties propagate through the non-linear, time-dependent problem.The solution to be found is then the mean value and the variance of the salt mass fraction, both are functions of space and time.An accurate estimate of the output uncertainties can facilitate a better understanding of the problem, better decisions, and improved control and design of the experiment.From an implementation point of view, we test how easy it is to couple multigrid and multilevel Monte Carlo solvers.
Salinisation of coastal aquifers.The intrusion of salt water occurs when the sea level rises and the salt water moves onto the land [71].This usually happens during storms, floods, droughts or when saltwater intrudes into freshwater aquifers and raises the water table.As groundwater is an essential resource for food and irrigation, its salinisation can be disastrous.Many hectares of farmland could be lost because it becomes too wet or too salty to grow crops.Therefore, accurate modelling of different saline flow scenarios is essential [1,71] to help farmers and researchers develop strategies to improve soil quality and reduce the impact of saltwater intrusion.
Modeling.The saltwater flow is density driven and described by a system of time-dependent nonlinear partial differential equations (PDEs).It is characterised by convection dominance and can exhibit very complicated behaviour [81].Uncertain parameters can have a strong influence on the flow and transport of salt.Random fields are used to model these parameters.
There are a number of studies where authors have modelled uncertainties in reservoirs (see [12,79]).The link between stochastic methods and hydrogeological applications was made in [8], where the authors analysed a collaboration between academics and water suppliers in Germany and made recommendations for optimisation and risk assessment.The basics of stochastic hydrogeology and an overview of stochastic tools and uncertainty management are described in [68].
The review [77] deals with hydrogeological applications of recent advances in uncertainty quantification, probabilistic risk assessment and decision making under uncertainty.The author reviews probabilistic risk assessment methods in hydrogeology under parametric, geological and model uncertainties.Density-driven vertical transport of saltwater through the freshwater lens on the island of Baltrum (Germany) is modelled in [62].
In [41], the authors examined the implications of transgression for a range of seawater intrusion scenarios based on simplified coastal freshwater aquifer settings.They stated that vertical intrusion during transgressions could involve density-driven convective processes, causing substantially greater amounts of seawater to enter the aquifer and create more extensive intrusion than horizontal seawater intrusion in the absence of transgression.
History and preliminary results.This work is an extension and continuation of our preliminary works [44,43].The original Henry saltwater intrusion problem was introduced by H.R. Henry in the 1960s (see [36]) and became a benchmark for numerical solvers for groundwater flow (see [81,73,72,20].In [66], the authors use the generalised polynomial chaos expansion approximation to investigate how incomplete knowledge of system properties affects the assessment of global quantities.In particular, they estimated the propagation of input uncertainties into a few dimensionless scalar parameters. Hydrogeological formations typically have complex and heterogeneous structures.These formations may consist of several layers of porous media with different porosity and permeability coefficients (cf.[64,70]).Measurements of the layer positions and their thicknesses are only possible with a certain degree of error, and average parameters are typically assumed for the materials within the layers.These layers are therefore excellent candidates for random field modelling.Furthermore, due to the non-linearities in the problem, computing with the averaged parameters does not necessarily lead to the mathematical expectation of the solution that can be obtained by averaging over all scenarios.
Methods.Many techniques can be used to quantify uncertainties.One classical method is Monte Carlo (MC) sampling.Although it is dimension independent, it converges very slowly, with a convergence rate of O( 1 √ N ).This method may not be affordable for time-consuming simulations.However, even modern techniques such as surrogate models and stochastic collocation require a few hundred to a few thousand time-consuming simulations and assume a certain smoothness of the quantity of interest (QoI).
Another class of methods is the class of perturbation methods [19].The idea is to decompose the QoI with respect to (w.r.t.) random parameters in a Taylor series.The higher order terms can be neglected for small perturbations, simplifying the analysis and numerical tests.These methods assume that the random perturbations are small.For larger perturbations, these methods usually do not work.
Other methods to compute the desired statistics of the QoI are direct integration methods, such as the quasi-MC (QMC) [11], collocation methods and surrogate-based (generalized polynomial chaos approximation [83,46,45,76,51] and stochastic Galerkin [3,30,23]) methods.Direct methods compute statistics directly by sampling uncertain input coefficients and solving the corresponding PDEs, whereas the surrogate-based method computes a cheap functional (polynomial, exponential, or trigonometrical) approximation of the QoI.Examples of the surrogate-based methods are radial basis functions [49,9,50,31], sparse polynomials [14,7,21], and polynomial chaos expansion [54,17,82].Surrogate methods are using some well-known functions such as the multivariate Legendre, Hermite, Chebyshev or Laguerre functions [58,83].Advantages of surrogate methods are: Once the model is constructed, it is easy to sample, and all the samples are almost free (much cheaper than sampling the original stochastic PDE).The non-trivial part of surrogate models is to define which polynomial order is needed and how accurately all coefficients should be computed.Another difficulty is that not every function can be well approximated by a polynomial.
The quantification of uncertainties in stochastic PDEs could be a significant challenge due to a) the large possible number of involved random variables and b) the high cost of each deterministic solution of the governed PDE.The MC quadrature and its variance-reduced variants have a dimension-independent error convergence rate O(N − 1 2 ), and the QMC has the worst-case rate O(log M (N )N −1 ), where N is the number of samples, and M indicates the dimension of the stochastic space [52].The MC method is not affected by the dimension of the integration domain, such as collocations on sparse or full grid methods [2,55].A numerical comparison of other QMC sequences is presented in [63].
This work has the following structure.Section 2 describes the Henry problem and numerical methods to solve it.The well-known MLMC method is reviewed in Section 3. Next, Section 4 details the numerical results, which include the numerical analysis of the Henry problem, the computation of different statistics, the performance of the MLMC method, and the practical performance of the parallel ug4 solver for the Henry problem [36,73] with uncertain coefficients.Finally, we conclude this work with a discussion in section 5.
Our contribution and main results: We coupled the MLMC method with the multigrid method to estimate the propagation of uncertainties in the Henry problem.We used random fields to model unknown porosity, permeability and recharge.In our model, porosity and permeability are spatially dependent, random, multi-scale and have two layers.The recharge is time dependent and uncertain.We have not simplified our equations with the Boussinesq approximation, i.e. we have considered a more general case.We used MLMC to compute the mean and variance of the mass fraction.In addition, several other quantities of interest were computed, such as the total freshwater integral, the salt integral, the mass fraction computed at a point and integrated over a small subdomain.We investigated which of these QoIs were more appropriate to be computed by the MLMC method.To further speed up the computational process, we run all the MLMC random simulations and each individual simulation in parallel.This work is not only of theoretical interest, but also of practical interest; our modelling and solution allow us to answer the following questions: 1. How long can a particular drinking water well be used (i.e. when will the mass fraction of salt exceed a critical threshold)?
2. Which part of the aquifer has particularly high uncertainty?
3. What is the probability that the salt concentration at a given location and time will exceed a threshold value?
4. What is the mean scenario (and its variations)?
5. What are the extreme scenarios for the concentration of the salt?

How do the uncertainties change with time?
To the best of our knowledge, we are not aware of any other studies where Henry's problem [36,73] has been solved in parallel using MLMC methods with uncertain porosity, permeability and recharge parameters.

Problem setting
In coastal aquifers, saline seawater intruding into the formation from one side (the sea side) displaces pure water from the other side due to water recharge from land sources and precipitation.Because of its higher density, seawater tends to infiltrate along the bottom of the aquifer.This process can reach a steady state, but may be time-dependent due to the periodicity of the recharge or the control of the pumping rate from the wells.Accurate simulation of salinisation is essential for predicting the availability of water resources.However, the accuracy of such predictions strongly depends on the hydrogeological parameters of the formation and the geometry of the computational domain, denoted by D.
The aquifer D ⊂ R d , d ∈ {2, 3}, can be modelled as an immobile porous matrix filled with liquid phasea solution of salt in water.Due to the inhomogeneous density distribution, gravity induces the motion of the liquid phase.This motion transports the salt, which would otherwise be subject to molecular diffusion.
A simple but very illustrative model of coastal aquifers is the so-called Henry problem, first considered in [36].In this two-dimensional setting, the aquifer is represented by a rectangular domain completely saturated with the liquid phase (Fig. 1).The salty seawater intrudes from the right and the pure water recharges from the left.The top and bottom are assumed to be impermeable.Analogous settings with partially saturated domains are considered in [75].
The mass conservation laws for the whole liquid phase and the salt give the following equations where p = p(t, x) : [0, +∞) × D → R is the hydrostatic pressure, µ = µ(c) denotes the viscosity of the liquid phase, and g = (0, . . ., 0, −9.8) T ∈ R d represents the gravity vector.Inserting (3) into (1-2) results in a system of two time-dependent PDEs in the unknowns c and p.This system should be closed with boundary conditions for c and p and an initial condition for c.Following the classical setting in [81], for this variant of the Henry problem (see also [72,73]), we set and D = ϕDI, with a scalar D ∈ R, and I ∈ R d×d the identity matrix. (5) Furthermore, we assume the isotropic permeability This setting is consistent with the problem setting in [81].However, we do not assume the Boussinesq approximation and keep the density variable for all terms.
For the initial conditions, we set c| t=0 = 0.
The boundary conditions are presented in Fig. 1 (left).On the right side of the domain, we impose Dirichlet conditions for c and p variables that represent the adjacent seawater aquifer: Parameter Values and Units Description scaling in the Kozeny-Carman Table 1: Parameters of the considered density-driven flow problem On the left side, we prescribe the inflow of fresh water: where e x = (1, 0) ⊤ , and qin is a constant.For the classical formulation of the Henry problem, this value was set to qin = 6.6 • 10 −2 kg/s in [81] or qin = 3.3 • 10 −2 kg/s in [73,72].The Neumann zero boundary conditions are imposed on the upper and lower sides of D.
An example of c(t, x) and the flow for the parameters from Table 1 (i.e.ϕ = φ and K = K) is shown in Fig. 1(right).The dark red colour corresponds to c = 1 (salt water) and the dark blue colour corresponds to c = 0 (pure water).Due to its higher density, the saltwater intrudes into the aquifer at the bottom right.It is pushed back by the lighter pure water coming from the left.This process creates a vortex in the flow at the bottom right of the domain.The salt water enters at the lower part of the right boundary and diverts upwards and to the right, back to the sea, forming a salt triangle.This flow does not transport the salt to the left part of the domain.The salt moves further to the left by diffusion and dispersion and is washed out by the recharge.In the classical formulation, this salt triangle initially grows with time but reaches a steady state (see [81,73,72]).However, the initial non-stationary phase can take considerable time.Investigation of this phase is particularly important to understand the behaviour of the system when the intensity of recharge varies, as may occur, for example, due to climate change.We consider integrals over the whole domain D, which describes the total amount of pure water and the total amount of salt (as in equations ( 35)-( 36)).We also compute local integrals over 15 small rectangular subdomains ∆ 1 , . . ., ∆ 15 , around preselected points defined below.The list of selected points follows: (see Fig. 2.) The reason we consider 15 points (and 15 small subdomains) is that not all points are "interesting", i.e., not all points have significant variation in c.MLMC reduces the variance, but if the initial variance is smaller according to the prescribed tolerance, there is no point in using MLMC.
The knowledge of c at these spatial points may help track salinity changes over time in groundwater wells and understand which areas in the aquifer are most vulnerable.

Modeling porosity, permeability, and recharge
The primary sources of uncertainty are the hydrogeological properties of the porous medium -porosity (ϕ) and permeability (K) fields of the solid phase -and the freshwater recharge flux qx through the left boundary.The QoIs are related to the mass fraction c, a function of ϕ, K, and the recharge.We model the uncertain ϕ using a random field and assume K to be isotropic and dependent on ϕ: The distribution of ϕ(x, ξ), x ∈ D, is determined by a set of stochastic parameters ξ = (ξ 1 , . . ., ξ M ).Each component ξ i is a random variable depending on a random event ω.For concision, we skip ω and write ξ := ξ(ω).
The dependence in (10) is specific for every material.We refer to [59,60,18] for a detailed discussion.In the proposed model, we use a Kozeny-Carman-like dependence where κ KC is a constant scaling factor.The inflow flux is kept constant across the left boundary but depends on the stochastic variable q in .We also assume that the inflow flux is independent of ϕ and K.
Remark 1 Equation ( 5) is a simplification typically adopted when the Henry problem is used as a benchmark, see [81,72,73].It is then assumed that D is not only the molecular diffusion coefficient but also includes some dispersivity.In a more accurate model, the dispersion tensor should depend on the Darcy velocity q, for example using the Bear-Scheidegger law (see [69,5]).However, we note that the parameters of mechanical dispersion depend on the configuration of the porous matrix (e.g. grain size) and therefore cannot be considered independent of porosity and permeability in stochastic modelling.The influence of dispersion on the propagation of uncertainties is therefore a non-trivial issue even from a modelling point of view.However, it should be addressed in the future.

Numerical methods for the deterministic problem
The system (1-2) is numerically solved in the domain D × [0, T ], where × denotes the Cartesian product.D is covered with a grid D h of quadrilaterals with mesh size h.Equations (1-2) are discretized on D h using a vertex-centered finite-volume scheme as presented in [24,25,26].The number of degrees of freedom associated with D h is denoted by n.There are two degrees of freedom per grid vertex in D h : one for the mass fraction c and another for the pressure p.Note that h = O(n −1/d ), d = 2.We use the implicit Euler method with a time step τ for time discretization.The number of the computed time steps is r = T /τ .Solution of the discretized system yields the approximation c h,τ for c on this grid.We use partial upwind for the convective terms (cf.[24]).Therefore, the discretization error is of the second order w.r.t. the spatial mesh size h.However, Equation ( 2) is essentially convection dominated.For the grids in the numerical experiments, the observed reduction of the discretization error after grid refinement corresponds to the first order.Furthermore, the Euler method provides the first-order discretization error w.r.t.τ .Thus, we assume the first-order dependence of the discretization error w.r.t.h and τ , i.e., as d = 2, which is consistent with our numerical experiments.The implicit time-stepping scheme provides unconditional stability but requires the solution of an extensive nonlinear algebraic system of the discretized equations with n unknowns in every time step.The Newton method is used to solve this system.Linear systems inside the Newton iteration are solved using the BiCGStab method (cf.[4]) preconditioned with the geometric multigrid method (V-cycle, cf.[33]).In the multigrid cycle, the ILU β -smoothers [34] and Gaussian elimination are used as the coarse grid solver.
To construct the spatial grid hierarchy D 0 , D 1 , . . ., D L , we start with a coarse grid D 0 consisting of 512 grid elements (quadrilaterals) and n 0 = 1122 degrees of freedom.(n ℓ denotes the number of degrees of freedom in D ℓ .)Grid D ℓ of freedom is regularly refined twice to obtain grid D ℓ+1 .We treat this as one spatial refinement step.After every spatial grid refinement step, the number of grid elements is multiplied by 16, i.e., n ℓ ≈ n 0 • 16 ℓ , see Table 2).We also construct the temporal grid hierarchy T 0 , T 1 , . . ., T L .The time step on each temporal grid is denoted by τ ℓ with τ ℓ+1 = 1 4 τ ℓ .The number of time steps on the ℓth grid (level) is r ℓ+1 = 4r ℓ , so that r ℓ = r 0 4 ℓ , where r 0 is the number of time steps in T 0 .On the ℓth level, the MLMC uses the grid D ℓ × T ℓ .According to (12), for the numerical solution c ℓ = c h,τ on D ℓ × T ℓ , we have: Up to four grid levels ℓ were used in the numerical experiments.
In the context of this work, it is critical to estimate the numerical complexity of the deterministic solver with respect to the grid level ℓ.The most time consuming part of the simulation is the solution of the discretised nonlinear system.It is usually difficult to predict the number of Newton iterations in each time step.But in our numerical experiments, two Newton iterations were sufficient to achieve the required accuracy independently on the grid level.Accordingly, the linear solver was called at most twice per time step.Furthermore, the convergence rate of the geometric multigrid method does not depend on the mesh size (cf.[34]).Thus the computational complexity of a time step is O(n ℓ ), where n ℓ is the number of degrees of freedom at grid level ℓ.Therefore, the total numerical cost of computing a scenario at grid level ℓ for r ℓ time steps is

Multilevel Monte Carlo
To reduce the total computing cost, we apply the MLMC method, which is a natural idea because the deterministic solver uses a multigrid method (see Section 2.3).The MLMC method efficiently combines samples from various levels.A more in-depth description of these techniques is found in [15,16,28,29,35,78,48].We let ξ(ω) and g(ξ) = g(ξ(ω)) represent a vector of random variables and the QoI, respectively, where ω is a random event.In this work, g is an integral over some of ∆ i 's or D with the integrand depending on c, s.Section 4. The MLMC method aims to approximate the expected value E [g] with an optimal computational cost.It constructs a telescoping sum, defined over a sequence of spatial and temporal meshes, ℓ = 0, . . ., L, as described next, to achieve this goal.The QoI g, numerically evaluated on level ℓ, is denoted by g h ℓ ,τ ℓ ,ℓ or, for simplicity, by just g ℓ , where h ℓ and τ ℓ are the discretization steps in space and time on level ℓ.Furthermore, we assume that E [g h,τ ] → E [g] as h → 0 and τ → 0.
Remark 2 Since the problem described in Subsection 2.1 is non-stationary, the quantities related to its solution c below depend on time, too.Examples of time-dependent quantities are g, g ℓ , To keep our notation simple, we do not explicitly specify the time argument for them.We shall only indicate it where it is essential.
Let s 0 be the computing cost to evaluate one realization of g 0 (the most expensive one from all computed realizations on 0th mesh).Similarly, s ℓ denotes the computing cost of evaluating g ℓ − g ℓ−1 .For simplicity, we assume that s ℓ for g ℓ − g ℓ−1 is almost the same as s ℓ for g ℓ .The number of iterations is variable; thus, the cost of computing a sample of g ℓ − g ℓ−1 may fluctuate for various realizations.
The MLMC method calculates E [g L ] ≈ E [g] using the following telescopic sum: In the above equation, level ℓ in the superscript (ℓ, i) indicates that independent samples are used at each correction level.As ℓ increases, the variance of g ℓ − g ℓ−1 decreases.Thus, the total computational cost can be reduced by taking fewer samples on finer meshes.
Furthermore, we denote V 0 = V [g 0 ] and for ℓ ≥ 1, let V ℓ be the variance of g ℓ − g ℓ−1 : The standard theory states the following facts for the mean and for the variance: The cost of the multilevel estimator Y is In the following, we repeat the well-known [29] results on the computation of the sequence m 0 , . . ., m L .For a fixed variance V [Y ] =: ε 2 /2, the cost S is minimized by choosing as m ℓ the solution of the optimization problem min m0,...,m L F (m 0 , . . ., m L ), where F (m 0 , . . ., m L ) := L ℓ=0 m ℓ s ℓ + µ 2 V ℓ m ℓ , µ 2 is a Lagrange multiplier.Thus, at the desired set of m ℓ , the derivatives of F w.r.t.m ℓ are equal to zero: Solving the system (22), we obtain Taking into account that the variation V [Y ] is fixed and substituting ( 23) into (19), i.e.
L ℓ=0 V ℓ /m ℓ = ε 2 /2, we obtain an equation for µ: From this equation, we get µ = 2ε −2 L ℓ=0 √ V ℓ s ℓ , and therefore For this set of m ℓ , the total computational cost of Y is For further analysis of this sum, see [29], p.4.The mean squared error (MSE) is used to measure the quality of the multilevel estimator: where Y is what we computed via MLMC, and E [g] what actually should be computed.To achieve for some prescribed tolerance ε, we ensure that and The bias error (27) corresponds to the discretization error (13) discussed in Subsection 2.3.Later, in the numerical Section we will see that The bias error can be made smaller than ε 2 /2 by choosing a sufficiently large L.Then, for this L, we can compute optimal m 0 , . . ., m L by formula in (24) to provide (28).
Combining this idea with a sequence of levels of the simulation grids in which the cost increases exponentially with ℓ while the weak error E [g L − g] and multilevel correction variance V ℓ decrease exponentially leads to the following theorem (cf.Theorem 1, p. 6 in [29]): Theorem 4 Consider a fixed t = t * .Suppose positive constants α, β, γ > 0 exist such that α ≥ 1 2 min(β, γ d), and Then, for any accuracy ε < e −1 , a constant c 4 > 0 and a sequence of realizations {m ℓ } L ℓ=0 exist, such that MSE < ε 2 , where MSE is defined in (26), and the computational cost is This theorem (see also [39,38,13,15,28]) indicates that, even in the worst-case scenario, the MLMC algorithm has a lower computational cost than that of the traditional (single-level) MC method, which scales as O(ε −2− dγ/α ).
Remark 5 In Theorem 4, the factors c 1 , c 2 , c 3 , c 4 as well as the exponents α, β and γ may generally depend on the time point t * .This makes L and m ℓ time-dependent, too.In this work, we assume that there exist upper bounds ĉi for the factors: c i ≤ ĉi , i = 1, . . ., 4, as well as the bounds for the exponents: α ≥ α, β ≥ β, γ ≤ γ over the whole time interval [0, T ] (later, we will see that this assumption is confirmed by our numerical tests).Then the corresponding L and m ℓ computed from (27) and (24) are constant in time, too.For simplicity, in what follows, we omit the hats and refer to these bounds.
Remark 6 In ( 27) and (28), ε is supposed to have the units of QoI g and to have the corresponding meaning.For example, if g is the mass of the salt in a subdomain, measured in kg, then ε bounds the error of this mass.This raises the question of appropriate scaling of ε. Let In this work, it makes sense to consider the error relatively to g.For this, in ( 27) and (28), we replace ε by ε • E 0 .Equivalently, we can divide (rescale) E [g L − g] by E 0 and V i by E 2 0 .Therefore, we get: Now, to satisfy (27), we want From this inequality, we can estimate L: Equations ( 24) and ( 25) attain the form Note that although E 0 is an inaccurate approximation of E [g], it is sufficient to be used in (32)-( 33) for scaling purposes.In contrast to E [g ℓ ] with ℓ > 0, E 0 can be better estimated by the MLMC method as a large number of samples are computed on the grid level 0.
Using preliminary numerical tests (see Fig. 10), we can estimate the convergence rates α for the mean (the so-called weak convergence) and β for the variance (the so-called strong convergence), as well as the constants c 1 and c 2 .In addition, α is strongly connected to the order of the discretization error (see Section 2.3), which equals 1.Note that precise estimates of parameters α and β are crucial to distribute the computational effort optimally.

Numerical Experiments
In this section, we perform numerical tests with the MLMC method described above.One of our aims is to compare the theoretical predictions with the obtained numerical results.We present two types of numerical experiments.The first type (denoted by A1 -A3) demonstrate the solution in the whole domain, in a point, and salt and fresh water integral values: where I(•) is the indicator function identifying a subdomain {x : c(t, x, ω) ≤ 0.012178}.For these integrals, we observe numerically α ≈ 2 (see (29)) that agrees with the theoretical value.The idea of these tests is to collect more information about the problem.
Remark 7 Our motivation for the threshold 0.012178 (in (35)) is technical: it was used in our simulations of coastal aquifers at the North Sea [71], where 1m 3 of saltwater contains 35kg of salt (corresponding to c = 1 in the model) and has a density of 1035 kg m 3 .Then the maximum recommended salinity of drinking water 412 mg l corresponds to the scaled mass fraction c = 0.012178.The choice of this value is not essential.
In the second type (denoted by B1-B4), as a QoI, we consider an integral over a small sub-regions around a point (see the list of points in ( 9)): The value of Q i is the mass of the salt in a subdomain ∆ i .The size of each ∆ i is small (0.2 2 = 0.04), compared to D. We choose the point with index 9 (see Fig. 2) and the corresponding domain ∆ 9 as an example.The numerical scheme provides only the first order of the accuracy to compute Q i (t, ω), i.e., the convergence rate α (weak convergence) should be ≈ 1 in (29a).
Uncertain porosity and recharge: Often unknown porosity is modelled by a random field.This random field is approximated by a truncated Karhunen-Loéve expansion (KLE).KLE requires knowledge of the covariance matrix, which is typically assumed to be from the very large Matérn class of functions.The KLE components are computed by solving an auxiliary eigenvalue problem and assuming a random distribution for the random variables.This approach is difficult and requires a separate article to describe it.We have implemented a simple alternative.We used L 2 orthogonal functions (sin() and cos()), which mimic the orthogonal functions used in KLE.The uniform random variables used are also very common, usually one of three options is chosen: uniform, Gaussian or log-normal.
Additionally, we assume the presence of two horizontal layers: y ∈ (−0.8, 0] (the upper layer) and y ∈ [−1, −0.8] (the lower layer).The porosity inside each layer is uncertain and is modeled as in (37): where Random fields generated in this way provide: a) periodicity of high and low porosity areas in the aquifer; b) multi-scale behaviour; c) spatial dependence.This approach has been agreed with the developers of the UG4 software, who have many years of experience in solving very realistic problems at an industrial level.The recharge flux is also uncertain and is equal to Here ξ 1 , ξ 2 , and ξ 3 are sampled independently and uniformly in [−1, 1].
Note that there are many ways to model recharge, and we chose to model it by a function of t with one random variable ξ 3 .Of course, one can use more random variables to model the recharge, but we are limited by the available computational resources (i.e., by the number of independent random variables which we can use).Therefore, we decided to use only one RV.Taking into account a periodic nature of precipitations, we decided to model it by a periodic function, for example, by sin(•).
In the settings of the classical Henry problem, the flux qin is simply prescribed.It models the flux of the water coming from the "land" (due to rivers, precipitation etc.) to the coast.This means in the real situation, on the one hand, qin should depend on the hydrogeological parameters of the medium further on the left, behind the left boundary.And for these parameters, the same stochastic modeling should be applied.Furthermore, this flux should depend on the flow inside the domain, i.e. on the porosity and the permeability in the problem setting, and not only near the left boundary.These arguments reflect the general difficulties with imposing realistic in-and outflow boundary conditions but their discussion is beyond the scope of our paper.This is a modeling issue that must be considered aside of the numerical methods.
Note that the porosity (37) and therefore the permeability (11) are highly oscillating functions of the geometric coordinates (x, y).Input interfaces of some codes assume precalculation of these values and storing them in e.g.vertex-or cell-centered rasters.During assembling of the discretized algebraic systems, these quantities are interpolated at the integration points of the finite-volume or finite-element schemes.For multilevel approaches, these points depend on the grids in the hierarchy, and the resolution of the rasters as well as the type of the interpolation can influence the properties of the methods.In constrast to that, our implementation does not involve precalculation of the porosity and the permeability fields, so that the values (37) and (11) are directly evaluated at the integration points during the assembling phase.The same holds for the recharge (41) as a function of time.

Software and parallelization:
The computations presented in this work were performed using the ug4 simulation software toolbox [37,65,80].This software has been applied for subsurface flow simulations of real-world aquifers (cf.[71]).The toolbox was parallelized using MPI, and the presented results were obtained on the Shaheen II cluster provided by the King Abdullah University of Science and Technology.Every sample was computed on 32 cores of a separate cluster node.Each simulation (scenario) was localized to one node to reduce the communication time between nodes.All scenarios were concurrently computed on different nodes.A similar approach was used in [45,46].Simulations times for one realization on different grid levels ℓ are presented in Table 2.

Computation via QMC
The following calculations are performed using the QMC method with 600 samples (Halton sequence).These numerical tests are not directly relevant to the MLMC method, but are used to better understand the uncertainties in different QoIs.
Test A1.The mean and variance of the mass fraction are provided in Fig. 4   We observed that the variability (uncertainty) of the mass fraction can vary from one grid point to another.At some points (dark blue regions) the solution does not change.At other points (white-yellow regions) the variability is low (light red regions) or high (dark red regions).In regions of high uncertainty, it is useful to refine the mesh and apply the MLMC method.
Test A2.As a further numerical test, we want to know the range in which the mean and variance of the QoI g ℓ change.In Fig. 5

Numerical tests for the MLMC method
The following numerical experiments are required to demonstrate the work and efficiency of the MLMC method.Test B1.The next QoI is the integral value as in (36).This integral can be estimates from above  In both figures, as we move to finer and finer meshes, we can see that the amplitude decreases, i.e. the accuracy increases.Figure 8 shows 100 realisations of g 1 − g 0 (left), g 2 − g 1 (center), g 3 − g 2 (right).Here QoI g ℓ is the integral value as in (42) computed over a subdomain ∆ 9 for t ∈ [τ, 48τ ].The amplitude of g 1 − g 0 (on the left) achieves ≈ 19, of g 2 − g 1 (in the middle) achieves ≈ 5.5, and the amplitude of g 3 − g 2 (on the right) achieves 1.2 (i.e.we observe decay by factor ≈ 4).Table 2 contains average computing times, which are necessary to estimate the number of samples m ℓ at each level ℓ.The fifth column contains the average computing time, and the sixth and seventh columns contain the shortest and longest computing times.The computing time for each simulation varies depending on the number of iterations, which depends on the porosity and permeability.We observed that, after ≈ 6016 s, the solution is almost unchanging; thus, we perform the experiment only for t ∈ [0, T ], where T = 6016.For example, if the number of time steps is r ℓ = 94 (Level 0 in Test B2.The QoI g ℓ is the integral Q 9 at t = 10τ (see (36)).In Figure 9 we compare the decay of E [g ℓ ] and E [g ℓ − g ℓ−1 ] (left) and of V [g ℓ ] and V [g ℓ − g ℓ−1 ] (right).We can see that the decay of g ℓ − g ℓ−1 is much faster.3 and in Fig. 11.We can see that MLMC outperforms MC by a factor of ≈ 30 − 3000 depending on the MSE accuracy ε 2 .We can also see the estimated number of required levels L and the numbers of samples {m 0 , m 1 , m 2 , m 3 } for each accuracy ε 2 .These results are also visualized in Fig. 11.Additionally, to the MLMC and MC curves, we added two ) (see (30) from Theorem 4 for β < dγ, with β ≈ 1.7, d = 3, γ = 1, α = 0.9).The second (yellow line) corresponds to ε −2 , and we labelled it "optimal".Thus, the MLMC slope is 2 + d•γ−β α = 2 + 3•1−1.7 0.9 = 3.44, and the MC slope ε −2−3•1/0.9 ) = O(ε −5.3 ).Here γ = 1 and d = 3 (spatial dimension 2 plus temporal dimension 1), and the required condition 0.9 ≥

Conclusion
We investigated the applicability and efficiency of the MLMC approach to the Henry-like problem with uncertain porosity, permeability and recharge.These uncertain parameters were modelled by random fields with three independent random variables.Permeability is a function of porosity.Both functions are time-dependent, have multi-scale behaviour and are defined for two layers.The numerical solution for each random realisation was obtained using the well-known ug4 parallel multigrid solver.The number of random samples required at each level was estimated by calculating the decay of the variances and the computational cost for each level.These estimates depend on the minimisation function defined in (21).
The MLMC method was used to compute the expected value and variance of several QoIs, such as the solution at a few preselected points (t, x), the solution integrated over a small subdomain, and the time evolution of the freshwater integral.We have found that some QoIs require only 2-3 mesh levels and samples from finer meshes would not significantly improve the result.Other QoIs require more grid levels.Note that a different type of porosity in (37) may lead to a different conclusion.
The numerical results confirm that the computational cost of the MLMC method is lower than that of the MC method.Therefore, sampling at different mesh levels is useful and helps to reduce the overall computational cost.Constraints.1.It may happen that the QoIs calculated at different mesh levels are very similar.In this case the default (Q)MC is sufficient.2. The time dependence is challenging: the optimal number of MLMC samples depends on the QoI.For example, this number may be large for the solution computed at one point (t, x) and small for the solution computed at another time.3.Each new QoI requires new estimates of all parameters α, β, γ, and the resulting MLMC graphs may be very different.Future work. 1.It would be beneficial to consider a more complicated/multiscale/realistic porosity and permeability with more random variables.2. A more advanced version of MLMC may give better estimates of the number of levels L and the number of samples on each level m ℓ .3. Known experimental and measured data of porosity, permeability, velocity or mass fraction could be used to identify unknown parameters and minimise uncertainties [42,47,53,67].

Figure 2 :
Figure 2: Positions of the 15 pre-selected points with small subdomains around them.The size of the subdomain around point (x i , y i ) is [x i − 0.1, x i + 0.1] × [y i − 0.1, y i + 0.1].
on the left and right, respectively.The expectation takes values from [0, 1].The value 0 corresponds to the blue colour and the value 1 to the dark red.The variance range is [0, 0.04].The areas with high variance (dark red) indicate regions with high variability/uncertainty.Such regions may need additional attention from specialists (e.g., placement of sensors).The right image displays four contour lines {x : Var[c](t, x) = 0.01 • i}, i = 1..4, t = T = 6016.

Figure 11 :
Figure 11: Comparison of MC and MLMC for different values of ε (x-axes).

Table 2 :
Number of degrees of freedom n ℓ , number of time steps r ℓ , step size in time τ ℓ , average, minimal, and maximal computing times on each level ℓ.

Table 3 :
Comparison of MC and MLMC and the number of samples on each level vs ε.
, D 1 , . . ., D L hierarchy of spatial meshes T 0 , T 1 , . . ., T L hierarchy of temporal meshes s ℓ , S complexity on level ℓ, total complexity h ℓ (or h), n ℓ spatial step size and number of spatial degrees of freedom on level ℓ τ ℓ (or τ ), r ℓ time step size and number of time steps on level ℓ m