A local-global multiscale mortar mixed finite element method for multiphase transport in heterogeneous media

TL;DR

This paper introduces a multiscale mortar mixed finite element method that efficiently models multiphase flow in heterogeneous media, combining multiscale basis functions with explicit finite volume methods for improved accuracy and computational performance.

Contribution

The paper presents a novel local-global multiscale mortar mixed finite element method that integrates multiscale basis functions with explicit finite volume schemes for multiphase transport.

Findings

Accurate pressure and velocity fields in heterogeneous media

Effective multiscale basis construction from global pressure

Validated performance on 2D and 3D benchmark models

Abstract

In this paper, we propose a local-global multiscale mortar mixed finite element method (MMMFEM) for multiphase transport in heterogeneous media. We consider the two-phase flow system, the pressure equation is solved via the multiscale mortar mixed finite element method, a mass conservation velocity field can be obtained, then we use explicit finite volume method to solve the saturation equation. We use polynomials and multiscale basis to form the coarse mortar space. The multiscale basis is the restriction of global pressure field obtained at previous time step on the coarse interface. We solve the pressure equation on the fine grid to initialize the simulation. Numerical experiments on some benchmark 2D and 3D heterogeneous models are provided to validate the performance of our method.

Figures

Tables

Table 1. Table 1: Accuracy and efficiency comparisons between MMMFEM and fine-scale simulations, model 2.

Simulator	Basis type [$Nb$]	$n$	Dim	$e_S$	CPU time (hours)
MMFEM	-	-	1209600	-	0.93
MMMFEM	$P_0$+multiscale [2]	10	24192	0.0495	0.29
MMMFEM	$P_1$+multiscale [4]	10	48384	0.0222	0.34
MMMFEM	$P_0$+multiscale [2]	5	96768	0.0246	0.32
MMMFEM	$P_1$+multiscale [4]	5	193536	0.0069	0.37

Table 2. Table 2: Accuracy and efficiency comparisons between MMMFEM and fine-scale simulations, model 3.

Simulator	Basis type [$Nb$]	$n$	Dim	$e_S$	CPU time (hours)
MMFEM	-	-	2007200	-	3.51
MMMFEM	$P_0$+multiscale [2]	10	40144	0.1139	0.55
MMMFEM	$P_1$+multiscale [4]	10	80288	0.0837	0.63
MMMFEM	$P_0$+multiscale [2]	5	160576	0.0270	0.70
MMMFEM	$P_1$+multiscale [4]	5	321152	0.0127	1.13

Table 3. Table 3: Accuracy and efficiency comparisons between MMMFEM and fine-scale simulations, model 2, source type II.

Simulator	Basis type [$Nb$]	$n$	Dim	$e_S$	CPU time (hours)
MMFEM	-	-	1209600	-	0.93
MMMFEM	$P_0$+multiscale [2]	10	24192	0.0489	0.30
MMMFEM	$P_1$+multiscale [4]	10	48384	0.0204	0.35
MMMFEM	$P_0$+multiscale [2]	5	96768	0.0289	0.33
MMMFEM	$P_1$+multiscale [4]	5	193536	0.0082	0.38

Table 4. Table 4: Accuracy and efficiency comparisons between MMMFEM and fine-scale simulations, model 3, source type II.

Simulator	Basis type [$Nb$]	$n$	Dim	$e_S$	CPU time (hours)
MMFEM	-	-	2007200	-	3.58
MMMFEM	$P_0$+multiscale [2]	10	40144	0.0474	0.57
MMMFEM	$P_1$+multiscale [4]	10	80288	0.0303	0.64
MMMFEM	$P_0$+multiscale [2]	5	160576	0.0321	0.69
MMMFEM	$P_1$+multiscale [4]	5	321152	0.0179	1.09

Equations

$$\varphi \frac{\partial S_l}{\partial t}+\nabla \cdot {\mathbf{u}}_l=q_l,$$

$${\mathbf{u}}_l=-{\lambda }_l(S_l)K\nabla p_l.$$

