Multi-grid Multi-Level Monte Carlo Method for Stokes-Darcy interface Model with Random Hydraulic ConductivityThis work is partially supported by NSF grants DMS-1418624 and DMS-1722647, National Science Foundation of China (91330104).

# Multi-grid Multi-Level Monte Carlo Method for Stokes-Darcy interface Model with Random Hydraulic Conductivity††thanks: This work is partially supported by NSF grants DMS-1418624 and DMS-1722647, National Science Foundation of China (91330104).

Zhipeng Yang Division of Applied and Computational Mathematics, Beijing Computational Science Research Center, Beijing 100094, P. R. China, yangzhp@csrc.ac.cn    Ju Ming School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, P. R. China, jming@hust.edu.cn, corresponding author.    Xiaoming He Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409, U. S. A., hex@mst.edu    Li Zhang School of mathematical Sciences, University of Electronic Science and Technology, Chengdu Sichuan 611731, P. R. China, lizhang_137363@163.com
###### Abstract

In this article we develop a multi-grid multi-level Monte Carlo (MGMLMC) method for the stochastic Stokes-Darcy interface model with random hydraulic conductivity both in the porous media domain and on the interface. Because the randomness through the interface affects the flow in the Stokes domain, we investigate the coupled stochastic Stokes-Darcy model to improve the fidelity as this model also considers the second and third porosity of the free flow. Then we prove the existence and uniqueness of the weak solution of the variational form. For the numerical solution, we adopt the Monte Carlo (MC) method and finite element method (FEM), for the discrete form in the probability space and physical space, respectively. In the traditional single-level Monte Carlo (SLMC) method, more accurate numerical approximate requires both larger number of samples in probability space and smaller mesh size in the physical space. Then the computational cost increase significantly, which is the product of the number of samples and the computational cost of each sample, as the mesh size becomes smaller for the more accurate numerical approximate. Therefore we adopt the multi-level Monte Carlo (MLMC) method to dramatically reduce the computational cost in the probability space, because the number of samples decays fast while the mesh size decreases. We also develop a strategy to calculate the number of samples needed in MLMC method for the stochastic Stokes-Darcy model. Furthermore MLMC naturally provides the hierarchial grids and sufficient information on these grids for multi-grid (MG) method, which can in turn improve the efficiency of MLMC. In order to fully make use of the dynamical interaction between this two methods, we propose the multi-grid multi-level Monte Carlo method for more efficiently solving the stochastic model, with additional efforts on the interface. Numerical examples are provided to verify and illustrate the proposed method and the theoretical conclusions.

Key words. stochastic Stokes-Darcy interface model, multi-level Monte Carlo, multi-grid method.

AMS subject classifications. 35R60, 65C05, 65M60 76S05.

## 1 Introduction

The Stokes-Darcy interface model has attracted significant attention from scientists and engineers due to its wide range of applications, such as interaction between surface and subsurface flows [18, 28, 29, 54, 61], industrial filtrations [34, 52], groundwater system in karst aquifers [17, 39, 50, 51], petroleum extraction [1, 3, 55], and many others [19, 21, 26, 27, 76, 81, 86]. Therefore it is not surprising that many different numerical methods have been proposed and analyzed for the Stokes-Darcy model, including domain decomposition methods [10, 15, 22, 30, 31, 49, 85], Lagrange multiplier methods [4, 40, 41, 56, 62], discontinuous Galerkin methods [24, 46, 57, 64, 73, 74], multi-grid methods [2, 12, 67], partitioned time stepping methods [59, 68, 79], coupled finite element methods [13, 58, 66, 77], and many others [5, 9, 23, 33, 37, 47, 56, 69, 78, 84, 89].

The above existing works only consider the deterministic Stokes-Darcy model, for which the problem data, including the model coefficients, the forcing terms, the domain geometry, the boundary conditions and the initial conditions, are assumed to be perfectly known. However, in reality there is a significant amount of uncertainty involved in determining these real-life data due to measurements and simplifications [25, 42, 71, 80].

There are some works on the uncertainties of the porous media flow by assuming the hydraulic conductivity of the porous media is a random field in the second order elliptic equation [32, 38, 63, 90]. But the Stokes-Darcy model has a much more complicated system for the uncertainties due to the flow interaction on the interface between the porous media flow and the free flow in conduits. Hence it is not trivial to study the effect of randomness of the hydraulic conductivity on the whole coupled flow performance, which is key component of this paper, especially around the interface.

