Multigrid MultiLevel Monte Carlo Method for StokesDarcy interface Model with Random Hydraulic Conductivity^{†}^{†}thanks: This work is partially supported by NSF grants DMS1418624 and DMS1722647, National Science Foundation of China (91330104).
Abstract
In this article we develop a multigrid multilevel Monte Carlo (MGMLMC) method for the stochastic StokesDarcy 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 StokesDarcy 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 singlelevel 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 multilevel 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 StokesDarcy model. Furthermore MLMC naturally provides the hierarchial grids and sufficient information on these grids for multigrid (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 multigrid multilevel 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 StokesDarcy interface model, multilevel Monte Carlo, multigrid method.
AMS subject classifications. 35R60, 65C05, 65M60 76S05.
1 Introduction
The StokesDarcy 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 StokesDarcy 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], multigrid 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 StokesDarcy 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 reallife 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 StokesDarcy 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 [75] 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 multilevel Monte Carlo method [7, 20, 43, 44, 60, 83] to solve the sophisticated stochastic StokesDarcy 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 multilevel Monte Carlo method only reduces the computational cost in the probability space, not in the physical space. Inspired by a fact that the multilevel Monte Carlo method already has a set of hierarchical grids for the multilevel 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 multigrid method [11, 65, 70, 87, 88], which can further improve the efficiency of the proposed multilevel Monte Carlo method. Meanwhile, the saved information of the multilevel Monte Carlo method on the set of hierarchical grids will also significantly reduce the computational cost of the multigrid method. Therefore, we combine the multilevel Monte Carlo method and the multigrid method on the same set of hierarchical grids to propose an even more costly efficient method, which is the multigrid multilevel Monte Carlo method.
The rest of the paper is organized as follows. In section 2, we briefly recall the deterministic StokesDarcy model. In section 3, we present the stochastic StokesDarcy interface model, the weak formulation of the stochastic StokesDarcy model and the proof of the wellposedness. In section 4, we recall the Monte Carlo method to approximate the numerical moments of the stochastic solutions, adopt the multilevel Monte Carlo method to reduce the computational cost in probability space, and then develop the multigrid multilevel 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 StokesDarcy 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 StokesDarcy 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 [6]
\hb@xt@.01(2.1)  
\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 secondorder form of the Darcy system
\hb@xt@.01(2.3) 
In the conduit domain , the flow is governed by the Stokes equations:
\hb@xt@.01(2.4)  
\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:
\hb@xt@.01(2.6)  
\hb@xt@.01(2.7)  
\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 BeaversJoseph condition [8, 14, 16, 17, 36, 53].
3 StokesDarcy 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 wellposedness 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 [35] 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 seminorm .
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
Here is the norm given as, ,
which is induced by following inner product, ,
For , we denote the Hilbert space and with the standard norm and seminorm . For , we denote with the standard norm . For , we denote and . For the vector , , 2norm of is
For simplicity, we define
3.2 Stochastic StokesDarcy interface equations
With the complete probability space , let , , be a random hydraulic conductivity tensor.
Then in the porous media domain, the stochastic secondorder form of Darcy equation with sink/source term is given as:
\hb@xt@.01(3.1) 
And the interface conditions are modified as:
\hb@xt@.01(3.2)  
\hb@xt@.01(3.3)  
\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
\hb@xt@.01(3.5)  
\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 velocitypressure spaces on the conduit domain as
and we denote the pressure space on the porous media as
For convenience, let , , , and , where , . The norms of are given as
\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
\hb@xt@.01(3.8) 
where
\hb@xt@.01(3.9)  
\hb@xt@.01(3.10)  
\hb@xt@.01(3.11)  
\hb@xt@.01(3.12)  
\hb@xt@.01(3.13)  
\hb@xt@.01(3.14)  
\hb@xt@.01(3.15)  
\hb@xt@.01(3.16) 
3.4 Wellposedness of the weak solution
The approach to analyze the wellposedness 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
\hb@xt@.01(3.17) 
the integrability condition: let and satisfy
\hb@xt@.01(3.18)
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 CauchySchwarz inequality, trace theorem and the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), we have
for . Thus the bilinear form is continuous on the space .
Lemma 3.3
The linear form is continuous on .
Proof. By using the CauchySchwarz inequality and trace theorem, we have
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 BeaversJoseph condition (LABEL:stochastic_stokes_Darcy_7) is small enough.
Proof. By using the Korn’s inequality, Poincar inequality, CauchySchwarz inequality, trace theorem and the Assumption (LABEL:K_assumption_1) or (LABEL:K_assumption_2), we have
where , for . Thus the bilinear form is coercive on when the coefficient in the BeaversJoseph (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 StokeDarcy interface problem (LABEL:Darcy_law_4_stochastic)(LABEL:stokes_2_stochastic) when the coefficient in the BeaversJoseph (LABEL:stochastic_stokes_Darcy_7) condition is small enough.
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 singlelevel Monte Carlo (SLMC) method is very high. Then the multilevel 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 multigrid (MG) method is used to reduce the computational cost in physical space. Thus the multigrid multilevel 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 [48], 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 lognormal distribution, i.e.,
\hb@xt@.01(4.1) 
where is a mean zero Gaussian random field on , with the continuous covariance function , , i.e.,
\hb@xt@.01(4.2)  
\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
\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
\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 [75] 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 quasiuniform 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:
\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:
\hb@xt@.01(4.7) 
Then the error of SLMC method with a given norm is bounded as
\hb@xt@.01(4.8) 
i.e., the accuracy of SLMC method is based on the sampling error and the FEM error.
4.3 Multilevel Monte Carlo methods
The total computational cost of singlelevel Monte Carlo is
\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 multilevel Monte Carlo (MLMC) method.
By the linearity of the expectation operator
\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 quasiuniform 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:
\hb@xt@.01(4.11) 
and the corresponding mean squared error of the MLCM method with norm is
\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
\hb@xt@.01(4.13) 
And the total computational cost is
\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:
\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