$$\mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\varphi \frac{\partial S_w}{\partial t}+\nabla \cdot ({f_w\mathbf{u}}_l)=q_w,$$

$$\mathbf{u}+\lambda K\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{u}\cdot \mathbf{n}$$

$${\mathbf{V}}_i=H(\text{div};K_i),\mathbf{V}=\underset{i=1}{\overset{N}{⨁}}{\mathbf{V}}_i,W_i=L_2(K_i),W=\underset{i=1}{\overset{N}{⨁}}{W}_i,M_i=H^{1/2}(E_i),M=\underset{i=1}{\overset{|{\mathcal{E}}_H|}{⨁}}{M}_i,$$

$$((\lambda K{)}^{-1}\mathbf{u},\mathbf{v}{)}_{K_i}-(p,\nabla \cdot \mathbf{v}{)}_{K_i}+{⟨\Lambda ,\mathbf{v}\cdot {\mathbf{n}}_i⟩}_{\partial K_i}$$

$$(\nabla \cdot \mathbf{u},w{)}_{K_i}$$

$$\sum _{i=1}^N{⟨\mathbf{u}\cdot {\mathbf{n}}_i,\mu ⟩}_{\partial K_i}$$

$$((\lambda K{)}^{-1}{\mathbf{u}}_h,{\mathbf{v}}_h{)}_{K_i}-(p_h,\nabla \cdot {\mathbf{v}}_h{)}_{K_i}+{⟨{\Lambda }_H,{\mathbf{v}}_h\cdot {\mathbf{n}}_i⟩}_{\partial K_i}$$

$$(\nabla \cdot {\mathbf{u}}_h,w_h{)}_{K_i}$$

$$\sum _{i=1}^N{⟨{\mathbf{u}}_h\cdot {\mathbf{n}}_i,{\mu }_H⟩}_{\partial K_i}$$

$$a_{H,i}(\Lambda ,\mu )=-{⟨{\mathbf{u}}_h^{\ast }(\Lambda )\cdot {\mathbf{n}}_i,\mu ⟩}_{\partial K_i},$$

$$a_H=\sum _{i=1}^N{a}_{H,i}(\Lambda ,\mu ),$$

$$\displaystyle\big{(}(\lambda{K})^{-1}\boldsymbol{u}_{h}^{\ast}(\Lambda),\boldsymbol{v}_{h}\big{)}_{K_{i}}-\big{(}p_{h}^{\ast}(\Lambda),\nabla\cdot\boldsymbol{v}_{h}\big{)}_{K_{i}}$$

$$\displaystyle\big{(}\nabla\cdot\boldsymbol{u}_{h}^{\ast}(\Lambda),w_{h}\big{)}_{K_{i}}$$

$$g_{H,i}(\mu )={⟨{\bar{\mathbf{u}}}_h\cdot {\mathbf{n}}_i,\mu ⟩}_{\partial K_i},$$

$$g_H(\mu )=\sum _{i=1}^N{g}_{H,i}(\mu ),$$

$$\displaystyle\big{(}(\lambda K)^{-1}\boldsymbol{\bar{u}}_{h},\boldsymbol{v}_{h}\big{)}_{K_{i}}-\big{(}\bar{p}_{h},\nabla\cdot\boldsymbol{v}_{h}\big{)}_{K_{i}}$$

$$\displaystyle\big{(}\nabla\cdot\boldsymbol{\bar{u}}_{h},w_{h}\big{)}_{K_{i}}$$

$$a_H({\Lambda }_H,\mu )=g_H(\mu ),\forall \text{leavevmode}\mu \in M_H.$$

$${\mathbf{u}}_h={\mathbf{u}}_h^{\ast }({\Lambda }_H)+{\bar{\mathbf{u}}}_h,p_h={\tilde{p}}_h-\frac{1}{|\Omega |}{\int }_{\Omega }{\tilde{p}}_h,$$

$$e_S(i):=\frac{\Vert S_{ms}(i)-S_{ref}(i){\Vert }_{(L^2,\Omega )}}{\Vert S_{ref}(i){\Vert }_{(L^2,\Omega )}}.$$

Full text

1 Introduction