On the other hand, in the numerical simulation area, the Monte Carlo method  has been a widely applied to solve the stochastic problems. The convergence of the Monte Carlo method is based on the number of the samplings. Unfortunately, for a high accuracy result, one usually needs a large number of samples, which significantly increases the computational cost. To develop an accurate and efficient numerical method for simulating the coupled stochastic porous media flow and free flow, we develop a multi-level Monte Carlo method [7, 20, 43, 44, 60, 83] to solve the sophisticated stochastic Stokes-Darcy interface model. This method is much more costly efficient by significantly reducing the number of samples on the fine meshes. But it is not trivial to determine how many samples should be used in each level to keep the global accuracy while minimizing the cost. Therefore, we develop a strategy based on a detailed analysis to overcome this difficulty.

Furthermore, the multi-level Monte Carlo method only reduces the computational cost in the probability space, not in the physical space. Inspired by a fact that the multi-level Monte Carlo method already has a set of hierarchical grids for the multi-level idea, it is a natural idea to fully make use of the same set of hierarchical grids to solve the discrete algebraic system by using the powerful multi-grid method [11, 65, 70, 87, 88], which can further improve the efficiency of the proposed multi-level Monte Carlo method. Meanwhile, the saved information of the multi-level Monte Carlo method on the set of hierarchical grids will also significantly reduce the computational cost of the multi-grid method. Therefore, we combine the multi-level Monte Carlo method and the multi-grid method on the same set of hierarchical grids to propose an even more costly efficient method, which is the multi-grid multi-level Monte Carlo method.

The rest of the paper is organized as follows. In section 2, we briefly recall the deterministic Stokes-Darcy model. In section 3, we present the stochastic Stokes-Darcy interface model, the weak formulation of the stochastic Stokes-Darcy model and the proof of the well-posedness. In section 4, we recall the Monte Carlo method to approximate the numerical moments of the stochastic solutions, adopt the multi-level Monte Carlo method to reduce the computational cost in probability space, and then develop the multi-grid multi-level Monte Carlo to further reduce the computational cost. In section 5 we provide numerical examples to verify the theoretical analysis and illustrate the features of the proposed methods.

## 2 Deterministic model for coupled fluid flow with porous media flow

The coupled Stokes-Darcy system describes the the free flow by Stokes equations in the conduit domain and the confined flow by Darcy system in the porous media domain. And three interface conditions displaced follow are used to couple the flows in these two domains. In this paper, we consider the coupled Stokes-Darcy system on a bounded domain , , where is the porous media domain and is the conduit domain. We decompose the boundary into two parts: , , and denote the interface as .

In the porous media domain , the flow is governed by the Darcy system 

 →um(x) = −K(x)∇ϕm(x)in Dm, \hb@xt@.01(2.1) ∇⋅→um(x) = fm(x)in Dm, \hb@xt@.01(2.2)

here, denotes the specific discharge in the porous media, is the hydraulic conductivity tensor of the porous media that is symmetric and positive definite in accordance with physical meaning, is the hydraulic head, and is the sink/source term.

By substituting (LABEL:Darcy_law_1) into (LABEL:Darcy_law_2), we obtain the second-order form of the Darcy system

 −∇⋅(K(x)∇ϕm(x))=fm(x)% in Dm. \hb@xt@.01(2.3)

In the conduit domain , the flow is governed by the Stokes equations:

 −∇⋅T(→us,ps) = →fsin Ds, \hb@xt@.01(2.4) ∇⋅→us = 0in Ds, \hb@xt@.01(2.5)

where denotes the fluid velocity, is the kinematic pressure, and is the external body force. is the stress tensor, defined as , where is the kinematic viscosity of the fluid and .

On the interface between the conduit and the porous media domain, we impose three interface conditions:

 →us⋅→ns = (K∇ϕm)⋅→nm%on ΓI, \hb@xt@.01(2.6) −→nTsT(→us,ps)→ns = g(ϕm−z)on ΓI, \hb@xt@.01(2.7) −τTjT(→us,ps)→ns = αν√d√trace(Π(x))τTj(→us+K∇ϕm)on%  ΓI, \hb@xt@.01(2.8)

