Multi-grid Multi-Level Monte Carlo Method for Stokes-Darcy interface Model with Random Hydraulic Conductivity

TL;DR

This paper introduces a novel multi-grid multi-level Monte Carlo method to efficiently solve the stochastic Stokes-Darcy interface model with random hydraulic conductivity, reducing computational costs while maintaining accuracy.

Contribution

The paper develops a combined MGMLMC approach for the stochastic Stokes-Darcy model, integrating multi-grid and multi-level Monte Carlo techniques for improved efficiency.

Findings

The MGMLMC method significantly reduces computational cost compared to traditional methods.

Numerical examples confirm the accuracy and efficiency of the proposed approach.

Theoretical analysis guarantees the existence and uniqueness of the solution.

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 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. Numerical examples are provided to verify and illustrate the proposed method and the theoretical conclusions.

Tables

Table 1. Table 1: Mean values of 𝜷 𝜷 \beta with different 𝝈 𝝈 \sigma and norm

$𝝈$	0.02	0.8	1.2
$\parallel \cdot {\parallel }_{{𝑳}_{\mathrm{𝟐}}}$	2.0204	2.0216	2.0209
$\parallel \cdot {\parallel }_{{𝑳}_{\infty }}$	1.6468	1.6487	1.6511
$\parallel \cdot {\parallel }_{{𝑯}_{\mathrm{𝟏}}}$	1.3030	1.3043	1.3081

Table 2. Table 2: Computational cost with different mesh size

h	$\mathrm{𝟏}/\mathrm{𝟒}$	$\mathrm{𝟏}/\mathrm{𝟖}$	$\mathrm{𝟏}/\mathrm{𝟏𝟔}$	$\mathrm{𝟏}/\mathrm{𝟑𝟐}$
cpu time $\left(\text{sec.}\right)$	0.48	1.84	9.64	245.65
tic-toc time $\left(\text{sec.}\right)$	0.44	1.43	8.38	240.46

Table 3. Table 3: Number of samples at every level

	${𝒆}_{\mathrm{𝟑}}^{𝑺𝑳}$	${𝒆}_{\mathrm{𝟑}}^{𝑴𝑳}$	${𝑵}_{\mathrm{𝟎}}^{𝑴𝑳}$	${𝑵}_{\mathrm{𝟏}}^{𝑴𝑳}$	${𝑵}_{\mathrm{𝟐}}^{𝑴𝑳}$	${𝑵}_{\mathrm{𝟑}}^{𝑴𝑳}$	${𝑵}_{\mathrm{𝟑}}^{𝑺𝑳}$
$\parallel \cdot {\parallel }_{{𝑳}^{\mathrm{𝟐}}}$	$\mathrm{𝟐}\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟕}}$	$1.95\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟕}}$	2127	504	83	14	122
$\parallel \cdot {\parallel }_{{𝑳}^{\infty }}$	$\mathrm{𝟐}\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	$1.98\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	2602	701	131	24	139
$\parallel \cdot {\parallel }_{{𝑯}^{\mathrm{𝟏}}}$	$\mathrm{𝟑}\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	$2.97\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$	3521	1071	225	47	146

Table 4. Table 4: Relative errors of solutions by SLMC and MGMLMC method

	${\varphi }_{𝒎}$	${𝒑}_{𝒔}$	${𝒖}_{𝒎}^{\mathrm{𝟏}}$	${𝒖}_{𝒎}^{\mathrm{𝟐}}$	${𝒖}_{𝒔}^{\mathrm{𝟏}}$	${𝒖}_{𝒔}^{\mathrm{𝟐}}$
$𝜷=2.02$	$3.39\%$	$0.02\%$	$0.02\%$	$0.02\%$	$4.31\%$	$3.24\%$
$𝜷=1.65$	$3.08\%$	$0.02\%$	$0.02\%$	$0.02\%$	$3.75\%$	$2.78\%$
$𝜷=1.30$	$2.5\%$	$0.01\%$	$0.016\%$	$0.015\%$	$2.65\%$	$2.07\%$

Table 5. Table 5: Efficiency of MLMC and MGMLMC

	${𝑻}_{𝒄}^{𝑺𝑳}$ (sec.)	${𝑻}_{𝒄}^{𝑴𝑳}$ (sec.)	${𝑻}_{𝒄}^{𝑴𝑮𝑴𝑳}$ (sec.)	${𝑻}_{𝒄}^{𝑴𝑳}/{𝑻}_{𝒄}^{𝑺𝑳}$	${𝑻}_{𝒄}^{𝑴𝑮𝑴𝑳}/{𝑻}_{𝒄}^{𝑺𝑳}$
$\parallel \cdot {\parallel }_{{𝑳}^{\infty }}$	30791	7315	1896	$23.76\%$	$6.16\%$
$\parallel \cdot {\parallel }_{{𝑳}^{\mathrm{𝟐}}}$	27025	4436	1226	$16.41\%$	$4.54\%$
$\parallel \cdot {\parallel }_{{𝑯}^{\mathrm{𝟏}}}$	32342	13701	3324	$42.36\%$	$10.28\%$

Equations

$${\vec{u}}_m(x)$$

$$\nabla \cdot {\vec{u}}_m(x)$$

$$-\nabla \cdot (\mathbb{K}(x)\nabla {\varphi }_m(x))=f_m(x)\text{in}{D}_m.$$

$$-\nabla \cdot \mathbb{T}({\vec{u}}_s,p_s)$$

$$\nabla \cdot {\vec{u}}_s$$

$${\vec{u}}_s\cdot {\vec{n}}_s$$

$$-{\vec{n}}_s^T\mathbb{T}({\vec{u}}_s,p_s){\vec{n}}_s$$

$$-{\mathbf{\tau }}_j^T\mathbb{T}({\vec{u}}_s,p_s){\vec{n}}_s$$

$$L^2\left(\Omega ;W^{r,q}(D)\right):=\{\varphi :\Omega \to W^{r,q}(D)|\varphi \text{strongly measurable},\Vert \varphi {\Vert }_{L^2\left(\Omega ;W^{r,q}(D)\right)}<\infty \}.$$

$$\Vert \varphi {\Vert }_{L^2\left(\Omega ;W^{r,q}(D)\right)}:={\left(\mathbb{E}\left[\Vert \varphi (\omega ,\cdot ){\Vert }_{W^{r,q}(D)}^2\right]\right)}^{1/2}:={\left({\int }_{\Omega }\Vert \varphi (\omega ,\cdot ){\Vert }_{W^{r,q}(D)}^{2}d\mathcal{P}(\omega )\right)}^{1/2},$$

$$[\varphi ,\psi {]}_{L^2(\Omega ;W^{r,q}(D))}:=\mathbb{E}\left[(\varphi ,\psi {)}_{W^{r,q}(D)}\right]:={\int }_{\Omega }(\varphi ,\psi {)}_{W^{r,q}(D)}d\mathcal{P}(\omega ).$$