Modeling flow in porous media has many important practical applications, such as nuclear waste storage, and oil and natural gas production. However, advanced modern reservoir characterization techniques can generate very detailed reservoir models with multiple scales, from core scales (centimeters) to geological scales (kilometers). Therefore, direct simulation is inevitably expensive and model reduction techniques such as upscaling and multiscale methods are quite necessary. The basis idea of upscaling [13, 26] is to homogenize the media properties based on some rules so that the problem can be solved in a reduced model. In multiscale methods [15, 17, 5, 12, 11, 22, 25, 7, 4, 20], one solves the problems on a coarse grid but with precomputed multiscale basis functions that includes small scale information of the media.

Among different multiscale methods, the mixed multiscale finite element method [11] and the multiscale finite volume method [22] are especially useful for the flow simulation since these methods are mass conservative. In [1], the author used the restriction of global velocity field in coarse interface as boundary condition to compute the multiscale basis functions. This idea (using limited global information) is extended to the framework of finite volume [16]. In [3, 10, 14], improvements such as adaptive update the multiscale basis functions are developed. Our goal here is to apply above mentioned ideas to the mortar mixed finite element method (MMFEM) [19, 6], which is a modification of the mixed finite element methods by introducing a Lagrange multiplier to impose the continuity of flux. This modification brings some additional advantages, such as more efficient implementation since it yields a a symmetric and positive definite bilinear form that defined only on the interfaces of the grid. Besides, it also allows nonconforming grid discretization while keeping the mass conservation property.

The multiscale mortar mixed finite element method (MMMFEM) was proposed in [7, 8] to reduce the dimension of the mortar mixed finite element method with the sacrifices of some accuracy. In MMMFEM, the construction of the mortar space that defined on coarse interfaces is cubical. Polynomials [7] or multiscale basis functions [8, 28] are used to form the mortar space. However, polynomials are only sufficient for very smooth media, while multiscale functions lack global information which fails to generate satisfiable flux for accurate simulation of coupled flow and transport problems if the media is highly heterogeneous[27]. Developing efficient domain decomposition preconditioners [8] to solve the fine scale problem is a potential way to tackle this kind of problem. However, if one can equip the mortar space with basis functions that can capture global information with the evolving of time, then we can avoid solving the fine problems. In [29], the authors developed residual based online multiscale basis functions which contain global information and successfully applied this method to two phase simulations, however, in some cases, one needs to iterate several times to get accurate solution, which makes the mortar space too enriched. Therefore, further reduction is necessary. We adopt the above mentioned ideas of using limited global information to mortar finite element. However, naive extension, using the restriction of global pressure on coarse interface to form the mortar space fails. Our strategy here is fill the mortar space with fixed local basis functions such as polynomials and the restriction of global pressure field obtained at previous time step on the coarse interface. Therefore, our method is a local-global method. To initialize the simulation, we first need to compute an accurate initial global pressure field, which can be achieved by solving the fine problem with preconditioners [8] or online multiscale method [29]. Numerical results show that with only 2 bases in 2D case and 4 bases in 3D cases, the MMMFEM can generate very accurate saturation profiles with much reduced computational cost.

The paper is organized as follows. In Section 2, we introduce the coupled two-phase flow system. In Section 3, we present MMMFEM method for the pressure equation. The construction of the local-global multiscale mortar space in introduced in Section 4. Numerical examples are given in Section 5, and conclusions are made in the last section.

2 Mathematical model

We consider a model problem for immiscible and incompressible two-phase flow in a reservoir domain (denoted by $\Omega$) with the assumption that the fluid displacement is driven by viscous effects, that is, we neglect compressibility and gravity for simplicity. We consider water and oil phases which are assumed to be immiscible. By the mass conservation law, we have following equations for each phase $l$:

[TABLE]

where the phase velocities ${\bf u}_{l}$ are given by Darcy’s law

[TABLE]