where , denote the unit outer normal to the conduit and the porous media regions at the interface , respectively, denote mutually orthogonal unit tangential vectors to the interface , is the hight, is the gravitational acceleration, and is the intrinsic permeability. The first interface condition (LABEL:BJ_1) is governed by the conservation of mass, the second interface condition (LABEL:BJ_2) represents the balance of the kinematic pressure in the matrix and the stress in the free flow at the normal direction along the interface, and the last interface condition (LABEL:BJ_3) is the famous Beavers-Joseph condition [8, 14, 16, 17, 36, 53].

## 3 Stokes-Darcy interface model with random permeability

To overcome the difficulty of measuring the exact permeability at every point in the porous media domain, we use an underlying random field to describe the intrinsic permeability tensor . Thus the hydraulic conductivity tensor is also a random field with the relationship . Then we obtain the stochastic partial differential equations to describe the coupled system with the random hydraulic conductivity, based on the deterministic model in the above section. We investigate the uncertainty in the porous domain and the uncertainty transferred to the conduit domain through the interface. Furthermore, we provide the weak formulation and prove the well-posedness of the weak solution of the coupled stochastic model.

### 3.1 Functional spaces and notations

Before the study of the stochastic coupled problem, we introduce some notations. Throughout this paper, we adopt the notations in  for the classical Sobolev spaces. Let be an open, connected, bounded, and convex subset of , , with polygonal and Lipschitz continuous boundary . Let , , and be a Sobolev space on with the standard norm and semi-norm .

Let be a complete probability space. Here is the set of outcomes, is the -algebra of events, and is a probability measure.

For the given probability space and the Sobolev space with the inner product and norm , we define the stochastic Sobolev space, which consists of strongly measurable, -summable mappings , by

 L2(Ω;Wr,q(D)):={ϕ:Ω→Wr,q(D) | ϕ strongly measurable, ∥ϕ∥L2(Ω;Wr,q(D))<∞}.

Here is the norm given as, ,

 ∥ϕ∥L2(Ω;Wr,q(D)):=(E[∥ϕ(ω,⋅)∥2Wr,q(D)])1/2:=(∫Ω∥ϕ(ω,⋅)∥2Wr,q(D)dP(ω))1/2,

which is induced by following inner product, ,

 [ϕ,ψ]L2(Ω;Wr,q(D)):=E[(ϕ,ψ)Wr,q(D)]:=∫Ω(ϕ,ψ)Wr,q(D)dP(ω).

For , we denote the Hilbert space and with the standard norm and semi-norm . For , we denote with the standard norm . For , we denote and . For the vector , , 2-norm of is

For simplicity, we define

 Lq(D) = L2(Ω;Lq(D)),\ with norm \ ∥⋅∥L(D)=∥⋅∥L2(Ω;Lq(D)), Hr(D) = L2(Ω;Hr(D)),\ with norm \ ∥⋅∥Hr(D)=∥⋅∥L2(Ω;Hr(D)), →Hr(D) = L2(Ω;Hr(D)),\ with norm \ ∥⋅∥→Hr(D)=∥⋅∥L2(Ω;Hr(D)).

### 3.2 Stochastic Stokes-Darcy interface equations

With the complete probability space , let , , be a random hydraulic conductivity tensor.

Then in the porous media domain, the stochastic second-order form of Darcy equation with sink/source term is given as:

 −∇⋅(K(ω,x)∇ϕm(ω,x))=fm(x),in Dm. \hb@xt@.01(3.1)

And the interface conditions are modified as:

 →us(ω,x)⋅→ns(x) = (K(ω,x)∇ϕm(ω,x))⋅→nm(x),on ΓI, \hb@xt@.01(3.2) −→n⊤sT(→us,ps)→ns = g(ϕm(ω,x)−z),on ΓI, \hb@xt@.01(3.3) −τ⊤jT(→us,ps)→ns = αν√d√trace(Π(ω,x))τ⊤j(→us(ω,x)+K(ω,x)∇ϕm(ω,x)),on ΓI. \hb@xt@.01(3.4)

Due to the randomness transferred from porous media domain through the interface conditions, the Stokes equations in the conduit domain become stochastic and are given as follows

 −∇⋅T(→us(ω,x),ps(ω,x)) = →fs(x),in Ds, \hb@xt@.01(3.5) ∇⋅→us(ω,x) = 0,in Ds. \hb@xt@.01(3.6)

