Uncertainty Quantification in Coastal Aquifers Using the Multilevel Monte Carlo Method

1. Introduction

Notation
QoI g	quantity of interest g
$\mathcal{D}$	computational spatial domain
${\mathcal{D}}_0,{\mathcal{D}}_1,\dots ,{\mathcal{D}}_L$	hierarchy of spatial meshes
${\mathcal{T}}_0,{\mathcal{T}}_1,\dots ,{\mathcal{T}}_L$	hierarchy of temporal meshes
L	number of levels
s	complexity
$h_{\ell }$ (or h), $n_{\ell }$	spatial step size and number of spatial degrees of freedom on level ℓ
${\tau }_{\ell }$ (or $\tau$), $r_{\ell }$	time step size and number of time steps on level ℓ
$m_{\ell }$	number of samples (scenarios) on level ℓ
$\mathbb{E}\left[·\right]$, $\mathbb{V}\left[·\right]$	expectation and variance
$\Theta$	multidimensional domain of integration in parametric space
$\omega$, $\mathbf{\xi }\left(\omega \right)$	random event and random vector
$\varphi (\mathbf{x},\omega )$	porosity random field
$\mathbf{K}(\mathbf{x},\omega )$	permeability random field
$\rho (\mathbf{x},\omega )$	density random field
$\mathbf{q}(t,\mathbf{x},\omega )$	volumetric velocity
$\mathbf{D}$	tensor field $\mathbf{D}=\mathbf{D}\left(\mathbf{q}\right)$: molecular diffusion and dispersion of salt
$\overline{\kappa }\left(\mathbf{x}\right)$	expectation of $\kappa (\mathbf{x},\omega )$
d	physical (spatial) dimension
$c=c(t,\mathbf{x},\omega )$	mass fraction of salt (solution of the problem)

Saltwater intrusion occurs when sea levels rise and saltwater moves onto the land. Usually, this occurs during storms, high tides, droughts, or when saltwater penetrates freshwater aquifers and raises the groundwater table. Since groundwater is an essential nutrition and irrigation resource, its salinization may lead to catastrophic consequences. Many acres of farmland may be lost because they can become too wet or salty to grow crops. Therefore, accurate modeling of different scenarios of saline flow is essential [1,2] to help farmers and researchers develop strategies to improve the soil quality and decrease saltwater intrusion effects.

Saline flow is density-driven and described by a system of time-dependent nonlinear partial differential equations (PDEs). It features convection dominance and can demonstrate very complicated behavior [3].