Here $\phi$ is the porosity, $S_{l}$ is the $l$-phase saturation (fraction of the void occupied by phase $l$), and $q_{l}$ is a source (or sink) term, ${K}$ is the permeability tensor which can be highly heterogeneous, $p_{l}$ is the phase pressure, and $\lambda_{l}$ is the phase mobility given by $\lambda_{l}(S_{l})=\frac{k_{rl}(S_{l})}{\mu_{l}}$ , where $k_{rl}$ and $\mu_{l}$ are the relative permeability and viscosity of phase $l$, respectively.

Since we neglect effects from capillary pressure so that $\nabla p_{o}=\nabla p_{w}$, and we denote $p_{o}=p_{w}=p$. Then we can derive the following pressure equation

[TABLE]

where $\bf u=\bf u_{o}+\bf u_{w}$ and $q=q_{o}+q_{w}$

We assume that the two phases occupy the void space completely, therefore $S_{o}+S_{w}=1$, let $\lambda(S)=\lambda_{w}(S_{w})+\lambda_{o}(1-S_{w})$ be the total mobility and the water fractional flow $f_{w}=\lambda_{w}/\lambda$, then we have the conservation equation for water as follows

[TABLE]

which is henceforth called the saturation equation, In the following we will, drop the w-subscripts of $S_{w}$ for ease of notation, We assume no flow at boundary and the initial saturation is given. Note the pressure equation and the saturation equation are coupled through the total mobility. Here, we follow the sequential formulation [1] to solve the above coupled system, that is we use saturation at the previous time step to solve the pressure equation to obtain the velocity $\bf u$. With the updated velocity, we solve the transport equation with explicit finite volume method and obtain new saturation, which is then used to solve the pressure equation again and the process is continued until the final time is reached.

3 Multiscale mortar mixed finite element method

In this section, we will introduce the framework of multiscale mortar mixed finite element method for the pressure equation. To fix our attention, We consider

[TABLE]

where $\Omega\subset\mathbb{R}^{d}(d=2,3)$ is a bounded polyhedral domain with outward unit normal vector $\boldsymbol{n}$ on the boundary, $q\in L^{2}(\Omega)$. $K(x)$ is a symmetric and uniformly positive definite permeability tensor with components in $L^{\infty}(\Omega)$.

3.1 Formulation of the method

We divide $\Omega$ into non-overlapping polygonal coarse elements $K_{i}$ so that $\overline{\Omega}=\cup_{i=1}^{N}\overline{K}_{i}$, where $N$ is the number of coarse elements. The decomposition of the domain can be nonconforming. We call $E_{H}$ a coarse interface of the coarse block $K_{i}$ if $E_{H}=\partial K_{i}\cap\partial K_{j}$ or $E_{H}=\partial K_{i}\cap\partial{\Omega}$. Let $\mathcal{E}_{H}(K_{i})$ be the set of all coarse interfaces on the boundary of the coarse element $K_{i}$, and $\mathcal{E}_{H}=\cup_{i=1}^{N}\mathcal{E}_{H}(K_{i})$ be the set of all coarse interfaces. For a coarse element $K_{i}$ and coarse interface $E_{i}$, let

[TABLE]

Then, a weak form of system (6a)-(6c) reads: find ${\bf u}\in{\bf V}$, $p\in W$ and $\Lambda\in M$ such that for each $1\leq i\leq N$,

[TABLE]

where $\boldsymbol{n}_{i}$ is the outer unit normal to $\partial K_{i}$. $\Lambda$ is a Lagrange multiplier introduced to impose the continuity of $\boldsymbol{v}\cdot\boldsymbol{n}$ on coarse interface.

3.2 The finite element approximation