For the boundary conditions, we assume the hydraulic head and the fluid velocity satisfy homogeneous Dirichlet boundary condition except on .

### 3.3 Weak formulation of the coupled problem

We denote the velocity-pressure spaces on the conduit domain as

 Xs = {→us∈→H1(Ds) | →us=0 on Γs}, X0s = {→us∈→H0(Ds) | →us=0 on Γs}, Xs,div = {→us∈Xs | ∇⋅→us=0 in Ds}, Qs = {qs∈L2(Ds)},

and we denote the pressure space on the porous media as

 Xm={ϕm∈H1(Dm) | ϕm=0 on Γm},X0m={ϕm∈H0(Dm) | ϕm=0 on Γm}.

For convenience, let , , , and , where , . The norms of are given as

 ∥u––∥Xr=(E[∥u––∥2Hr(Ds)×Hr(Dm)])1/2=(∥→us∥2→Hr(Ds)+∥ϕm∥2Hr(Dm))1/2,r=0, 1. \hb@xt@.01(3.7)

The projection onto the local tangential plane of the vector is denoted as . Then using the boundary conditions (LABEL:stochastic_stokes_Darcy_5)-(LABEL:stochastic_stokes_Darcy_7), we obtain the following weak formulation: find , such that

 {A(u––,v–)−B(v–,ps)=F(v–),∀v–=(→vs,ψm)∈X,B(u––,qs)=0,∀qs∈Qs, \hb@xt@.01(3.8)

where

 A(u––,v–) = E[a(u––,v–)]=∫Ωa(u––,v–)dω, \hb@xt@.01(3.9) a(u––,v–) = ∫Ds2νD(→us):D(→vs)dx+g∫Dm(K∇ϕm)⋅∇ψmdx \hb@xt@.01(3.10) + g∫ΓIϕm→vs⋅→nsdΓI+∫ΓIαν√d√trace(Π)Pτ(→us)⋅→vsdΓI \hb@xt@.01(3.11) − g∫ΓI(→us⋅→ns)ψmdΓI+∫ΓIαν√d√trace(Π)Pτ(K∇ϕm)⋅→vsdΓI, \hb@xt@.01(3.12) B(v–,ps) = E[b(v–,ps)]=∫Ωb(v–,ps)dω, \hb@xt@.01(3.13) b(v–,ps) = ∫Dsps∇⋅→vsdx, \hb@xt@.01(3.14) F(v–) = E[f(v–)]=∫Ωf(–v)dω, \hb@xt@.01(3.15) f(v–) = ∫Ds→fs⋅→vsdx+g∫Dmfmψmdx+∫ΓIgz→vs⋅→nsdΓI. \hb@xt@.01(3.16)

### 3.4 Well-posedness of the weak solution

The approach to analyze the well-posedness in our paper is inspired by the ideas in [7, 17, 45, 72]. One of the following two assumptions is needed to ensure the existence and uniqueness of the weak solution.

###### Assumption 3.1

Let be a diagonal matrix as .

• the strong elliptic condition: there are positive lower and upper bounds , such that

 0
• the integrability condition: let and satisfy

 0

Under the above two assumptions, we derive some properties of the weak formulation.

###### Lemma 3.2

Under the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), the bilinear form is continuous on .

Proof. By using the Cauchy-Schwarz inequality, trace theorem and the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), we have

 A(u––,v–) ≤2ν∥→us∥→H1(Ds))∥→vs∥→H1(Ds))+gdKmax∥ϕm∥H1(Dm)∥ψm∥H1(Dm) +g∥ϕm∥H1(Dm)∥→vs∥→H1(Ds)+g∥ψm∥H1(Dm)∥→us∥→H1(Ds) +α√gν√Kmin∥→us∥→H1(Ds)∥→vs∥L2(→H1(Ds)+αdKmax√gν√Kmin∥ϕm∥H1(Dm)∥→vs∥→H1(Ds),

for . Thus the bilinear form is continuous on the space .

###### Lemma 3.3

The linear form is continuous on .

Proof. By using the Cauchy-Schwarz inequality and trace theorem, we have

 F(v–)≤∥→fs∥→H1(Ds)∥→vs∥→H1(Ds)+g∥fm∥L2(Dm)∥ψm∥H1(Dm)+gz∥→vs∥→H1(Ds),

for . Thus the linear form is continuous on .

###### Lemma 3.4

Under the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), the bilinear form is coercive on when the coefficient in the Beavers-Joseph condition (LABEL:stochastic_stokes_Darcy_7) is small enough.