As a specific model, we consider a Henry-like problem with uncertain permeability and porosity. These parameters may strongly affect the flow and transport of salt. The original Henry saltwater intrusion problem was introduced by H.R. Henry in the 1960s (cf. [4]). The Henry problem became a benchmark for numerical solvers for groundwater flow (cf. [3,5,6,7]. In [8], the authors use the generalized polynomial chaos expansion approximation to investigate how incomplete knowledge of the system properties influences the assessment of global quantities. Particularly, they estimated the propagation of input uncertainties into a few dimensionless scalar parameters.

The hydrogeological formations typically have complicated and heterogeneous structures. These formations may consist of a few layers of porous media with various porosity and permeability coefficients (cf. [9,10]). Measurements of the layer positions and their thicknesses are only possible up to some error, and for the materials inside the layers, the average parameters are typically assumed. Thus, these layers are excellent candidates to be modeled by random fields. Further, due to the nonlinearities in the problem, averaging the parameters does not necessarily lead to the correct mathematical expectation of the solution.

To model uncertainties, we use random fields. Uncertainties in the input data propagate through the model and make the solution (e.g., the mass fraction) uncertain. An accurate estimation of the output uncertainties can facilitate a better understanding of the problem, better decisions, and improved control and design of the experiment.

The following questions can be answered:

• How long can a particular drinking water well be used (i.e., when will the mass fraction of the salt exceed a critical threshold)?
• What regions have especially high uncertainty?
• What is the probability that the salt concentration is higher than a threshold at a certain spatial location and time point?
• What is the average scenario (and its variations)?
• What are the extreme scenarios?
• How do the uncertainties change over time?

Many techniques can quantify uncertainties. A classical method is Monte Carlo (MC) sampling. Although it is dimension-independent, it converges very slowly and requires many samples. This method may not be affordable for time-consuming simulations. Nevertheless, even up-to-date techniques, such as surrogate models and stochastic collocation, require a few hundred to a few thousand time-consuming simulations and assume a certain smoothness by the quantity of interest (QoI).

Another class of methods is the class of perturbation methods [11]. 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 numerics. These methods assume that random perturbations are small (e.g., up to 5% of the mean, depending on the problem). For larger perturbations, these methods usually do not work.

There are quite a few studies where authors model uncertainties in reservoirs (cf. [12,13]). Reconnecting stochastic methods with hydrogeological applications was accomplished in [14], where the authors analyzed a collaboration between academics and water suppliers in Germany and made recommendations regarding optimization and risk assessment. The fundamentals of stochastic hydrogeology and an overview of stochastic tools and accounting for uncertainty are described in [15].

The review [16] deals with hydrogeologic applications of recent advances in uncertainty quantification, probabilistic risk assessment, and decision-making under uncertainty. The author reviewed probabilistic risk assessment methods in hydrogeology under parametric, geologic, and model uncertainties. Density-driven vertical transport of saltwater through the freshwater lens on the island of Baltrum (Germany) is modeled in [17].

In [18], 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.

The methods to compute the desired statistics of the QoI are direct integration methods, such as the MC, quasi-MC (QMC) and collocation methods and surrogate-based (generalized polynomial chaos approximation and stochastic Galerkin [19,20,21,22]) 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 [23,24,25,26], sparse polynomials [27,28,29], and polynomial chaos expansion [30,31,32]. Sparse grid methods to integrate high-dimensional integrals are considered in [31,33,34,35,36,37,38,39,40]. An idea to generate goal-oriented adaptive spatial grids and use them in the multilevel MC (MLMC) framework was presented in [41,42].

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 $\mathcal{O}\left(N^{-\frac{1}{2}}\right)$, and the QMC has the worst-case rate $\mathcal{O}\left({\log }^M\left(N\right)N^{-1}\right)$, where N is the number of samples, and M indicates the dimension of the stochastic space [43]. The MC method is not affected by the dimension of the integration domain, such as collocations on sparse or full grid methods [44,45]. A numerical comparison of other QMC sequences is presented in [46].

Construction of a cheap generalized polynomial chaos expansion-based surrogate model [47,48,49] is an alternative to the MC method. Some well-known functions, such as the multivariate Legendre, Hermite, Chebyshev, or Laguerre functions, have been taken as a basis [47,50]. Surrogate models have pros and cons. The pros are that the model can be easily sampled once constructed, and all samples are almost free (much cheaper than sampling the original stochastic PDE). For some problem settings, sampling is unnecessary because the solution can be computed analytically (e.g., computing an integral of a polynomial). The nontrivial part of surrogate models is to define how many coefficients are needed and how accurately they should be computed. Another difficulty is that not every function can be approximated well by a polynomial. The MLMC methods do not have such limitations.

This work is structured as follows. 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, computing different statistics, the performance of the MLMC method, and the practical performance of the parallel ug4 solver for the Henry problem [4,5] with uncertain coefficients. Finally, we conclude this work with a discussion in Section 5.

Our contribution: We investigate the propagation of uncertainties in the Henry-like problem. Assuming the porosity, permeability, and recharge are uncertain, we estimate the uncertainties in the density-driven flow. To reduce the high computing complexity, we applied the existing MLMC technique. We use the multigrid ug4 software library as a black-box solver, allowing us to solve the Henry problem and others (see more in [2]). We run all MLMC random simulations in parallel. To the best of our knowledge, we are unaware of any other studies where Henry’s problem [4,5] was solved using MLMC methods with uncertain porosity, permeability, and recharge parameters.

2. Henry Problem with Uncertain Porosity and Permeability

2.1. Problem setting

In coastal aquifers, salty seawater intruding on the formation on one side (the seaside) displaces the pure water due to water recharge from land sources and precipitation from the other side. Due to its higher density, seawater mainly penetrates along the bottom of the aquifer. This process can achieve a steady state but may be time-dependent due to the periodicity of the recharge or controlling the pumping rate from the wells. An accurate simulation of the salinization is vital for the prediction of water resource availability. However, the accuracy of such predictions strongly depends on the hydrogeological parameters of the formation and the geometry of the computational domain, denoted by $\mathcal{D}$.

The aquifer $\mathcal{D}\subset {\mathbb{R}}^d$, $d\in \{2,3\}$, can be modeled as an immobile porous matrix filled with liquid phase—a solution of salt in water. Due to the nonhomogeneous density distribution, gravitation induces the motion of the liquid phase. This motion transports the salt, which is otherwise subject to molecular diffusion.

A straightforward but very demonstrative model of coastal aquifers is the so-called Henry problem, first considered in [4]. In this two-dimensional setting, the aquifer is represented by a rectangular domain $\mathcal{D}=[0,2]\times [-1,0]$$\left[{\mathrm{m}}^2\right]$ entirely saturated with the liquid phase (Figure 1). The salty seawater intrudes from the right side, and pure water is recharged from the left. The top and bottom are considered impermeable. Analogous settings with partially saturated domains are considered in [51].

The mass conservation laws for the entire liquid phase and salt yield the following equations

$$\begin{array}{ccc}\hfill {\partial }_t\left(\varphi \rho \right)& +& \nabla ·\left(\rho \mathbf{q}\right)=0,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\partial }_t\left(\varphi \rho c\right)& +& \nabla ·(\rho c\mathbf{q}-\rho \mathbf{D}\nabla c)=0,\hfill \end{array}$$

(2)

where $\varphi :\mathcal{D}\to \mathbb{R}$ denotes the porosity, $\mathbf{K}:\mathcal{D}\to {\mathbb{R}}^{d\times d}$ represents the permeability, $c(t,\mathbf{x}):[0,+\infty )\times \mathcal{D}\to [0,1]$ is the mass fraction of the salt (or of the brine) in the solution, $\rho =\rho \left(c\right)$ indicates the density of the liquid phase, and $\mathbf{D}(t,\mathbf{x}):[0,+\infty )\times \mathcal{D}\to {\mathbb{R}}^{d\times d}$ denotes the molecular diffusion and mechanical dispersion tensor. For the velocity $\mathbf{q}(t,\mathbf{x}):[0,+\infty )\times \mathcal{D}\to {\mathbb{R}}^d$, we assume Darcy’s law:

$$\begin{array}{c}\hfill \mathbf{q}=-\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g}),\end{array}$$

(3)