$${\mathcal{L}}^q(D)$$

$${\mathcal{H}}^r(D)$$

$${\vec{\mathcal{H}}}^r(D)$$

$$-\nabla\cdot\big{(}\mathbb{K}(\omega,x)\nabla\phi_{m}(\omega,x)\big{)}=f_{m}(x),\hskip 11.38092pt\text{in }D_{m}.$$

$${\vec{u}}_s(\omega ,x)\cdot {\vec{n}}_s(x)$$

$$-{\vec{n}}_s^{\top }\mathbb{T}({\vec{u}}_s,p_s){\vec{n}}_s$$

$$-{\mathbf{\tau }}_j^{\top }\mathbb{T}({\vec{u}}_s,p_s){\vec{n}}_s$$

$$-\nabla \cdot \mathbb{T}({\vec{u}}_s(\omega ,x),p_s(\omega ,x))$$

$$\nabla \cdot {\vec{u}}_s(\omega ,x)$$

$$X_s$$

$$X_s^0$$

$$X_{s,div}$$

$$Q_s$$

$$\begin{gathered} X_m=\{{\varphi }_m\in {\mathcal{H}}^1(D_m)|{\varphi }_m=0\text{on}{\Gamma }_m\}, \\ {X}_m^0=\{{\varphi }_m\in {\mathcal{H}}^0(D_m)|{\varphi }_m=0\text{on}{\Gamma }_m\}. \end{gathered}$$

$$\Vert \underline{u}{\Vert }_{X^r}={\left(\mathbb{E}\left[\Vert \underline{u}{\Vert }_{{\mathbf{H}}^r(D_s)\times H^r(D_m)}^2\right]\right)}^{1/2}={\left(\Vert {\vec{u}}_s{\Vert }_{{\vec{\mathcal{H}}}^r(D_s)}^2+\Vert {\varphi }_m{\Vert }_{{\mathcal{H}}^r(D_m)}^2\right)}^{1/2},r=0,1.$$

$$\{\begin{array}{rl}A(\underline{u},\underline{v})-B(\underline{v},p_s)& =F(\underline{v}),\forall \underline{v}=({\vec{v}}_s,{\psi }_m)\in X,\\ B(\underline{u},q_s)& =0,\forall q_s\in Q_s,\end{array}$$

$$A(\underline{u},\underline{v})$$

$$a(\underline{u},\underline{v})$$

$$+$$

$$-$$

$$B(\underline{v},p_s)$$

$$b(\underline{v},p_s)$$

$$F(\underline{v})$$

$$f(\underline{v})$$

$$0<K_{min}\leq\big{\{}K_{ii}(\omega,x)\big{\}}^{d}_{i=1}\leq K_{max}<\infty,\text{\ \ for\ \ }(\omega,x)\in\Omega\times\bar{D}_{m};$$

$$0<K_{min}(\omega )\text{ and }\frac{1}{K_{min}(\omega )},K_{max}(\omega )\in L^{\infty }(\Omega )\text{ for a.e. }\omega \in \Omega .$$

$$A(\underline{u},\underline{v})$$

$$+g\Vert {\varphi }_m{\Vert }_{{\mathcal{H}}^1(D_m)}\Vert {\vec{v}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)}+g\Vert {\psi }_m{\Vert }_{{\mathcal{H}}^1(D_m)}\Vert {\vec{u}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)}$$

$$+\frac{\alpha \sqrt{g\nu }}{\sqrt{K_{min}}}\Vert {\vec{u}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)}\Vert {\vec{v}}_s{\Vert }_{L^2({\vec{\mathcal{H}}}^1(D_s)}+\frac{\alpha d{K}_{max}\sqrt{g\nu }}{\sqrt{K_{min}}}\Vert {\varphi }_m{\Vert }_{{\mathcal{H}}^1(D_m)}\Vert {\vec{v}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)},$$

$$F(\underline{v})\le \Vert {\vec{f}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)}\Vert {\vec{v}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)}+g\Vert f_m{\Vert }_{{\mathcal{L}}^2(D_m)}\Vert {\psi }_m{\Vert }_{{\mathcal{H}}^1(D_m)}+gz\Vert {\vec{v}}_s{\Vert }_{{\vec{\mathcal{H}}}^1(D_s)},$$