Proof. By using the Korn’s inequality, Poincar inequality, Cauchy-Schwarz inequality, trace theorem and the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), we have

 A(u––,u––)= ∫Ω∫Ds2νD(→us):D(→us)dDsdΩ+g∫Ω∫Dm(K∇ϕm)⋅(∇ϕm)dDmdΩ + ∫Ω∫ΓIαν√d√trace(Π)(Pτ(→us)⋅→us+Pτ(K∇ϕm)⋅→us)dΓIdΩ ≥2C1ν∥→us∥2→H1(Ds)+C2gKmin∥ϕm∥2H1(Dm)−αdKmax√gν√Kmin∥ϕm∥H1(Dm)∥→us∥→H1(Ds) ≥C1ν∥→us∥2→H1(Ds)+12C2gKmin∥ϕm∥2H1(Dm),

where , for . Thus the bilinear form is coercive on when the coefficient in the Beavers-Joseph (LABEL:stochastic_stokes_Darcy_7) condition is small enough.

###### Theorem 3.5

Under the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), there exists a unique weak solution and up to an additive constant for the weak formulation (LABEL:weak_formulation) of stochastic Stoke-Darcy interface problem (LABEL:Darcy_law_4_stochastic)-(LABEL:stokes_2_stochastic) when the coefficient in the Beavers-Joseph (LABEL:stochastic_stokes_Darcy_7) condition is small enough.

Proof. Based on the Lemma LABEL:lemma_A_1, Lemma LABEL:lemma_F and Lemma LABEL:lemma_A_2, there exists a unique weak solution by the Lax-Milgram Lemma. Then the assertion about is clear, which drives form the conclusions in the deterministic scenario [45, 62, 72].

## 4 Numerical solution for the stochastic coupled problem

Since the moments are the characteristic functions of the stochastic solution, the object is to design a numerical method to calculate the moments of the stochastic solution. The main difficulty in this design is how to represent the stochastic solution by a discrete form in the probability space and the physical space. For the discrete form in the probability space, we choose the ensemble representations in sampling methods, e.g., Monte Carlo (MC) method in this paper. But the total computational cost of the traditional single-level Monte Carlo (SLMC) method is very high. Then the multi-level Monte Carlo (MLMC) method is adopted to reduce the total computational cost in the probability space. For the discrete form in the physical space, the finite element method (FEM) is chosen. Furthermore the multi-grid (MG) method is used to reduce the computational cost in physical space. Thus the multi-grid multi-level Monte Carlo (MGMLMC) method is developed to reduce the computational cost both in the probability space and the physical space.

### 4.1 Realizations of the random hydraulic conductivity

The realizations of the random hydraulic conductivity in a discrete form on the spatial domain and the random fields are the basises of the numerical method. We adopt the grid based method in , because this method represents the random field exactly at the discrete points without any truncation.

For simplification, we assume is a diagonal matrix. The process to generate the realizations of is displayed as follows, which is as same as the processes to generate the realizations of .

Because is physical positive, we assume is a log-normal distribution, i.e.,

 K(ω,x)=eZ(ω,x),  ω∈Ω, x∈¯Dm, \hb@xt@.01(4.1)

where is a mean zero Gaussian random field on , with the continuous covariance function , , i.e.,

 E[Z(ω,x)] = 0,∀x∈¯Dm, \hb@xt@.01(4.2) E[Z(ω,x),Z(ω,y)] = r(x,y),∀x,y∈¯Dm. \hb@xt@.01(4.3)

For , , the vector represents all the discrete spatial points in , on which is provided as , . By the covariance function (LABEL:variance_Z), a positive definite matrix is given

 R=E[Z(ω,→x),Z(ω,→x)⊤]=(r(xi,xj))Mi,j=1. \hb@xt@.01(4.4)

Let be the Cholesky factorization of as . Then we can generate the realizations of at the discrete points without any truncation by

 Z(ω,→x)=ΘY, \hb@xt@.01(4.5)

where is a vector of independent identically distributed standard Gaussian random variables. It is easy to verify that , and . And the realizations of are generated by the formulation (LABEL:log_normal_K). Some samples of the realizations of the random hydraulic conductivity will be displayed in the latter section.

### 4.2 Monte Carlo methods