We further partition each coarse element $K_{i}$ into a finer mesh with mesh size $h_{i}$. Let $\mathcal{T}_{h}=\cup_{i=1}^{N}\mathcal{T}_{h}(K_{i})$ be the union of all these partitions, which forms a fine partition of the domain $\Omega$. We use $h=\text{max}_{1\leq i\leq n}h_{i}$ to denote the mesh size of $\mathcal{T}_{h}$. In addition, we let $\mathcal{E}_{h}(K_{i})$ be the set of all edges (or faces in three dimensions) of the partition $\mathcal{T}_{h}(K_{i})$ and let $\mathcal{E}_{h}=\cup_{i=1}^{N}\mathcal{E}_{h}(K_{i})$ be the set of all edges (or faces) in the partition $\mathcal{T}_{h}$. Let ${\bf V}_{h,i}\times W_{h,i}\subset{\bf V}_{i}\times W_{i}$ be any of the mixed finite element spaces satisfying the inf-sup condition such as the Raviart-Thomas spaces ([24]). Define ${\bf V}_{h}=\oplus_{i=1}^{N}{\bf V}_{h,i}$ and $W_{h}=\oplus_{i=1}^{N}W_{h,i}/\mathbb{R}$ for the global discrete flux and pressure. Let $M_{H,i}\subset L_{2}(E_{i})$ be the local coarse mortar finite space, and $M_{H}=\oplus_{1\leq i\leq{|\mathcal{E}_{H}|}}M_{H,i}$ be the entire coarse mortar finite element space. Then, the multiscale mortar mixed finite element approximation of (7a)-(7c) is to seek ${\bf u}_{h}\in{\bf V_{h}}$, $p_{h}\in W_{h}$ and $\Lambda_{H}\in M_{H}$ such that for each $1\leq i\leq N$,

[TABLE]

Local conservation within each coarse element can be obtained by setting $w_{h}=1$ in Equation (8b), weak continuity of the flux is enforced across the coarse interfaces via Equation (8c). To make the above approximation be well posed, the two scales must be chosen such that the mortar space is not too rich compared to the normal traces of the velocity spaces in coarse element. For the existence of the solution to above system, we refer [7] and the reference therein.

3.3 Interface problem

The main feature of the MMMFEM is that it could be implemented by just solving a small dimension system defined on the coarse interface together with the solutions of some local problems in each coarse element.

Define bilinear forms $a_{H,i}:M_{H}\times M_{H}\rightarrow\mathbb{R},i=1,\cdots,N$ by

[TABLE]

and $a_{H}:M_{H}\times M_{H}\rightarrow\mathbb{R}$ by

[TABLE]

where \big{(}{\bf u}_{h}^{\ast}(\Lambda),p_{h}^{\ast}(\Lambda)\big{)}\in{\bf V}_{h}\times W_{h} solves ($\Lambda$ given, $q=0$)

[TABLE]

for each $1\leq i\leq N.$

Define linear functionals $g_{H,i}:M_{H}\to\mathbb{R}$ by

[TABLE]

and $g_{H}:M_{H}\to\mathbb{R}$ by

[TABLE]

where $(\boldsymbol{\bar{u}}_{h},\bar{p}_{h})\in{\bf V}_{h}\times W_{h}$ solves ($\Lambda=0,q$ given) for $1\leq i\leq N$

[TABLE]

Then, the system 8 reduce to find $\Lambda_{H}\in M_{H}$ such that

[TABLE]

It can be shown (see [6]) that the interface problem (11) produces the solution of (8a)-(8c) via

[TABLE]

where $\tilde{p}_{h}=p_{h}^{\ast}({\Lambda_{H}})+\bar{p}_{h}.$

The bilinear form $a_{H}(\cdot,\cdot)$ defined on interface is symmetric and positive semi-definite on $M_{H}$. There are several ways to implement the MMMFEM such as using a multiscale flux basis [18]. In their implementation, the matrix corresponding to the bilinear form will not be formed explicitly since they solve the linear system iteratively, after the multiscale flux basis is computed, only linear combinations are needed. Here, we use another method, that is, we assemble the matrix explicitly and apply direct solver to solve the linear system since its dimension is very small. Parallel computed can be easily adopted for the local problems. We remark that the MMMFEM shares some similarities with the Multiscale Hybrid-Mixed finite element methods [4, 20], both methods involve solving some local problems and a small dimension global problem. However, the local problems in the Multiscale Hybrid-Mixed finite element methods are Neumann boundary condition problems.

4 Basis space