Full text

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 [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 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 $D_{ms}=D_{m}\cup D_{s}\subset{\mathbb{R}}^{d}$, $d=2,3$, where $D_{m}$ is the porous media domain and $D_{s}$ is the conduit domain. We decompose the boundary $\partial D$ into two parts: $\Gamma_{m}=\partial D_{m}\backslash\Gamma_{I}$, $\Gamma_{s}=\partial D_{s}\backslash\Gamma_{I}$, and denote the interface as $\Gamma_{I}=\partial D_{m}\cap\partial D_{s}$.

In the porous media domain $D_{m}$, the flow is governed by the Darcy system [6]

[TABLE]

here, $\vec{u}_{m}$ denotes the specific discharge in the porous media, $\mathbb{K}$ is the hydraulic conductivity tensor of the porous media that is symmetric and positive definite in accordance with physical meaning, $\phi_{m}$ is the hydraulic head, and $f_{m}$ is the sink/source term.

By substituting (1) into (2), we obtain the second-order form of the Darcy system

[TABLE]

In the conduit domain $D_{s}$, the flow is governed by the Stokes equations:

[TABLE]

where $\vec{u}_{s}$ denotes the fluid velocity, $p_{s}$ is the kinematic pressure, and $\vec{f_{s}}$ is the external body force. $\mathbb{T}$ is the stress tensor, defined as $\mathbb{T}(\vec{u}_{s},p_{s})=2\nu\mathbb{D}(\vec{u}_{s})-p_{s}\mathbb{I}$, where $\nu$ is the kinematic viscosity of the fluid and $\mathbb{D}(\vec{u}_{s})=\frac{1}{2}(\nabla\vec{u}_{s}+\left({\nabla}\vec{u}_{s}\right)^{T})$.

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

[TABLE]

where $\vec{n}_{s}$, $\vec{n}_{m}$ denote the unit outer normal to the conduit and the porous media regions at the interface $\Gamma_{I}$, respectively, \tau$${}_{j}(j=1,...,d-1) denote mutually orthogonal unit tangential vectors to the interface $\Gamma_{I}$, $z$ is the hight, $g$ is the gravitational acceleration, and $\Pi(x)=\frac{\mathbb{K}(x)\nu}{g}$ is the intrinsic permeability. The first interface condition (6) is governed by the conservation of mass, the second interface condition (7) 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 (8) 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 $\Pi$. Thus the hydraulic conductivity tensor $\mathbb{K}(x)$ is also a random field with the relationship $\Pi=\frac{\mathbb{K}\nu}{g}$. 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 [35] for the classical Sobolev spaces. Let $D$ be an open, connected, bounded, and convex subset of $\mathbb{R}^{d}$, $d=2,3$, with polygonal and Lipschitz continuous boundary $\partial D$. Let $r\in\mathbb{R}$, $q\in\mathbb{Z}$, and $W^{r,q}(D)$ be a Sobolev space on $D$ with the standard norm $\|\cdot\|_{W^{r,q}(D)}$ and semi-norm $\mid\cdot\mid_{W^{r,q}(D)}$.

Let $(\Omega,\mathcal{F},\mathcal{P})$ be a complete probability space. Here $\Omega$ is the set of outcomes, $\mathcal{F}$ is the $\sigma$-algebra of events, and $\mathcal{P}:\mathcal{F}\rightarrow[0,1]$ is a probability measure.

For the given probability space $(\Omega,\mathcal{F},\mathcal{P})$ and the Sobolev space $W^{r,q}(D)$ with the inner product $(\cdot,\cdot)_{W^{r,q}(D)}$ and norm $\|\cdot\|_{W^{r,q}(D)}$, we define the stochastic Sobolev space, which consists of strongly measurable, $r$-summable mappings $\phi:\Omega\rightarrow W^{r,q}(D)$, by

[TABLE]

Here $\|\cdot\|_{L^{2}\left(\Omega;H^{r}(D)\right)}$ is the norm given as, $\forall\phi\in L^{2}\left(\Omega;W^{r,q}(D)\right)$,

[TABLE]

which is induced by following inner product, $\forall\phi,\psi\in L^{2}(\Omega;W^{r,q}(D))$,

[TABLE]

For $q=2$, we denote the Hilbert space $H^{r}(D):=W^{r,2}(D)$ and $H^{r}_{0}(D):=\{u:u\in H^{r}(D),u\mid_{\partial D}=0\}$ with the standard norm $\|\cdot\|_{H^{r}(D)}$ and semi-norm $\mid\cdot\mid_{H^{r}(D)}$. For $r=2$, we denote $L^{q}(D):=W^{2,q}(D)$ with the standard norm $\|\cdot\|_{L^{2}(D)}$. For $d=2,3$, we denote $\mathbf{H}^{r}(D):=\left[H^{r}(D)\right]^{d}$ and $\mathbf{L}^{q}(D):=\left[L^{q}(D)\right]^{d}$. For the vector $\vec{v}=(v_{1},v_{2},\cdots,v_{n})^{\top}$, $n\in\mathbb{N}^{+}$, 2-norm $\|\vec{v}\|_{2}$ of $\vec{v}$ is $\|\vec{v}\|_{2}=\left(v_{1}^{2}+v_{2}^{2}+\cdots+v_{n}^{2}\right)^{1/2}.$

For simplicity, we define

[TABLE]

3.2 Stochastic Stokes-Darcy interface equations

With the complete probability space $(\Omega,\mathcal{F},\mathcal{P})$, let $\mathbb{K}(\omega,x)$, $\omega\in\Omega$, $x\in D_{m}$ be a random hydraulic conductivity tensor.

Then in the porous media domain, the stochastic second-order form of Darcy equation with sink/source term $f_{m}(x)$ is given as:

[TABLE]

And the interface conditions are modified as:

[TABLE]

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

[TABLE]

For the boundary conditions, we assume the hydraulic head $\phi_{m}$ and the fluid velocity $\vec{u}_{s}$ satisfy homogeneous Dirichlet boundary condition except on $\Gamma_{I}$.

3.3 Weak formulation of the coupled problem

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

[TABLE]

and we denote the pressure space on the porous media as

[TABLE]

For convenience, let $X^{1}=X=X_{s}\times X_{m}$, $X_{div}=X_{s,div}\times X_{m}$, $X^{0}=X^{0}_{s}\times X^{0}_{m}$, and $\underline{u}=(\vec{u}_{s},\phi_{m})\in X$, where $\vec{u}_{s}\in X_{s}$, $\phi_{m}\in X_{m}$. The norms of $X^{r},r=0,1$ are given as

[TABLE]

The projection onto the local tangential plane of the vector $\vec{u}$ is denoted as $P_{\tau}(\vec{u})=\vec{u}-(\vec{u}\cdot\vec{n}_{s})\vec{n}_{s}$. Then using the boundary conditions (10)-(12), we obtain the following weak formulation: find $(\underline{u},p_{s})\in X\times Q_{s}$, such that

[TABLE]

where

[TABLE]

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 1**.**

Let $\mathbb{K}(\omega,x)$ be a diagonal matrix as diag\big{(}K_{11}(\omega,x),\cdots,K_{dd}(\omega,x)\big{)},\omega\in\Omega,x\in D_{m},d=2,3.

• •

the strong elliptic condition: there are positive lower and upper bounds $K_{min}$, $K_{max}$ such that

[TABLE]

• •

the integrability condition: let K_{min}(\omega):=\min\limits_{x\in\bar{D}_{m}}\big{\{}K_{ii}(\omega,x)\big{\}}^{d}_{i=1} and K_{max}(\omega):=\max\limits_{x\in\bar{D}_{m}}\big{\{}K_{ii}\\ (\omega,x)\big{\}}^{d}_{i=1} satisfy

[TABLE]

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

Lemma 2**.**

Under the Assumption (25) or (26), the bilinear form $A(\cdot,\cdot)$ is continuous on $X_{div}\times X_{div}$.

Proof.

By using the Cauchy-Schwarz inequality, trace theorem and the Assumption (25) or (26), we have

[TABLE]

for $\forall\underline{u},\underline{v}\in X_{div}$. Thus the bilinear form $A(\cdot,\cdot)$ is continuous on the space $X_{div}\times X_{div}$. ∎

Lemma 3**.**

The linear form $F(\cdot)$ is continuous on $X_{div}$.

Proof.

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

[TABLE]

for $\forall\underline{v}\in X_{div}$. Thus the linear form $F(\cdot)$ is continuous on $X_{div}$. ∎

Lemma 4**.**

Under the Assumption (25) or (26), the bilinear form $A(\cdot,\cdot)$ is coercive on $X_{div}\times X_{div}$ when the coefficient $\alpha$ in the Beavers-Joseph condition (12) is small enough.

Proof.

By using the Korn’s inequality, Poincar$\acute{\text{e}}$ inequality, Cauchy-Schwarz inequality, trace theorem and the Assumption (25) or (26), we have

[TABLE]

where $\alpha^{2}\leq\frac{2C_{1}C_{2}K^{2}_{min}}{d^{2}K^{2}_{max}}$, for $\forall\underline{u}\in X_{div}$. Thus the bilinear form $A(\cdot,\cdot)$ is coercive on $X_{div}\times X_{div}$ when the coefficient $\alpha$ in the Beavers-Joseph (12) condition is small enough. ∎

Theorem 5**.**

Under the Assumption (25) or (26), there exists a unique weak solution $\underline{u}=(\vec{u}_{s},\phi_{m})\in X$ and $p_{s}$ up to an additive constant for the weak formulation (16) of stochastic Stoke-Darcy interface problem (9)-(14) when the coefficient $\alpha$ in the Beavers-Joseph (12) condition is small enough.

Proof.

Based on the Lemma 2, Lemma 3 and Lemma 4, there exists a unique weak solution $\underline{u}$ by the Lax-Milgram Lemma. Then the assertion about $p_{s}$ 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 $\mathbb{K}(\omega,x)$ in a discrete form on the spatial domain $D_{m}$ and the random fields $\Omega$ 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 $x_{1},\cdots,x_{M}\in D_{m}$ without any truncation.

For simplification, we assume \mathbb{K}(\omega,x)=diag\big{(}K_{11}(\omega,x),\cdots,K_{dd}(\omega,x)\big{)},\omega\in\Omega,x\in D_{m},d=2,3 is a diagonal matrix. The process to generate the realizations of $K(\omega,x)=K_{11}(\omega,x)$ is displayed as follows, which is as same as the processes to generate the realizations of $K_{ii}(\omega,x),i=2,3$.

Because $K(\omega,x)$ is physical positive, we assume $K(\omega,x)$ is a log-normal distribution, i.e.,

[TABLE]

where $Z(\omega,x)$ is a mean zero Gaussian random field on $\bar{D}_{m}$, with the continuous covariance function $r(x,y)$, $x,y\in\bar{D}_{m}$, i.e.,

[TABLE]

For $x_{i}\in\bar{D}_{m}$, $i=1,2,\cdots,M$, the vector $\vec{x}=(x_{1},x_{2},\cdots,x_{M})^{\top}$ represents all the discrete spatial points in $\bar{D}_{m}$, on which $Z(\omega,x)$ is provided as $Z(\omega,\vec{x})=(Z_{1},Z_{2},\cdots,Z_{M})^{\top}$, $Z_{i}=Z(\omega,x_{i})$. By the covariance function (29), a $M\times M$ positive definite matrix $R$ is given

[TABLE]

Let $\Theta$ be the Cholesky factorization of $R$ as $R=\Theta{\Theta}^{\top}$. Then we can generate the realizations of $Z(\omega,\vec{x})$ at the discrete points $\vec{x}$ without any truncation by

[TABLE]

where $Y:=(Y_{1}(w),...,Y_{M}(w))^{\top}$ is a $M\times 1$ vector of independent identically distributed standard Gaussian random variables. It is easy to verify that $\mathbb{E}[Z(\omega,x)]=\mathbb{E}[\Theta Y]=\Theta\mathbb{E}[Y]=\mathbf{0}$, and $\mathbb{E}[Z(\omega,\vec{x}),Z(\omega,\vec{x})^{\top}]=\mathbb{E}[(\Theta Y)(\Theta Y)^{\top}]=\Theta\mathbb{E}[YY^{\top}]\Theta^{\top}=\Theta\Theta^{\top}=R$. And the realizations of $K(\omega,x)$ are generated by the formulation (27). Some samples of the realizations of the random hydraulic conductivity $K$ 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 $\phi_{m}$, $\vec{u}_{s}$ and $p_{s}$, which is easy to be used for the high order of moments.

For simplification, the symbol $Q$ is used to substitute the quantity of interest (QoI) of $\phi_{m}$, $\vec{u}_{s}$ and $p_{s}$. Let $Q_{\ell}(\omega,x)$ denote the finite element approximation of $Q(\omega,x)$ on the quasi-uniform triangulation mesh $\mathcal{T}_{\ell}$ with the mesh size $h_{\ell}$, and $Q^{i}_{\ell}(x)$ denote the realization of $Q_{\ell}(\omega,x)$ with the sample $\mathbb{K}(\omega^{i},x)$. Then the approximation $\hat{Q}^{SL}_{\ell}(x)$ of the expected value of $Q$ by SLMC method with $N_{\ell}^{SL}$ samples $\{\mathbb{K}(\omega^{i},x)\}^{N^{\ell}_{L}}_{i=1}$ is given as:

[TABLE]

When no ambiguity arises, we may omit $x$ in $Q_{\ell}(x)$, $Q_{\ell}^{i}(x)$ and $\hat{Q}_{\ell}(x)$ for convenience.

The mean squared error of the SLMC method is:

[TABLE]

Then the error of SLMC method with a given norm $\|\cdot\|$ is bounded as

[TABLE]

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 $T^{SL}_{c}$ of single-level Monte Carlo is

[TABLE]

where $C_{L}$ is the computational cost of one sample with mesh size $h_{L}$. $T_{c}^{SL}$ would be very high when $N^{SL}_{L}$ and $C_{L}$ are both very large. By the accuracy formulation (34) of SLMC method, the sampling error and the FEM error should be both small enough, if a small mean squared error is required. Thus $N^{SL}_{L}$ should be larger while the mesh size $h_{L}$ becomes smaller. On the other hand, $C_{L}$ increase exponentially as the mesh size $h_{L}$ becomes smaller. Thus the total computational cost increases very fast as mesh size $h_{L}$ 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

[TABLE]

Then we can use the hierarchical meshes to construct the MLMC method to generate the expect value of $Q$. Let $\{\mathcal{T}_{\ell}\}_{\ell=0}^{L}$ be a sequence of quasi-uniform triangulation meshes with the mesh sizes $\{h_{\ell}\}_{\ell=0}^{L}$. These mesh sizes satisfy $h_{\ell}=h_{0}c_{h}^{-\ell}$, $\ell=0,1,2,\cdots,L$. And $\{N^{ML}_{\ell}\}_{\ell=0}^{L}$ are the numbers of samples with the mesh sizes $\{h_{\ell}\}_{\ell=0}^{L}$. Then the approximation $\hat{Q}_{L}^{ML}$ of the expected value by the MLMC method is given by:

[TABLE]

and the corresponding mean squared error of the MLCM method with norm $\|\cdot\|$ is

[TABLE]

For simplicity, let $Q_{-1}=0$, $h_{-1}=0$, $v_{\ell}=\|\mathbb{V}[Q_{\ell}-Q_{\ell-1}]\|$, $\ell=0,1,2,\cdots,L$, and $C_{\ell}$ be the computational cost of generating one sample of $Q_{\ell}-Q_{\ell-1}$, $\ell=0,1,\cdots,L$. Then the mean squared error is rewrote as

[TABLE]

And the total computational cost $T^{ML}_{c}$ is

[TABLE]

By the mean squared error of SLMC method (34) and MLMC method (39), the accuracy of approximation of expected value is based on two parts, i.e., the sampling error and FEM error. The FEM error $\|(\mathbb{E}[Q_{L}]-\mathbb{E}[Q])^{2}\|$ is fixed when the mesh size $h_{L}$ is given. Thus the sampling error should be small enough with the given mesh size $h_{L}$. We substitute the sampling errors in SLMC method and MLMC method by:

[TABLE]

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

[TABLE]

Then we show our strategy to generate the key parameters for MLMC method: the total number of levels $L$, and the number of samples at every level $\{N^{ML}_{\ell}\}_{\ell=0}^{L}$.

The total number of levels $L=\log_{c_{h}}^{h_{0}/h_{L}}$ depends on three variables: the mesh size decrease parameter $c_{h}$, the largest mesh size $h_{0}$ and the smallest mesh size $h_{L}$. The largest mesh size $h_{0}$ is constrained by the size of the physical area. The smallest mesh size $h_{L}$ depends on the accuracy of FEM as the practical problem required. Then $L$ is given after the setting $c_{h}=2$.

The guideline in designing the number of samples at every level is minimizing computational cost under the given sampling error. Thus we introduce the optimization problem as follow:

[TABLE]

This optimization problem is solved by the method of Lagrangian multipliers:

[TABLE]

Then the equations for $\{N_{\ell}\}_{\ell=0}^{L}$ are

[TABLE]

Then the number of samples at the every level is

[TABLE]

and the optimal computational cost is

[TABLE]

In the application of this strategy, the parameter $e^{ML}_{L}$ is given by the formula (42).

We assume $v_{\ell}=\mathcal{O}(h_{\ell}^{\beta})$ by the virtue of experience, and $C_{\ell}=\mathcal{O}(h_{\ell}^{-\gamma})$ because the number of information be calculated increase exponentially while the mesh size becomes smaller. Under the choice $c_{h}=2$, i.e., $h_{\ell}=h_{0}2^{-\ell}$, $\ell=0,1,\cdots,L$, by the formula (46), for any $j>i$

[TABLE]

Thus the number of samples be calculated becomes smaller while mesh size becomes smaller. The decrease of $N_{\ell}$ is the reason why the MLMC method can reduce the total computational cost.

Since the computational cost of every sample with the mesh size $h_{0}$ is low, $v_{0}$ is easy to calculate by Monte Carlo method with low computational cost. Then $v_{\ell}$, $\ell=1,\cdots,L$, can be given by $v_{\ell}=\mathcal{O}(h_{\ell}^{\beta})$ with the corresponding parameter $\beta$. Thus how to determine parameter $\beta$ is a key problem for MLMC method. Our strategy is provided in the following section.

4.4 Multi-grid methods

The total computational cost depends on the number of samples and the computational cost of every sample. Since we have reduced the total computational cost in probability space by using the MLMC method to reduce the number of samples, it is a heuristic problem that can we also reduce the total computational cost in physical space by reducing computational cost of every sample. Inspired by the hierarchical meshes used in the MLMC method, we adopt the multi-grid (MG) method to reduce the computational cost in physical space.

In the physical space, the finite element method (FEM) is chosen to construct the discrete form of weak formulation (16) under the given samples of hydraulic conductivity. We adopt the Taylor-Hood element in the conduit domain, and the quadratic element in the porous media domain. Then for every given sample of hydraulic conductivity $\mathbb{K}(\omega,x)$, the weak formulation (16) is discretized into the following matrix-vector form

[TABLE]

where $A_{m}$ is the discretization of $g\int_{D_{m}}(\mathbb{K}\nabla\phi_{m})\cdot\nabla\psi_{m}dx$, $B_{1}$ is the discretization of $-g\int_{\Gamma_{I}}(\vec{u}_{s}\cdot\vec{n}_{s})\psi_{m}d\Gamma_{I}$, $B_{2}$ is the discretization of $\int_{\Gamma_{I}}g\phi_{m}\vec{v}_{s}\cdot\vec{n}_{s}+\frac{\alpha\nu\sqrt{d}}{\sqrt{\text{trace}(\Pi)}}P_{\tau}(\mathbb{K}\nabla\phi_{m})\cdot\vec{v}_{s}d\Gamma_{I}$, $A_{s}$ is the discretization of $\int_{D_{s}}2\nu\mathbb{D}(\vec{u}_{s}):\mathbb{D}(\vec{v}_{s})dx+\int_{\Gamma_{I}}\frac{\alpha\nu\sqrt{d}}{\sqrt{\text{trace}(\Pi)}}P_{\tau}(\vec{u}_{s})\cdot\vec{v}_{s}d\Gamma_{I}$, $B^{\prime}_{p}$ is the discretization of $\int_{D_{s}}p_{s}\nabla\cdot\vec{v}_{s}dx$, $b_{m}$ is the discretization of $g\int_{D_{m}}f_{m}\psi_{m}dx$, and $\vec{b}_{s}$ is the discretization of $\int_{D_{s}}\vec{f}_{s}\cdot\vec{v}_{s}dx+\int_{\Gamma_{I}}gz\vec{v}_{s}\cdot\vec{n}_{s}d\Gamma_{I}$.

Since there exists a zeros block in the diagonal of stiffness matrix $\mathbf{L}$, we can not directly solve the algebraic equations $\eqref{algebraic_equations}$ by iterative method such as Gauss-Seidel method. Inspired by the multi-grid method for Stokes equations, we adopt the efficient least squares commutator distributive Gauss-Seidel (LSC-DGS) relaxation [87, 88] in this paper. The right-side operator $\mathbf{M}$ is given as:

[TABLE]

Multiplying $\mathbf{L}$ with $\mathbf{M}$ yields

[TABLE]

By $\mathbf{S}:=\mathbf{L}\mathbf{M}$ and $\mathbf{y}:=\mathbf{M^{-1}x}$, the equivalent algebraic equations are given as

[TABLE]

The standard Gauss-Seidel method is proposed to solve the equivalent algebraic equations (51). And the following $\mathcal{V}$-cycle multi-grid method is applied to reduce the computational cost in physical space. As same as in the MLMC method, the hierarchical quasi-uniform triangulation meshes are $\mathcal{T}_{\ell}$ with the mesh sizes $h_{\ell}=h_{0}c_{h}^{-\ell}$, $\ell=0,1,2,\cdots,L$. Then the $\mathcal{V}$-cycle multi-grid method on the mesh $\mathcal{T}_{\ell}$ with the mesh size $h_{\ell}$ is given as:

Algorithm 6**.**

$\mathbf{y}\leftarrow$* $\mathcal{V}$-cycle$(\mathbf{S},\mathbf{b},\ell)$*

1. (1)

Relax $\lambda_{1}$ times on the fine mesh $h=h_{\ell}$ with the initial gauss $\mathbf{y}$ to reach $\mathbf{y}^{h}$. 2. (2)

Obtain the residual on the fine mesh as

[TABLE]

and restrict the residual from the fine mesh $h$ to the coarse mesh $H=h_{\ell-1}$ by $r^{H}=\mathcal{R}^{H}_{h}r^{h}$, where $\mathcal{R}^{H}_{h}$ is the restriction matrix. 3. (3)

Solve the corrected error from the residual equation on the coarse mesh $H$:

• •

If $\ell=1$, use a direct or fast iterative method to solve $\mathbf{S}^{H}e^{H}=r^{H}$;

• •

If $\ell>1$, use the $\ell$-grid method to solve $\mathbf{S}^{H}e^{H}=r^{H}$ from a zero initial gauss on the mesh $\mathcal{T}_{\ell-1}$ by $e^{H}\leftarrow$ $\mathcal{V}$-cycle$(\mathbf{S}^{H},r^{H},\ell-1)$;

where $\mathbf{S}^{H}$ is the approximation of $\mathbf{S}$ on the coarse mesh. 4. (4)

Prolongate the corrected error form coarse mesh $H$ to the fine mesh $h$ by $e^{h}=\mathcal{I}^{h}_{H}e^{H}$, where $\mathcal{I}^{h}_{H}$ is the interpolation matrix. And correct the approximation by

[TABLE] 5. (5)

Relax $\lambda_{2}$ times on the fine mesh $h$ with the initial gauss $\mathbf{y}^{new}$.

We can replace the $\mathcal{V}$-cycle by $\mathcal{W}$-cycle or $\mathcal{F}$-cycle. Furthermore, the solutions be calculated on the coarse mesh in MLMC method could be used as the initial gauss on the fine mesh in MG method to further reduce the computational cost. Then the following multi-gird multi-level Monte Carlo (MGMLMC) method is developed to reduce the computational cost both in the probability space and in physical space.

Algorithm 7**.**

multi-grid multi-level Monte Carlo method

1. (1)

On the mesh grid $\mathcal{T}_{0}$ with the mesh size $h_{0}$, for the 1st to the $N^{ML}_{0}$ sample of hydraulic conductivity $\mathbb{K}(\omega,x)$, solve the numerical approximations $Q_{0}^{i}$ by standard Gauss-Seidel with the initial gauss $\mathbf{0}$, $i=1,2,\cdots,N_{0}^{ML}$; 2. (2)

On the mesh grid $\mathcal{T}_{1}$ with the mesh size $h_{1}$, for the 1st to the $N^{ML}_{1}$ sample of hydraulic conductivity $\mathbb{K}(\omega,x)$, solve the numerical approximations $Q_{1}^{i}$ by $\mathcal{V}$-cycle $2$-grid method with the initial gauss $\mathcal{I}^{h}_{H}Q_{0}^{i}$, $i=1,2,\cdots,N_{1}^{ML}$, where $h=h_{1}$ and $H=h_{0}$; 3. (3)

$\cdots\cdots$** 4. (4)

On the mesh grid $\mathcal{T}_{\ell}$ with the mesh size $h_{\ell}$, for the 1st to the $N^{ML}_{\ell}$ sample of hydraulic conductivity $\mathbb{K}(\omega,x)$, solve the numerical approximations $Q_{\ell}^{i}$ by $\mathcal{V}$-cycle $(\ell+1)$-grid method with the initial gauss $\mathcal{I}^{h}_{H}Q_{\ell-1}^{i}$, $i=1,2,\cdots,N_{\ell}^{ML}$, where $h=h_{\ell}$ and $H=h_{\ell-1}$; 5. (5)

$\cdots\cdots$** 6. (6)

On the mesh grid $\mathcal{T}_{L}$ with the mesh size $h_{L}$, for the 1st to the $N^{ML}_{L}$ sample of hydraulic conductivity $\mathbb{K}(\omega,x)$, solve the numerical approximations $Q_{L}^{i}$ by $\mathcal{V}$-cycle $(L+1)$-grid method with the initial gauss $\mathcal{I}^{h}_{H}Q_{L-1}^{i}$, $i=1,2,\cdots,N_{L}^{ML}$, where $h=h_{L}$ and $H=h_{L-1}$.

4.5 Computational cost of SLMC method and MGML method

The numerical error is estimated as follow:

Proposition 8**.**

Under the assumption (25) or (26), the solutions of problem (16) satisfy the following discrete error estimate

[TABLE]

where $h_{\ell}$ is the mesh size of the given quasi-uniform triangulation mesh $\mathcal{T}_{\ell}$, $r=0,1$, and $\|\cdot\|_{X^{r}}$ is the norm of $\underline{u}$ defined in (15).

Proof.

Based on the analysis in [16, 17, 62], we have

[TABLE]

with $r=0,1$. Then the assertion follows with the above conclusion by the the definition of the norm $\|\cdot\|_{X^{r}}$ in (15). ∎

Then the numerical errors of SLMC method and MGML method are bounded by the mesh size $h_{\ell}$ and the number of samples $N_{\ell}$.

Lemma 9**.**

Under the assumption (25) or (26), the error bounds of SLMC method (32) and MLMC (37) for the problem (16) are given as follows

[TABLE]

where $C$ depends on $\underline{u}$ and $p_{s}$, $h_{\ell}$ is the mesh size of the quasi-uniform triangulation mesh $\mathcal{T}_{\ell}$, and $\hat{\underline{u}}^{SL}_{\ell}$, $\hat{p}^{SL}_{s,\ell}$, $\hat{\underline{u}}^{ML}_{\ell}$, $\hat{p}^{ML}_{s,\ell}$ are the approximations of expect value by SLMC method, MLMC method, $\ell=0,1,\cdots,L$.

Proof.

For simplification, let the symbol $Q$ substitute the variables $\phi_{m}$, $u_{s}$ or $p_{s}$. And let $\mathcal{L}(V)$ denote the corresponding space of $Q$. Thus $\mathcal{L}(V)$ may denote $\mathcal{H}^{1}(D_{m})$, $\mathcal{L}^{2}(D_{s})$ or $\vec{\mathcal{H}}^{1}(D_{s})$, i.e., $V$ may be $H^{1}(D_{m})$, $L^{2}(D_{s})$ or $\mathbb{H}^{1}(D_{s})$, which depends on the choice of $Q$. Then we can analyse the error of the approximation of expect value of $\phi_{m}$, $u_{s}$ or $p_{s}$ by analysing $\mathbb{E}[Q]-\hat{Q}^{SL}_{\ell}$ with the norm $\|\cdot\|_{\mathcal{L}(V)}$ as follow.

[TABLE]

For $\|\mathbb{E}[Q]-\mathbb{E}[Q_{\ell}]\|_{\mathcal{L}(V)}$, we have

[TABLE]

For $\|\mathbb{E}[Q_{\ell}]-\hat{Q}^{SL}_{\ell}\|_{\mathcal{L}(V)}$, we have

[TABLE]

The last inequality is based on $\mathbb{E}[(\mathbb{E}[Q_{\ell}]-Q_{\ell})^{2}]=\mathbb{E}[(Q_{\ell})^{2}]-(\mathbb{E}[Q_{\ell}])^{2}\leq\mathbb{E}[(Q_{\ell})^{2}]$.

Thus we obtain

[TABLE]

Then by the Proposition 8, we have

[TABLE]

where $C(\underline{u},p_{s})$ depends on $\|\underline{u}\|_{X^{0}}$, $\|\underline{u}\|_{X^{1}}$ and $\|p_{s}\|_{\mathcal{L}^{2}(D_{s})}$.

Because the idea to prove the assertion of MLMC method is as same as that in the proof of the assertion of SLMC method, we skip it. ∎

By equilibrating the sampling error in probability space and the FEM error in physical space, we have the following two conclusions based on the conclusions in Lemma 9.

[TABLE]

The formula (59) is the relationship between the numbers of samples $N_{L}$ and the mesh sizes $h_{L}$ in the SLMC method, which is based on the conclusion (53). And the formula (59) is the relationship between the numbers of samples $\{N_{\ell}\}_{\ell=0}^{L}$ and the mesh sizes $\{h_{\ell}\}_{\ell=0}^{L}$ in the MLMC method, which is based on the conclusion (54).

In the SLMC method, by (59), it is easy to see that the number of samples $N_{L}$ on the finest mesh is determined by the mesh size $h_{L}$, then the computational cost is distinct. In the MLMC method, the number of samples $N_{\ell}$ on every level is determined by the formula (46), since the sampling error is bounded by the formula (60), then the computational cost is also distinct.

Theorem 10**.**

Under the assumption (25) or (26), for the problem (16), if we choose the SLMC method (32) on the triangulation mesh $\mathcal{T}_{L}$ with the mesh size $h_{L}$, or the MGMLMC method Algorithm 7 on the hierarchical quasi-uniform triangulation meshes $\{\mathcal{T}_{\ell}\}_{\ell=0}^{L}$ with mesh sizes $h_{\ell}=h_{0}2^{-\ell},\ell=0,1,\cdots,L$ to solve the approximations of expect value, we can evaluate the computational cost as follows:

[TABLE]

where $M_{L}$ is the number of information be calculated for one sample on the mesh $\mathcal{T}_{L}$ with the mesh size $h_{L}$, $d$ is the dimension of the physical space, and $\beta$ is the decrease rate of the variance. Furthermore, the ratio of the computational cost of SLMC method and MGMLMC method is given as

[TABLE]

Proof.

Under the assumption M_{\ell}=\mathcal{O}\big{(}h_{\ell}^{-d}\big{)}, M_{\ell}=\mathcal{O}\big{(}M_{L}2^{(\ell-L)d}\big{)} is given by the setting $h_{\ell}=h_{0}2^{-\ell}$. Since the standard Gauss-Seidel method is chosen to solve the algebraic equations (51) in SLMC method, the computational cost $C_{\ell}$ with the mesh size $h_{\ell}$ is C_{\ell}=\mathcal{O}\big{(}M^{2}_{L}2^{2(\ell-L)d}\big{)}.

For the SLMC method, the number of samples on mesh $\mathcal{T}_{L}$ with mesh size $h_{L}$ is N_{L}^{SL}=\mathcal{O}\big{(}h_{L}^{-2}\big{)}=\mathcal{O}\big{(}M_{L}^{2/d}\big{)}, by the bound of sampling error in (59). Then the computational cost of SLMC method is T^{SL}_{c}=N_{L}^{SL}C_{L}=\mathcal{O}\Big{(}M_{L}^{2+\frac{2}{d}}\Big{)}.

For the MLMC method, the bound of sampling error is e_{L}^{ML}=\mathcal{O}\big{(}h_{L}^{2}\big{)}=\mathcal{O}\big{(}M_{L}^{-2/d}\big{)} by (60). Then by (46) the number of samples at the initial level is

[TABLE]

And by the relationship between the numbers of samples at two different levels (48), the number of samples on mesh $\mathcal{T}_{\ell}$ with mesh size $h_{\ell}$ is

[TABLE]

In this paper, we adopt the $\mathcal{V}$-cycle multi-grid methods. Then the computational cost [11, 82] on mesh $\mathcal{T}_{\ell}$ with mesh size $h_{\ell}$ is

[TABLE]

Then the computational cost of MGMLMC method is

[TABLE]

By $h_{L}=h_{0}2^{-L}$ and M_{L}=\mathcal{O}\big{(}h_{L}^{-d}\big{)}, we have

[TABLE]

∎

5 Numerical experiments

In this section, we use numerical experiments to demonstrate both the features of the MGMLMC method and the theoretical conclusion. The first part is to generate the realizations of random hydraulic conductivity $\mathbb{K}$ by the grid based method. The second part is to determine the parameters $\beta$, which will be used to calculate the $\{N^{ML}_{\ell}\}_{\ell=0}^{L}$ in MLMC method. The last part is to provide the numerical results in detail.

We assume that the domain $D_{ms}$ consists of two rectangles, the upper rectangle is the porous media domain $D_{m}=(0,1)\times(0,0.75)$, and the other rectangle is the conduit domain $D_{s}=(0,1)\times(-0.25,0)$, shown as the Figure 1. The whole domain $D_{ms}=D_{m}\cup D_{s}$ with the interface $\Gamma_{I}=(0,1)\times\{0\}$. The boundary are $\Gamma_{m}=\{0,1\}\times(0,0.75)\cup(0,1)\times\{0.75\}$ and $\Gamma_{s}=\Gamma_{s_{1}}\cup\Gamma_{s_{2}}\cup\Gamma_{s_{3}}$, where $\Gamma_{s_{1}}=\{0\}\times(-0.25,0),\Gamma_{s_{2}}=(0,1)\times\{-0.25\},\Gamma_{s_{3}}=\{1\}\times(-0.25,0)$. For simplicity, let $g=1$, $z=0$, $\alpha=1$, $\nu=1$ and $\mathbb{K}(\omega,x)=e^{Z(\omega,x)}\mathbb{I}$. The covariance function of $Z$ is r(x,y)=r\big{(}(x_{1},x_{2}),(y_{1},y_{2})\big{)}=0.1e^{-\frac{|x_{1}-y_{1}|}{0.2}-\frac{|x_{2}-y_{2}|}{0.2}}.

5.1 The realizations of random hydraulic conductivity

Because the diagonal matrix is given as $\mathbb{K}(\omega,x)=e^{Z(\omega,x)}\mathbb{I}$, it is natural to generate the realizations of $K(\omega,x)=e^{Z(\omega,x)}$, and then copy the realizations $d$ times to construct $\mathbb{K}(\omega,x)$. As the hierarchical meshes are used in MLMC method, the realizations of $K(\omega,x)=e^{Z(\omega,x)}$ could be first generated on the finest mesh by the grid based method. And then the realizations on the coarse mesh can be chosen as a subset of the realizations on finest mesh. Then the consistency of hydraulic conductivity $K(\omega,x)$ on every mesh could be ensured.

Because the Gauss quadrature points are the key points in the performance of finite element method, we calculate the value of the approximation of the hydraulic conductivity $K(\omega,x)$ on the Gauss quadrature points by the grid based method. In every triangle of the triangulation mesh $\mathcal{T}_{L}$, Gauss quadrature rule is applied with seven points and degree of precision three.

Four realizations of $K(\omega,x)$ are illustrated in Figure 2, by which the randomness of hydraulic conductivity is exhibited.

5.2 Determination of the parameter $\beta$ in MLMC method

As we have discussed in the analysis for the number of samples at every level, how to determine the parameter $\beta$ is a key problem in the numerical implementation of MLMC method. The $\beta$ can be approximated by using the variances $\{v_{\ell}\}_{\ell=0}^{L}$, which are calculated by SLMC method with mesh size $\{h_{\ell}\}_{\ell=0}^{L}$. Thus the computational cost of calculating the variances $\beta$ is more expensive than that of SLMC method with mesh size $h_{L}$, which contradicts the purpose of MLMC method. Then a practical method is needed to determine the parameter $\beta$ without calculating the variances $\{v_{\ell}\}_{\ell=0}^{L}$ by SLMC method at every level. In this paper, we develop the following method to calculate the parameter $\beta$.

Since the random hydraulic conductivity is only a parameter in the Darcy domain and on the interface, the Stokes equations in the coupled problem can be regarded as a boundary condition for a stochastic Darcy problem. We assume that the $\beta$ in the stochastic Darcy problem is an approximation for the $\beta$ in the stochastic Stokes-Darcy problem with the same random hydraulic conductivity. Compared with the coupled stochastic Stokes-Darcy problem, the domain of the stochastic Darcy problem is smaller and the computational cost of every sample is cheaper. Furthermore, we can also use the multi-grid method to reduce the computational cost in generating the approximation of $\beta$.

The stochastic Darcy problem is given by

[TABLE]

where $f(\omega,x)$ is a piecewise constant approximation of white noise, i.e.,

[TABLE]

Here $\sigma$ is a given constant, $V_{i}$ is the volume of non-overlapping tessellation $\{D_{i}\}_{i=1}^{I}$ as $D_{m}=\cup_{i=1}^{I}D_{i}$, $\chi_{i}(x)$ is the indicator function corresponding to $D_{i}$, and $\{X_{i}(\omega)\}_{i=1}^{I}$ is a given set of independent identically distributed standard Gaussian random variables.

Given the $\sigma$ and $\{X_{i}(\omega)\}_{i=1}^{I}$ in (70), it is easy to calculate $\beta$ with a given norm $\|\cdot\|$. For each $\sigma=0.02,0.8,1.2$, we choose 40 samples of $f$. For every given $\sigma$ and one sample of $f$, three $\beta$ are calculated with $\|\cdot\|_{L_{2}}$, $\|\cdot\|_{L_{\infty}}$ and $\|\cdot\|_{H_{1}}$ norms. The results of $\beta$ with each sample of $f$ and the choice of $\sigma$ are exhibited in Figure 3. And the mean values of $\beta$ with 40 samples of $f$ are shown in Table 1. One can see that the mean value of the $\beta$ changes only in a small range when $\sigma$ becomes larger. Thus the parameter $\beta$ is given as 2.02, 1.65, 1.30 w.r.t. $\|\cdot\|_{L_{2}}$, $\|\cdot\|_{L_{\infty}}$, $\|\cdot\|_{H_{1}}$ norms.

5.3 Main numerical results

Let $f_{m}=0$, $\vec{f}_{s}=0$, $\psi_{0}=0$, on $\Gamma_{m}$, $\vec{u}_{s}=(1,0)^{T}$, on $\Gamma_{s_{1}}$, $\vec{u}_{s}=(0,0)^{T}$, on $\Gamma_{s_{2}}$, and $\vec{u}_{s}=(1,0)^{T}$, on $\Gamma_{s_{3}}$.

For exhibiting the stochastic property of our problem, four samples of numerical solutions on the mesh $h_{L}=1/32$ with four different samples of $\mathbb{K}(\omega,x)$ are shown in the Figure 4.

For the hierarchical quasi-uniform triangulation mesh $\{\mathcal{T}_{\ell}\}_{\ell=0}^{L}$, four levels are chosen, i.e., $h_{\ell}=\frac{2^{-\ell}}{4},\ell=0,1,2,3$ with $h_{0}=1/4$. An explicit numerical method is needed to determine the parameter $\gamma$ in calculating $\{N^{ML}_{\ell}\}_{\ell=0}^{3}$, which is needed for MGMLMC method. Based on $C_{\ell}=\mathcal{O}(h_{\ell}^{-\gamma}),\ell=0,1,2,3$, we can compute $\gamma$ after the computational cost $\{C_{\ell}\}_{\ell=0}^{3}$ of a few samples at every level are recorded. The cpu time and tic-toc time with different mesh size are shown in the Table 2, and the corresponding $\gamma$ are 2.0536, 2.4549, which are illustrated in the Figure 5. In this paper, we choose $\gamma=2.4549$.

If the variance $v_{0}$ at the first level is known, the variance $\{v_{\ell}\}_{\ell=0}^{L}$ at the every level could be calculated by $v_{\ell}=\mathcal{O}(h_{\ell}^{\beta})$, while the parameter $\beta$ is approximated by the $\beta$ of stochastic Darcy problem. The variance $v_{0}$ at the first level is easy to be calculated with low computational cost. Then using the formula (46) with the parameters $\beta$ and $\gamma$ we have gained, we can obtain the number of samples at every level based on the optimization problem (43) with the given sampling error $e_{L}$. The numbers of samples on every level with given sampling error $e_{L}$ are shown in the Table (3).

To verify the accuracy of MGMLMC method, the relative errors between solutions of SLMC method and those of MGMLMC method are shown in the Table 4, and the numerical approximations of expectation of speed on mesh $h_{L}=1/32$ are compared with those two methods in the Figure 6. To illustrate the efficiency of MLMC method and MGMLMC method, the computational cost of SLMC method, MLMC method and MGMLMC method are compared in the Table 5. Based on these results, it is easy to see that the MGMLMC method significantly reduce the computational cost with the same accuracy as SLMC method.

6 Conclusion

In this paper, for the stochastic Stokes-Darcy interface problem, we proved the well-posedness of weak solution, and developed an accurate and efficient multi-grid multi-level Monte Carlo method to solve the numerical approximations. In the proof of the well-posedness, we overcame the difficulties caused by the random hydraulic conductivity both in the porous media domain and on the interface. For the MLMC method, we provided a strategy to calculate the number of samples on every level. We verified the features of the numerical method and the theoretical conclusions.