where $p=p(t,\mathbf{x}):[0,+\infty )\times \mathcal{D}\to \mathbb{R}$ is the hydrostatic pressure, $\mu =\mu \left(c\right)$ denotes the viscosity of the liquid phase, and $\mathbf{g}={(0,\dots ,0,-g)}^T\in {\mathbb{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 [4], for this variant of the Henry problem, we set

$$\begin{array}{c}\hfill \rho \left(c\right)={\rho }_0+({\rho }_1-{\rho }_0)c,\mu =\mathrm{const}\end{array}$$

(4)

and

$$\begin{array}{c}\hfill \mathbf{D}=\varphi D\mathbf{I}\end{array}$$

(5)

with a constant scalar $D\in \mathbb{R}$, and the identity matrix $\mathbf{I}\in {\mathbb{R}}^{d\times d}$. Furthermore, we assume the isotropic permeability

$$\begin{array}{c}\hfill \mathbf{K}=K\mathbf{I},K\in \mathbb{R}.\end{array}$$

This setting is consistent with the problem setting in [3]. However, we do not assume the Boussinesq approximation and keep the density variable for all terms. For the initial conditions, we set

$${\left.c\right|}_{t=0}=0.$$

(6)

The boundary conditions are presented in Figure 1(a). On the right side of the domain, we impose Dirichlet conditions for the c and p variables that represent the adjacent seawater aquifer:

$${\left.c\right|}_{x=2}=1,{\left.p\right|}_{x=2}=-{\rho }_{1}gy.$$

(7)

On the left side, we prescribe the inflow of fresh water:

$${\left.c\right|}_{x=0}=0,{\left.\rho \mathbf{q}·{\mathbf{e}}_x\right|}_{x=0}={\widehat{q}}_{\mathrm{in}},$$

(8)

where ${\mathbf{e}}_x={(1,0)}^{\top }$, and ${\widehat{q}}_{\mathrm{in}}$ is a constant. For the classical formulation of the Henry problem, this value was set to ${\widehat{q}}_{\mathrm{in}}=6.6·{10}^{-2}$$\mathrm{kg}/\mathrm{s}$ in [3] or ${\widehat{q}}_{\mathrm{in}}=3.3·{10}^{-2}$$\mathrm{kg}/\mathrm{s}$ in [5,6]. The Neumann zero boundary conditions are imposed on the upper and lower sides of $\mathcal{D}$.

An example of $c(t,\mathbf{x})$ and $\mathbf{q}(t,\mathbf{x})$ for the parameters from Table 1 is presented in Figure 1(right). The dark red color corresponds to $c=1$, and dark blue corresponds to $c=0$. Due to its higher density, the saltwater intrudes into the aquifer in the lower right part. It is pushed back by the lighter pure water coming from the left. This process induces a vortex in the flow in the lower right corner of the domain. The saltwater flows in at the lower part of the right boundary and deviates to the top and right, back to the seaside, forming a salt triangle. This flow does not transport the salt to the left part of the domain. The salt propagates further to the left due to diffusion and dispersion and is washed out by the recharge. In the classical formulation, this salt triangle initially increases over time but achieves a steady state (cf. [3,5,6]). However, the initial nonstationary phase may take significant time. Investigating this phase is especially important to understand the system behavior when changing the recharge. For this, in addition to the mean and variance, we consider the mass fraction at 12 points (listed below) and an integral value—the total amount of pure water (as in Eq. 23). The list of chosen points follows:

$$\begin{array}{cc}\hfill \left\{{(x,y)}_{i=1,\dots ,12}\right\}& =\left\{\right(1.10,-0.95),(1.35,-0.95),(1.60,-0.95),(1.85,-0.95),(1.10,-0.75),(1.35,-0.75),\hfill \\ \hfill & (1.60,-0.75),(1.85,-0.75),(1.10,-0.50),(1.35,-0.50),(1.60,-0.50),(1.85,-0.50).\}\hfill \end{array}$$

(9)

The motivation is to consider points where the concentration variation is considerable. In addition, the mass fraction c at each point $\mathbf{x}$ is a function of time.

These spatial points may help track salinity changes over time in groundwater wells and understand which areas in the aquifer are most vulnerable. Farmers can use this information to take action, such as decreasing salinity or adapting strategies by planting salt-tolerant crops.

2.2. Modeling porosity, permeability, and mass fraction

The primary sources of uncertainty are the hydrogeological properties of the porous medium—porosity ($\varphi$) and permeability ($\mathbf{K}$) fields of the solid phase—and the freshwater recharge flux ${\widehat{q}}_x$ through the left boundary. The QoIs are related to the mass fraction c, a function of $\varphi$, $\mathbf{K}$, and the recharge. We model the uncertain $\varphi$ using a random field and assume $\mathbf{K}$ to be isotropic and dependent on $\varphi$:

$$\begin{array}{c}\hfill \mathbf{K}=K\mathbf{I},K=K\left(\varphi \right)\in \mathbb{R}.\end{array}$$

(10)

The distribution of $\varphi (\mathbf{x},\mathbf{\xi })$, $\mathbf{x}\in \mathcal{D}$, is determined by a set of stochastic parameters $\mathbf{\xi }=({\xi }_1,\dots ,{\xi }_M,...)$. Each component ${\xi }_i$ is a random variable depending on a random event $\omega$. For concision, we skip $\omega$ and write $\mathbf{\xi }:=\mathbf{\xi }\left(\omega \right)$.

The dependence in Eq. (10) is specific for every material. We refer to [52,53,54] for a detailed discussion. In the proposed model, we use a Kozeny–Carman-like dependence

$$\begin{array}{c}\hfill K\left(\varphi \right)={\kappa }_{KC}·{\displaystyle \frac{{\varphi }^3}{1-{\varphi }^2}},\end{array}$$

(11)

where the scaling factor ${\kappa }_{KC}$ is chosen to satisfy the equality $K\left(\mathbb{E}\left[\varphi \right]\right)\mathbf{I}=E\left(\mathbf{K}\right)$, resembling the parameters of the standard Henry problem. The inflow flux is kept constant across the left boundary but depends on the stochastic variable $q_{\mathrm{in}}$. We also assume that the inflow flux is independent of $\varphi$ and $\mathbf{K}$.

2.3. Numerical methods for the deterministic problem

The system (1–2) is numerically solved in the domain $\mathcal{D}\times [0,T]$, where the symbol × denotes the Cartesian product. After the discretization of $\mathcal{D}$ using quadrilaterals of size h, we obtain ${\mathcal{D}}_h$. Equations (1–2) are discretized using a vertex-centered finite-volume scheme with a “consistent velocity” for the approximation of Darcy’s law (3), as presented in [55,56,57]. The degrees of freedom associated with ${\mathcal{D}}_h$ are denoted by n. There are two degrees of freedom per grid vertex: one for the mass fraction and another for the pressure. We use the implicit Euler method with a fixed time step $\tau$ for time discretization. The number of the computed time steps is $r=T/\tau$.

We use partial upwind for the convective terms (cf. [55]). Therefore, the discretization error is of the second order w.r.t. the spatial mesh size h. However, the diffusion in (2) is minimal compared with the velocity. For the grids in the numerical experiments, the observed reduction of the discretization error after grid refinement corresponds to the first order. Thus, we assume the first-order dependence of the discretization error w.r.t. h, which is consistent with the numerical experiments. Furthermore, the Euler method provides the first-order dependence of the discretization error w.r.t. $\tau$.

The implicit time-stepping scheme provides unconditional stability but requires the solution to 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. [58]) preconditioned with the geometric multigrid method (V-cycle, cf. [59]). In the multigrid cycle, the ILU_β-smoothers [60] and Gaussian elimination are used as the coarse grid solver.

To construct the spatial grid hierarchy ${\mathcal{D}}_0,{\mathcal{D}}_1,\dots ,{\mathcal{D}}_L$, we start with a coarse grid consisting of 512 grid elements (quadrilaterals) and $n_0=1122$ degrees of freedom. This grid is regularly refined to obtain all other grid levels. After every spatial grid refinement, the number of grid elements is multiplied by a factor of four. Consequently, the number of degrees of freedom is increased by a factor of four (i.e., $n_{\ell }\approx n_0·2^{d\ell }$, $d=2$; see Table 2). This hierarchy is used in the geometric multigrid preconditioner and MLMC method. We also construct the temporal grid hierarchy ${\mathcal{T}}_0,{\mathcal{T}}_1,\dots ,{\mathcal{T}}_L$. The time step on each temporal grid is denoted by ${\tau }_{\ell }$ with ${\tau }_{\ell +1}=\frac{1}{2}{\tau }_{\ell }$. The number of time steps on the ℓth grid (level) is $r_{\ell +1}=2r_{\ell }$ and $r_{\ell }=r_0{2}^{\ell }$, where $r_0$ is the number of grid points on ${\mathcal{T}}_0$. On the ℓth level, the MLMC uses the grid ${\mathcal{D}}_{\ell }\times {\mathcal{T}}_{\ell }$. Up to six spatial and time grids were used in the numerical experiments.

In the context of this work, it is critical to estimate the numerical complexity of the deterministic solver w.r.t. the grid level ℓ. The most time-consuming part of the simulation is the solution of the discretized nonlinear system. Typically, it is challenging to predict the number of Newton iterations in every time step, but in the numerical experiments, two iterations were sufficient to achieve the prescribed accuracy. Accordingly, the linear solver was called at most two times per time step. Furthermore, the convergence rate of the geometric multigrid method does not depend on the mesh size (cf. [60]). Hence, the computation complexity of one time step is $\mathcal{O}\left(n_{\ell }\right)$, where $n_{\ell }$ is the number of the degrees of freedom on the grid level ℓ. Therefore, the overall numerical cost of the computation of one scenario on grid level ℓ for $r_{\ell }$ time steps is

$$s_{\ell }=\mathcal{O}\left(n_{\ell }{r}_{\ell }\right),s_{\ell }\propto s_{\ell -1}{2}^{(d+1)},d=2.$$

(12)

3. Multilevel Monte Carlo

Various numerical methods can quantify uncertainty, and every method has pros and cons. For example, the classical MC method converges slowly and requires numerous samples. 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. Further, we repeat the main idea of the MLMC method. A more in-depth description of these techniques is found in [61,62,63,64,65,66,67].

We let $\mathbf{\xi }\left(\omega \right)$ and $g\left(\mathbf{\xi }\right)=g\left(\mathbf{\xi }\right(\omega \left)\right)$ represent a vector of random variables and the QoI, respectively, where $\omega$ is a random event. The MLMC method aims to approximate the expected value $\mathbb{E}\left[g\right]$ with an optimal computational cost. In this work, g could be $c(t,\mathbf{x},\mathbf{\xi })$ in the whole domain or at a point $(t,\mathbf{x})$ or an integral over a subdomain. The MLMC method constructs a telescoping sum, defined over a sequence of spatial and temporal meshes, $\ell =0,\dots ,L$, as described next, to achieve this goal. Moreover, g, numerically evaluated on level ℓ, is denoted by $g_{h_{\ell },{\tau }_{\ell },\ell }$ or, for simplicity, by just $g_{\ell }$, where $h_{\ell }$ and ${\tau }_{\ell }$ are the discretization steps in space and time on level ℓ. Further, we assume that $\mathbb{E}\left[g_{h,\tau }\right]\to \mathbb{E}\left[g\right]$ as $h\to 0$ and $\tau \to 0$.

Furthermore, $s_0$ is the computing cost to evaluate one realization of $g_0$ (the most expensive one from all realizations). Similarly, $s_{\ell }$ denotes the computing cost of evaluating $g_{\ell }-g_{\ell -1}$. For simplicity, we assume that $s_{\ell }$ for $g_{\ell }-g_{\ell -1}$ is almost the same as $s_{\ell }$ for $g_{\ell }$. The number of iterations is variable; thus, the cost of computing a sample of $g_{\ell }-g_{\ell -1}$ may fluctuate for various realizations.

For a better understanding, we consider a two-level MLMC (cf. [64]) and estimate the optimal number of needed samples on both levels. The two-level MLMC has only two meshes: a coarse one and a fine one. The QoI $\mathbb{E}\left[g\right]$ can be approximated on the fine mesh by $\mathbb{E}\left[g_1\right]$ and on the coarse mesh by $\mathbb{E}\left[g_0\right]$. Furthermore,

$$\mathbb{E}\left[g_1\right]=\mathbb{E}\left[g_0\right]+\mathbb{E}\left[g_1-g_0\right]\approx m_0^{-1}\sum _{i=1}^{m_0}{g}_0^{\left(i\right)}+m_1^{-1}\sum _{j=1}^{m_1}(g_1^{\left(j\right)}-g_0^{\left(j\right)}),$$

(13)

where $g_1^{\left(j\right)}-g_0^{\left(j\right)}:=g_1\left({\mathbf{\xi }}_j\right)-g_0\left({\mathbf{\xi }}_j\right)$, ${\mathbf{\xi }}_j$ is a random vector, and $m_0$ and $m_1$ represent the numbers of quadrature points (numbers of samples/realizations) on the coarse and fine meshes, respectively. The total computational cost of evaluation (13) is $m_0{s}_0+m_1{s}_1$. The variances of $g_0$ and $g_1-g_0$ are denoted by $V_0$ and $V_1$, and the total variance is obtained by $V_0/m_0+V_1/m_1$, assuming that $g_0^{\left(i\right)}$ and $g_1^{\left(j\right)}-g_0^{\left(j\right)}$ use independent samples. By solving an auxiliary minimization problem, the variance is minimal if $m_1=m_0·\frac{\sqrt{V_1/s_1}}{\sqrt{V_0/s_0}}$. Thus, with the estimates of the variances and $m_0$, we can estimate $m_1$.

The idea presented above can be extended to a case with multiple levels. Thus, we can find (quasi-) optimal numbers of samples $m_0,m_1,\dots ,m_L$. The MLMC method calculates $\mathbb{E}\left[g_L\right]\approx \mathbb{E}\left[g\right]$ using the following telescopic sum:

$$\begin{array}{cc}\hfill \mathbb{E}\left[g_L\right]& =\mathbb{E}\left[g_0\right]+\sum _{\ell =1}^L\mathbb{E}\left[g_{\ell }-g_{\ell -1}\right]\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & \approx m_0^{-1}\sum _{i=1}^{m_0}{g}_0^{(0,i)}+\sum _{\ell =1}^L\left(m_{\ell }^{-1}\sum _{i=1}^{m_{\ell }}(g_{\ell }^{(\ell ,i)}-g_{\ell -1}^{(\ell ,i)})\right).\hfill \end{array}$$

(15)

In the above equation, level ℓ in the superscript $(\ell ,i)$ indicates that independent samples are used at each correction level. As ℓ increases, the variance of $g_{\ell }-g_{\ell -1}$ decreases. Thus, the total computational cost can be reduced by taking fewer samples on finer meshes.

We recall that $h_{\ell }=h_0·2^{-2\ell }$ and ${\tau }_{\ell }={\tau }_0·2^{-\ell }$. We assume that the average cost of generating one sample of $g_{\ell }$ (the cost of one deterministic simulation for one random realization) is

$$s_{\ell }=\mathcal{O}\left(n_{\ell }{r}_{\ell }\right)=\mathcal{O}\left(h_{\ell }^{-1}{\tau }_{\ell }^{-1}\right)=\mathcal{O}\left(\frac{1}{h_0{\tau }_0}{2}^{2\ell }{2}^{\ell }\right)=\mathcal{O}\left(\frac{1}{h_0{\tau }_0}{2}^{3\ell }\right)=\mathcal{O}\left(\frac{1}{h_0{\tau }_0}{2}^{(d+1)\ell \gamma }\right),$$

(16)

where $d=2$ is the spatial dimension, and $\gamma =1$ is determined by the computational complexity of the deterministic solver (ug4).

We let $V_{\ell }$ be the variance of one sample of $g_{\ell }-g_{\ell -1}$. Then, the total cost and variance of the multilevel estimator in Eq. (14) are ${\sum }_{\ell =0}^L{m}_{\ell }{s}_{\ell }$ and ${\sum }_{\ell =0}^L{V}_{\ell }/m_{\ell }$, respectively. For a fixed variance, the cost is minimized by choosing $m_{\ell }$ to minimize the following functional for some value of the Lagrange multiplier ${\mu }^2$:

$$F(m_0,\dots ,m_L):=\sum _{\ell =0}^L{m}_{\ell }{s}_{\ell }+{\mu }^2\frac{V_{\ell }}{m_{\ell }}.$$

(17)

To determine $m_{\ell }$, we take the derivatives w.r.t. $m_{\ell }$ and set them equal to zero:

$$\frac{\partial F(m_0,\dots ,m_L)}{\partial m_{\ell }}:=s_{\ell }-{\mu }^2\frac{V_{\ell }}{m_{\ell }^2}=0.$$

After solving the obtained equations, we obtain

$$m_{\ell }^2={\mu }^2\frac{V_{\ell }}{s_{\ell }}\mathrm{and}{m}_{\ell }=\mu \sqrt{\frac{V_{\ell }}{s_{\ell }}}.$$

To achieve an overall variance of ${\epsilon }^2$, that is,

$$\sum _{\ell =0}^L{V}_{\ell }/m_{\ell }={\epsilon }^2,$$

we substitute $m_{\ell }$ with the computed $m_{\ell }=\mu \sqrt{\frac{V_{\ell }}{s_{\ell }}}$, and obtain

$$\sum _{\ell =0}^L\frac{V_{\ell }}{\mu \sqrt{\frac{V_{\ell }}{s_{\ell }}}}={\epsilon }^2.$$

From the last equation, we obtain

$$\mu ={\epsilon }^{-2}\sum _{\ell =0}^L\sqrt{V_{\ell }{s}_{\ell }},\mathrm{and}$$

$$m_{\ell }={\epsilon }^{-2}\sqrt{\frac{V_{\ell }}{s_{\ell }}}\sum _{i=0}^L\sqrt{V_i{s}_i}.$$

(18)

The total computational cost is $S:={\epsilon }^{-2}{\left({\sum }_{\ell =0}^L\sqrt{V_{\ell }{s}_{\ell }}\right)}^2$ (for further analysis of this sum, see [64], p.4).

Definition 1. 

We let

$$\begin{array}{c}\hfill \mathbb{E}\left[Y_{\ell }\right]:=\left\{\begin{array}{cc}\mathbb{E}\left[g_0\right],\hfill & \ell =0\hfill \\ \mathbb{E}\left[g_{\ell }-g_{\ell -1}\right],\hfill & \ell >0\hfill \end{array}\right..\end{array}$$

(19)

In addition, $Y:={\sum }_{\ell =0}^L{Y}_{\ell }$ denotes a multilevel estimator of $\mathbb{E}\left[g\right]$ based on $L+1$ levels and $m_{\ell }$ independent samples on level ℓ, where $\ell =0,\dots ,L$. Moreover, $Y_{\ell }=m_{\ell }^{-1}{\sum }_{i=1}^{m_{\ell }}(g_{\ell }^{(\ell ,i)}-g_{\ell -1}^{(\ell ,i)})$, where $g_{-1}\equiv 0$.

The standard theory indicates that $\mathbb{E}\left[Y\right]=\mathbb{E}\left[g_L\right]$, $\mathbb{V}\left[Y_{\ell }\right]={\sum }_{\ell =0}^L{m}_{\ell }^{-1}{V}_{\ell }$, and $V_{\ell }\equiv \mathbb{V}\left[g_{\ell }-g_{\ell -1}\right]$.

The mean squared error (MSE) is used to measure the quality of the multilevel estimator:

$$\mathrm{MSE}:=\mathbb{E}\left[{(Y-\mathbb{E}\left[g\right])}^2\right]=\mathbb{V}\left[Y\right]+{\left(\mathbb{E}\left[Y\right]-\mathbb{E}\left[g\right]\right)}^2.$$

(20)

To obtain an MSE smaller than ${\epsilon }^2$, we ensure that both $\mathbb{V}\left[Y\right]$ and ${\left(\mathbb{E}\left[Y\right]-\mathbb{E}\left[g\right]\right)}^2={\left(\mathbb{E}\left[g_L-g\right]\right)}^2$ are smaller than ${\epsilon }^2/2$. Combining this idea with a geometric sequence of levels in which the cost increases exponentially with the level while the weak error $\mathbb{E}\left[g_L-g\right]$ and multilevel correction variance $V_{\ell }$ decrease exponentially leads to the following theorem (cf. Theorem 1, p. 6 in [64]):

Theorem 2. 

We let d denote the problem dimension. Suppose positive constants $\alpha ,\beta ,\gamma >0$ exist such that $\alpha \ge \frac{1}{2}min(\beta ,\gamma d)$, and

$$\begin{array}{cc}\hfill |\mathbb{E}\left[g_{\ell }-g\right]|& \le c_1{2}^{-\alpha \ell }\hfill \end{array}$$

(21a)

$$\begin{array}{cc}\hfill V_{\ell }& \le c_2{2}^{-\beta \ell }\hfill \end{array}$$

(21b)

$$\begin{array}{cc}\hfill S_{\ell }& \le c_3{2}^{d\gamma \ell }.\hfill \end{array}$$

(21c)

Then, for any accuracy $\epsilon <e^{-1}$, a constant $c_4>0$ and a sequence of realizations ${\left\{m_{\ell }\right\}}_{\ell =0}^L$ exist, such that

$$\mathrm{MSE}:=\mathbb{E}\left[{(Y-\mathbb{E}\left[g\right])}^2\right]<{\epsilon }^2,$$

and the computational cost is

$$\begin{array}{c}\hfill S=\left\{\begin{array}{cc}\mathcal{O}\left({\epsilon }^{-2}\right),\hfill & \beta >d\gamma \hfill \\ \mathcal{O}({\epsilon }^{-2}){\left(\log \left(\epsilon \right)\right)}^2,\hfill & \beta =d\gamma \hfill \\ \mathcal{O}\left({\epsilon }^{-\left(2+\frac{d\gamma -\beta }{\alpha }\right)}\right),\hfill & \beta <d\gamma .\hfill \end{array}\right.\end{array}$$

(22)

This theorem (see also [61,63,68,69,70]) 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 $\mathcal{O}\left({\epsilon }^{-2-d\gamma /\alpha }\right)$. Furthermore, in the best-case scenario presented above, the computational cost of the MLMC algorithm scales as $\mathcal{O}\left({\epsilon }^{-2}\right)$.

Using preliminary tests, we can estimate the convergence rates $\alpha$ for the mean (the so-called weak convergence) and $\beta$ for the variance (the so-called strong convergence). In addition, $\alpha$ is strongly connected to the order of the discretization error (see Section 2.3), which equals 1, and precise estimates of parameters $\alpha$ and $\beta$ are crucial to distribute the computational effort optimally.

4. Numerical Experiments

The goal is to reduce the total computational cost of stochastic simulations. We use the MLMC method to compute the mean value of various QoIs, such as c in the whole domain, c at a point, or an integral value (we call it the freshwater integral):

$$Q_{FW}(t,\omega ):={\int }_{\mathbf{x}\in \mathcal{D}}I(c(t,\mathbf{x},\omega )\le 0.012178)d\mathbf{x},$$

(23)

where $I(·)$ is the indicator function identifying a subdomain $\{\mathbf{x}:c(t,\mathbf{x},\omega )\le 0.012178\}$, meaning the mass of the fresh water at a time t. Each simulation may contain up to $n=0.5·{10}^6$ spatial mesh points and a few thousand time steps ($r=6016$ on the finest mesh).

Software and parallelization: The computations presented in this work were performed using the ug4 simulation software toolbox (https://github.com/ug4/ughub.wiki.git) [71,72]. This software has been applied for subsurface flow simulations of real-world aquifers (cf. [2]). 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 [48,49]. Simulations were performed on different meshes; thus, the computation time of each simulation varied over a wide range (see Table 2).

Porosity and recharge: We assume two horizontal layers: $y\in (-0.75,0]$ (the upper layer) and $y\in [-1,-0.75]$ (the lower layer). The porosity inside each layer is uncertain and is modeled as in Eq. (24):

$$\begin{array}{cc}\hfill \varphi (\mathbf{x},\mathbf{\xi })& =0.35·(1+0.15({\xi }_2\cos (\pi x/2)+{\xi }_2\sin \left(2\pi y\right)+{\xi }_1\cos \left(2\pi x\right)))·C_0\left({\xi }_1\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \mathrm{where}{C}_0\left({\xi }_1\right)& =\left\{\begin{array}{cc}1+0.2{\xi }_1\hfill & \mathrm{if}y<-0.75\hfill \\ 1-0.2{\xi }_1\hfill & \mathrm{if}y\ge -0.75,\hfill \end{array}\right.\hfill \end{array}$$

(25)

Additionally, the recharge flux is also uncertain and is equal to

$$\begin{array}{c}\hfill {\widehat{q}}_{\mathrm{i}n}=-6.6·{10}^{-2}(1+0.5·{\xi }_3),\end{array}$$

(26)

where ${\xi }_1$, ${\xi }_2$, and ${\xi }_3$ are sampled independently and uniformly in $[-1,1]$. Figure 2 depicts a random realization of the porosity random field $\varphi \left(\mathbf{\xi }\right)$ (left) and the corresponding solution $c(t,\mathbf{x},\mathbf{\xi })=c(t,\varphi (\mathbf{\xi }\left)\right)$ at $t=T$ (right). Additionally, four isolines $\{\mathbf{x}:|c(t,\varphi \left(\mathbf{\xi }\right))-\overline{c}\left(t\right)|=0.1·i\}$, $i=1,2,3,4$, are presented on the right. The isolines demonstrate the absolute value of the difference between the computed realization $c(t,\varphi (\mathbf{\xi }\left)\right)$ and the expected value $\overline{c}\left(t\right)$. These computations were performed for $\mathbf{\xi }={\mathbf{\xi }}^{*}=(-0.5898,-0.7257,-0.9616)$ and $t=T=6016$ s.

The mean and variance of the mass fraction are provided in Figure 3 on the left and right, respectively. The expectation takes values from $[0,1]$, and the variance range is $[0,0.05]$. 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 additional sensors). Additionally, the right image displays five contour lines $\{\mathbf{x}:\mathrm{Var}-2mu\left[c\right](t,\mathbf{x})=0.01·i\}$, $i=1..5$, $t=T=6016$.

We observed that the variability (uncertainty) of the mass fraction might 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 very low or high (dark red regions). In regions with high uncertainty, refining the mesh and applying the MLMC method make sense.

Before we run the MLMC method, we first examine the solution $c(t,\mathbf{x})$ at 12 preselected points (see Eq. (9)). Figure 4 includes 12 subfigures. Each subfigure presents 600 QMC realizations of $c(t,\mathbf{x})$ and five quantiles depicted by dotted lines. The dotted line at the bottom indicates the quantile $0.025$. The following dotted line is the quantile $0.25$, and the dotted line on the top indicates the quantile $0.975$. All five quantiles from the bottom to the top are 0.025, 0.25, 0.50, 0.75, and 0.975, respectively. We observe that c at the final point $t=T$ varies considerably.

Example. In Figure 5, we demonstrate the probability density function (pdf) of $t^{*}\left(\omega \right)={min}_t\{t:Q_{FW}(t,\omega )<1.2\}$ (left), and the pdf of $t^{*}\left(\omega \right)={min}_t\{t:Q_{FW}(t,\omega )<1.7\}$ (right). On average, after approximately 29 time steps (on the left) and six time steps (on the right), the volume of the fresh water becomes less than 1.2 and 1.7, respectively. The initial volume of the fresh water was 2.0.

All 600 realiations of $Q_{FW}\left(t\right)$ are depicted in Figure 6 The time is along the x-axis, $t\in [\tau ,48\tau ]$. Additionally, five quantiles are represented by dotted curves from the bottom to the top and are 0.025, 0.25, 0.50, 0.75, and 0.975, respectively.

Example. Figure 7 (left) displays the evolution of the pdf of $c(t,\mathbf{x},\omega )$ at a fixed point $\mathbf{x}=(1.85,-0.95)$ in time $t=\{3\tau ,\dots ,48\tau \}$. From left to right, the farthest left (blue) pdf corresponds to $t=3\tau$, the second curve from the left (red) corresponds to $t=4\tau$, and so on. In the beginning, $t=3\tau$, and the mass fraction c is low, about 0.15 on average. Then, with time, c increases and, at $t=T=48\tau$, is approximately equal to 1. Example. The next QoI is the earliest time moment when $c(t,\mathbf{x})$, at fixed $\mathbf{x}=(1.85,-0.95)$, becomes smaller than the threshold value $0.9$ (maximum is 1.0). Figure 7 (right) presents its pdf. On average, after $t\approx 10$ time steps, the mass fraction becomes smaller than $0.9$, but 40 time steps are needed in some scenarios.

Next, we research how $g_{\ell }-g_{\ell -1}$ depends on the time and level. All graphics in Figure 8 display 100 realizations of the differences between solutions computed on two neighbor meshes for every time point $t_i$, $i=1\dots 48$ (along the x-axis). The top left graphic indicates the differences between the mass fractions computed on Levels 1 and 0. The other graphics reveal the same, but for Levels 2 and 1, 3 and 2, 4 and 3, and 5 and 4, respectively. The largest value decreases from $2.5·{10}^{-2}$ (top left) to $5·{10}^{-4}$. Considerable variability is observed for $t\in [3,7]$ and $t\in [8,25]$. Starting with $t\approx 30$, the variability between solutions decreases and stabilizes. From these five graphics, we can estimate that the maximal amplitude decreases by a factor $\approx 2$, at 0.015, 0.008, 0.004, 0.0015, and 0.0008. However, it is challenging to make a similar statement about each time point t. This observation makes it difficult to estimate the weak and strong convergence rates and the optimal number of samples correctly on each mesh level. They are different for each time t (and for each $\mathbf{x}$). For some time points, the solution is smooth and requires only a few levels and a few samples on each level. For other points with substantial dynamics, the numbers of levels and samples are higher.

Because $g_{\ell }-g_{\ell -1}$ is random, we visualize its mean and variance. Figure 9 demonstrates the mean (left) and variance (right) of the differences in concentrations $g_{\ell }-g_{\ell -1}$, $\ell =1,\dots ,5$. On the left, the amplitude decreases when ℓ increases. A slight exception is the blue line for $t\approx 9,10,11$ (right). A possible explanation is that the solutions $g_0$ or $g_1$ are insufficiently accurate. The right image presents how the amplitude of the variances decays. This decay is necessary for the successful work of the MLMC method. We also observe a potential issue; the weak and strong convergence rates vary for various time points t. Thus, determining the optimal number of samples $m_{\ell }$ for each level is not possible (only suboptimal).

At the beginning $t=0$, the variability is zero and starts to increase. We observe changes during a specific time interval, and then the process starts to stabilize after $\approx 45$ time steps. The variability is either unchanging from level to level or decreases.

Table 2 contains average computing times, which are necessary to estimate the number of samples $m_{\ell }$ at each level ℓ. The fourth column contains the average computing time, and the fifth and sixth 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 $\approx 6016$ s, the solution is almost unchanging; thus, we restrict this to only $t\in [0,T]$, where $T=6016$. For example, if the number of time steps is $r_{\ell }=188$ (Level 0 in Table 2), then the time step $\tau =\frac{T}{r_{\ell }}=\frac{6016}{188}=32$ s.

The time step $\tau$ is adaptive and changing from $\tau =\frac{6016}{128}=32$ s (very coarse mesh) to $\tau =\frac{6016}{6016}=1$ s (finest mesh). Starting with level $\ell =2$, the average time increases by a factor of eight. These numerical tests confirm the theory in Eq. (12), stating that the numerical solver is linear w.r.t. $n_{\ell }$ and $r_{\ell }$.

Table 2. Number of degrees of freedom $n_{\ell }$, number of time steps $r_{\ell }$, step size in time ${\tau }_{\ell }$, average, minimal, and maximal computing times on each level ℓ.

Level ℓ	${\mathit{n}}_{\mathit{\ell }}$	${\mathit{r}}_{\mathit{\ell }}$	${\mathit{\tau }}_{\mathit{\ell }}=6016/{\mathit{r}}_{\mathit{\ell }}$	Computing times (${\mathit{s}}_{\mathit{\ell }}$)
				average	min.	max.
0	1122	188	32	1.15	0.88	1.33
1	4290	376	16	4.1	3.4	4.87
2	16770	752	8	19.6	17.6	22
3	66306	1504	4	136.0	128	144
4	263682	3008	2	1004.0	891	1032
5	1051650	6016	1	8138.0	6430	8480

Table 2. Number of degrees of freedom $n_{\ell }$, number of time steps $r_{\ell }$, step size in time ${\tau }_{\ell }$, average, minimal, and maximal computing times on each level ℓ.

Level ℓ	${\mathit{n}}_{\mathit{\ell }}$	${\mathit{r}}_{\mathit{\ell }}$	${\mathit{\tau }}_{\mathit{\ell }}=6016/{\mathit{r}}_{\mathit{\ell }}$	Computing times (${\mathit{s}}_{\mathit{\ell }}$)
				average	min.	max.
0	1122	188	32	1.15	0.88	1.33
1	4290	376	16	4.1	3.4	4.87
2	16770	752	8	19.6	17.6	22
3	66306	1504	4	136.0	128	144
4	263682	3008	2	1004.0	891	1032
5	1051650	6016	1	8138.0	6430	8480

With estimates for each level, we can estimate the rates of $\alpha$ and $\beta$ (Eqs. (21a)-(21b)) in weak and strong convergences.

The slope in Figure 10 can be used to estimate the rates of the weak (left) and strong (right) convergences. The differences are indicated on the horizontal axis.

We use computed variances $V_{\ell }$ and computing times (work) $s_{\ell }$ from Table 2 to estimate the optimal number of samples $m_{\ell }$ and compute the telescopic sum from Eq. (15) to approximate the expectation.

Table 3 lists $m_{\ell }$ for a given total variance ${\epsilon }^2$:

After the telescopic sum is computed, we can compare the results with the QMC results. Figure 11 depicts the decay of the absolute (left) and relative (right) errors vs. levels along the x-axes. The ’true’ solution was computed using the QMC method on the finest mesh level $L=5$.

5. Conclusion

We investigated the applicability and efficiency of the MLMC approach for the Henry-like problem with uncertain porosity, permeability, and recharge. These uncertain parameters were modeled by random fields with three independent random variables. The numerical solution for each random realization was obtained using the well-known ug4 parallel multigrid solver. The number of required random samples on each level was estimated by computing the decay of the variances and computational costs for each level. These estimates depend on the minimization function in Eq. (17).

We also computed the expected value and variance of the mass fraction in the whole domain, the evolution of the pdfs, the solutions at a few preselected points $(t,\mathbf{x})$, and the time evolution of the freshwater integral value. We have found that some QoIs require only 2-3 of the coarsest mesh levels, and samples from finer meshes would not significantly improve the result. Note that a different type of porosity in Eq. (24) may lead to a different conclusion.

The results show that the MLMC method is faster than the QMC method at the finest mesh. Thus, sampling at different mesh levels makes sense and helps to reduce the overall computational cost.

• Limitations. 1. It may happen that the QoIs computed on different grid levels are the same (for the given random input parameters). In this case the standard (Q)MC on a coarse mesh will be sufficient. 2. The time dependence is challenging. The optimal number of samples depends on the point $(t,\mathbf{x})$ and may be small for some points and large for others. 3. Twenty-four hours may not be sufficient to compute the solution at the sixth mesh level.
• Future work. Our model of porosity in Eq. (24) is quite simple. It would be beneficial to consider a more complicated/multiscale/realistic model with more random variables. 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_{\ell }$. Another hot topic is data assimilation and the identification of unknown parameters [73,74,75,76]. Known experimental data and measurements of porosity, permeability, velocity or mass fraction could be used to minimise uncertainties.