The coarse space $M_{H}$ is quite important to the accuracy and efficiency of the MMMFEM. If the mortar space is small, then the resulting linear system can be solved relatively easily, however the accuracy can not be guaranteed since the continuity of the flux only be imposed very weakly. If the mortar space is too enriched, then the method will still be expensive. Various types of mortar were proposed in literature. In [7], the authors considered the polynomials for the mortar space and proved some convergence results. Several media dependent multiscale mortar space are also studied, for example, in [28], an enriched multiscale space was designed for the high-contrast flow. Homogenized multiscale basis (see [27]) was constructed for the flow in highly heterogeneous media, which can yields a moderate accurate velocity field. However this velocity field is not accurate enough to be used to transport some substance (see [9]). Residual driven multiscale online basis that includes global information was introduced in [29], successful applications for the two-phase flow simulations were presented. However, for some tough cases, the mortar space may be large ($10\sim 15$ local bases are needed in three dimensions). Our goal here is to further reduce the size of $M_{H}$ while keeping accuracy. We adopt the idea of using limited global information (see [3, 10, 14, 16, 2]) to construct the mortar space. However, only use the restriction of the global pressure field on the coarse interface to form $M_{H}$ is not enough. In particular, for each coarse interface, we will use two types functions to fill local $M_{H,i}$ at time step $t$, polynomial $P_{r}$ ($r$ is the polynomial order) and the restriction of $\Lambda_{H}$ obtained with the MMMFEM at time step $t-1$. Therefore, the space is different at different time steps. We use the MMFEM to solve the initial pressure equation. In practice, $r=0$ or $1$ is enough for simulation, therefore the local mortar space $M_{H,i}$ only consists $2\sim 4$ bases.

5 Numerical examples

In this section, we present several representative numerical experiments to show the performance of our method. We consider three highly heterogeneous permeability fields $K$ depicted in Figure 1. All these permeability fields are extracted from the tenth SPE comparative solution project (SPE10) [21], which is commonly used benchmark permeability field to test different upscaling techniques and multiscale methods. The first model is the last layer of the SPE10 dataset, as it is shown, it is highly heterogeneous and contains long channels. The second model is the first 30 layers of the SPE10 dataset, which represents a prograding near-shore environment. The third model is extracted from the last 50 layers of the SPE10 dataset which represents an Upper Ness fluvial formation, with clearly visible channels. The model 2 is smoother than the model 3. The precision of these models are $60\times 220$, $60\times 220\times 30$ and $60\times 220\times 50$ respectively. We will consider two coarsening ratio, i.e, $n=5$ and $n=10$. In the reported tables below, we use "$Nb$" to represent the number of basis functions used in the local mortar space, and "Dim" to represent the dimension of the resulting coarse linear system of the MMMFEM.

We use the quadratic relative permeability curves $k_{ro}=(1-S)^{2}$ and $k_{rw}=S^{2}$ and, along with $\mu_{w}=1$ and $\mu_{o}=5$ for the fluid viscosities. We set the initial water saturation as zero in the whole reservoir. We include 4 injectors in the corners of the reservoir and 1 producer in the center of the reservoir, we will record the oil production in this producer. The water cut (or fractional flow) is defined as the fraction of water in the produced fluid and is given by $\frac{q_{w}}{q_{t}}$, where $q_{t}=q_{o}+q_{w}$, with $q_{o}$ and $q_{w}$ being the flow rates of oil and water at the producer. We apply 10 times cheap Jacobi iterations [23] to smooth the multiscale solution. The total simulation time is 2000, we update the velocity at every 50 time steps, i.e,. the pressure equation is solved 40 times. We define the relative saturation at time step $i$ $(i=1,2\cdots,40)$ as

[TABLE]

We also denote $e_{S}$ as the average saturation error history over the time instants. For the 2D problem, we will consider two mortar spaces, i.e, multiscale basis and multiscale basis together with constant function (P0). For the 3D problems, we will consider using multiscale basis with constant function or first order polynomials (P1). Since two coarsening ratio cases will be tested, therefore, 4 multiscale mortar solutions will be computed for all test models. The reference solution is computed with the MMFEM on the fine grid.

5.1 2D test