The Monte Carlo method  is a classical method to calculate the numerical approximation of moments. In this paper, we only investigate the process to generate the expected value of , and , which is easy to be used for the high order of moments.

For simplification, the symbol is used to substitute the quantity of interest (QoI) of , and . Let denote the finite element approximation of on the quasi-uniform triangulation mesh with the mesh size , and denote the realization of with the sample . Then the approximation of the expected value of by SLMC method with samples is given as:

 ^QSLℓ(x)=1NSLℓNSLℓ∑i=1Qiℓ(x). \hb@xt@.01(4.6)

When no ambiguity arises, we may omit in , and for convenience.

The mean squared error of the SLMC method is:

 MSE(^QSLL)=E[(^QSLL−E[Q])2]=E[(^QSLL−E[QL]+E[QL]−E[Q])2]≤2E[(^QSLL−E[QL])2]+2E[(E[QL]−E[Q])2]=2V[QL]NSLL+2(E[QL]−E[Q])2. \hb@xt@.01(4.7)

Then the error of SLMC method with a given norm is bounded as

 ∥MSE(^QSLL)∥≤2∥V[QL]∥NSLL+2∥(E[QL]−E[Q])2∥, \hb@xt@.01(4.8)

i.e., the accuracy of SLMC method is based on the sampling error and the FEM error.

### 4.3 Multi-level Monte Carlo methods

The total computational cost of single-level Monte Carlo is

 TSLc=NSLLCL, \hb@xt@.01(4.9)

where is the computational cost of one sample with mesh size . would be very high when and are both very large. By the accuracy formulation (LABEL:normal_Error_Monte_Carlo_FE_for_QoI) of SLMC method, the sampling error and the FEM error should be both small enough, if a small mean squared error is required. Thus should be larger while the mesh size becomes smaller. On the other hand, increase exponentially as the mesh size becomes smaller. Thus the total computational cost increases very fast as mesh size become smaller. An efficient algorithm is needed to reduce the total computational cost. We adopt the multi-level Monte Carlo (MLMC) method.

By the linearity of the expectation operator

 E[QL]=E[Q0]+L∑ℓ=1E[Qℓ]−[Qℓ−1]=E[Q0]+L∑ℓ=1E[Qℓ−Qℓ−1]. \hb@xt@.01(4.10)

Then we can use the hierarchical meshes to construct the MLMC method to generate the expect value of . Let be a sequence of quasi-uniform triangulation meshes with the mesh sizes . These mesh sizes satisfy , . And are the numbers of samples with the mesh sizes . Then the approximation of the expected value by the MLMC method is given by:

 ^QMLL=1NML0NML0∑i=1Qi0+L∑ℓ=11NMLℓNMLℓ∑i=1(Qiℓ−Qiℓ−1), \hb@xt@.01(4.11)

and the corresponding mean squared error of the MLCM method with norm is

 ∥MSE(^QMLL)∥=∥E[(^QMLL−E[Q])2]∥≤2∥V[Q0]∥NML0+2L∑ℓ=1∥V[Qℓ−Qℓ−1]∥NMLℓ+2∥(E[QL]−E[Q])2∥. \hb@xt@.01(4.12)

For simplicity, let , , , , and be the computational cost of generating one sample of , . Then the mean squared error is rewrote as

 ∥MSE(^QMLL)∥≤2L∑ℓ=0vℓNMLℓ+2∥(E[QL]−E[Q])2∥. \hb@xt@.01(4.13)

And the total computational cost is

 TMLc=L∑ℓ=0NMLℓCℓ . \hb@xt@.01(4.14)

By the mean squared error of SLMC method (LABEL:normal_Error_Monte_Carlo_FE_for_QoI) and MLMC method (LABEL:MSE_MLMC), the accuracy of approximation of expected value is based on two parts, i.e., the sampling error and FEM error. The FEM error is fixed when the mesh size is given. Thus the sampling error should be small enough with the given mesh size . We substitute the sampling errors in SLMC method and MLMC method by:

 eSLL=∥V[QL]∥NL,\ \ and \ \ eMLL=L∑ℓ=0vℓNMLℓ . \hb@xt@.01(4.15)

For guaranteeing the accuracy of MLMC method is as same as SLMC method, the following relationship between two sampling errors should be ensured

 eMLL≤eSLL,\ i.e., \ v0NML0+v1NML1+