Figure 2 shows the saturation error history of the multiscale solution with different test setting with respect to the reference solution for test model 1. We can see that if only the multiscale basis is used, the error is large especially at earlier times. if we add constant functions to the mortar space, the accuracy of the multiscale solution is improved significantly, In Figure 3, we plot the watercut comparison between different multiscale solution and the reference solution for test model 1. We also observe that if two basis functions are used, the multiscale solution can approximate the reference solution pretty well. Actually, it is hard to distinguish the line that represents the multiscale watercut when $Nb=2$ and $n=5$ and the reference watercut line. Figure 4 shows the saturation profile comparison between the multiscale solution and the reference solution. We find that when only 1 basis is used and $n=10$, the saturation profile (middle figure in 4) generated with the MMMFEM fails to capture all the details of reference saturation profile (top figure in Figure 4), however if we use 2 bases and increase the coarse-grid size, the MMMFEM solution (bottom figure in Figure 4) is almost identical as the reference solution.

5.2 3D tests

The $L^{2}$ error of the saturation is presented in Figure 5. We can observe that for each case and all time instants, the error is less than $10\%$ which clearly demonstrates the accuracy of the proposed method. Using more basis functions and increase coarse grid size will definitely reduce the error. We plot the water cut comparison for model 2 in Figure 6. The lines that represent the MMMFEM solution is very close to the reference line. All lines are almost identical at later times. Figure 7 shows the saturation profiles at $t=1000$ with different methods, we observe the results obtained using the MMMFEM look exactly the same as the fine-scale reference solution. We report the accuracy and efficiency comparisons between MMMFEM and fine-scale simulations for model 2 in Table 1. We find that for all MMMFEM solutions, the average $L^{2}$ error of the saturation is less than $5\%$. The computational time of the MMMFEM decreased by about $60\%$ to $70\%$ compared with the reference solution.

For model 3, the $L^{2}$ error and watercut comparison are shown in Figures 8 and 9 respectively. We find that if the coarsening ratio is 10, the $L^{2}$ error is $\mathcal{O}(10^{-1})$, which is not satisfiable. However, if we set $n=5$ then the accuracy of the MMMFEM solution will be improved a lot. The saturation profile comparisons at $t=1000$ are provided in Figure 10, we observe that the MMMFEM solution can capture almost all the fine-scale features. We compare the accuracy and efficiency of MMMFEM and fine-scale simulations for model 3, the results are presented in Table 2. The CPU time savings ranges from $68\%$ to $84\%$, while the average error is about $1\%\sim 11\%$. We can expect more computational cost savings for model larger than the test models here. We would like to note the relative saturation error at earlier time is larger than later times, this is because the saturation profile has sharp fronts and thus the permeability field changes a lot compared with the initial permeability field in some coarse elements. By comparing the results for model 2 and model 3, we find that the performance of the proposed method is better for smoother media, therefore, adaptive assign the basis number for different coarse interfaces based on the media’s heterogeneity may save more computational cost. We will explore this in our future research.

We also considered another type of source (source type II) for the 3D problems, that is, we inject the water in the center of the reservoir and place 4 producers in the corners of the reservoir. The saturation error history for model 2 and 3 are shown in Figures 11 and 12 respectively. Again, we find that the MMMFEM is quite accurate especially for model 3. Tables 3 and 4 show the performance comparisons of the MMMFEM and the reference solution. The average error are all less than $5\%$, similar CPU savings as previous test cases are observed.

6 Conclusions

We have developed a local-global multiscale mortar mixed finite element method for the two-phase flow problems in heterogeneous porous media. We use the multiscale mortar mixed finite element method to solve the pressure equation, and explicit finite volume method to solve the saturation equation. We propose a dynamic mortar space that consists of polynomials and multiscale basis, which is the restriction of global pressure field obtained at previous time step on the coarse interface. Fine scale solution of the pressure is obtained to initialize the simulation. We provide numerical experiments on some benchmark 2D and 3D heterogeneous models to demonstrate the performance of the proposed method.

Acknowledgment

Shubin Fu wants to thank Bingqi Yi (Sun Yat-sen University) for providing computing resources.