Splitting-based domain decomposition methods for two-phase flow with different rock types

TL;DR

This paper introduces two robust domain decomposition schemes for simulating two-phase flow in heterogeneous rocks, utilizing operator-splitting and multirate time stepping to improve efficiency and accuracy.

Contribution

The paper develops novel splitting-based domain decomposition methods with multirate time stepping for two-phase flow in different rock types, enhancing computational robustness and flexibility.

Findings

The proposed schemes are practical and accurate in numerical experiments.

Multirate time stepping improves efficiency for multi-scale problems.

Different DD methods show varying convergence behaviors.

Abstract

In this paper, we are concerned with the global pressure formulation of immiscible incompressible two-phase flow between different rock types. We develop for this problem two robust schemes based on domain decomposition (DD) methods and operator-splitting techniques. The first scheme follows a sequential procedure in which the (global) pressure, the saturation-advection and the saturation-diffusion problems are fully decoupled. In this scheme, each problem is treated individually using various DD approaches and specialized numerical methods. The coupling between the different problems is explicit and the time-marching is with no iterations. To adapt to different time scales of problem components and different rock types, the novel scheme uses a multirate time stepping strategy, by taking multiple finer time steps for saturation-advection within one coarse time step for saturation-diffusion and pressure, and permits independent time steps for the advection step in the different rocks. In the second scheme, we review the classical Implicit Pressure--Explicit Saturation (IMPES) method (by decoupling only pressure and saturation) in the context of multirate coupling schemes and nonconforming-in-time DD approaches. For the discretization, the saturation-advection problem is approximated with the explicit Euler method in time, and in space with the cell-centered finite volume method of first order of Godunov type. The saturation-diffusion problem is approximated in time with the implicit Euler method and in space with the mixed finite element method, as in the pressure problem. Finally, in a series of numerical experiments, we investigate the practicality of the proposed schemes, the accuracy-in-time of the multirate and nonconforming time strategies, and compare the convergence of various DD methods within each approach.

Tables

Table 1. Table 1: Relative tolerances.

Interf. solver	CG	OSWR	Newton	GMRes
Tol.	1E-6	1E-6	1E-4	1E-3
Subd. solver	Newton		GMRes
Tol	1E-4		1E-4

Table 2. Table 2: Test case 1: computational iterations of the Newton-GMRes method.

Newton-GMRes method
Interf. Newton	Interf. GMRes	Subd. Newton
without precond.
Tot. Avg.	Tot. Avg.	Tot. Avg.
318  6.36	835  16.72	1660 13.23
with precond.
Tot. Avg.	Tot. Avg.	Tot. Avg.
206  4.11	705  11.2	1510 9.21

Table 3. Table 3: Test case 1: computational iterations of the OSWR method.

OSWR method
Interf. OSWR	Subd. Newton
Tot. Avg.	Tot. Avg.
403  8.15	1695 13.19

Table 4. Table 4: Test case 3 (multiple rocks): the physical properties.

	$K$ (md)	$\varphi$	$\alpha$
Rock 1	6	0.5	0.3
Rock 2	3	0.5	0.3
Rock 3	0.6	0.3	5
Rock 4	0.3	0.3	5

Equations

$$\mathrm{\Phi }\frac{\partial s_{\ell }}{\partial t}+\nabla \cdot \left(-\mathbf{K}{k}_{\ell }(s_{\ell })(\nabla p_{\ell }-{\rho }_{\ell }{\mathbf{u}}_g)\right)=0,$$

$$s_w+s_n=1,$$

$$p_n-p_w=\pi (s_w),$$

$$s_w(\cdot ,0)=s_0,\text{ in }\mathrm{\Omega }.$$

$$\mathbf{K}{k}_{\ell }(s_{\ell })(\nabla p_{\ell }-{\rho }_{\ell }{\mathbf{u}}_g)\cdot \mathbf{n}=0,\text{ on }\partial \mathrm{\Omega }\times (0,T).$$

$$p:=\frac{1}{2}(p_w+p_n)+\beta (s),$$

$$\alpha (s):={\int }_0^s-\frac{k_w(a)k_n(a)}{k_w(a)+k_n(a)}{\pi }^{\prime }(a)\text{d}a,\forall s\in (0,1).$$

$$\mathrm{\Phi }\frac{\partial s}{\partial t}+\nabla \cdot {\mathbf{u}}_w=0,$$

$${\mathbf{u}}_w=f(s)\mathbf{u}+f_g(s){\mathbf{u}}_g-\mathbf{K}\nabla \alpha (s),$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}=-\mathbf{M}(s)\left(\nabla p-\rho (s){\mathbf{u}}_g\right),$$

$${\mathbf{u}}_w\cdot \mathbf{n}=0,\text{and}{\mathbf{u}}_i\cdot \mathbf{n}=0,$$

$$s(\cdot ,0)=s^0,$$

$$f_g(s):=\mathbf{K}\frac{k_w(s)k_n(s)}{k_w(s)+k_n(s)}({\rho }_w-{\rho }_n),\text{and }\rho (s):=\frac{k_w(s){\rho }_w+k_n(s){\rho }_n}{k_w(s)+k_n(s)}.$$

$${\mathrm{\Phi }}_i\frac{\partial s_i}{\partial t}+\nabla \cdot {\mathbf{u}}_{wi}=0,$$

$${\mathbf{u}}_{wi}=f_i(s_i){\mathbf{u}}_i+f_{gi}(s_i){\mathbf{u}}_g-{\mathbf{K}}_i\nabla {\alpha }_i(s_i),$$

$$\nabla \cdot {\mathbf{u}}_i=0,$$

$${\mathbf{u}}_i=-{\mathbf{M}}_i(s_i)\left(\nabla p_i-{\rho }_i(s_i){\mathbf{u}}_g\right),$$

$${\mathbf{u}}_{wi}\cdot {\mathbf{n}}_i=0,\text{and}{\mathbf{u}}_i\cdot {\mathbf{n}}_i=0,$$

$$s_i(\cdot ,0)=s^0,$$

$$[[\pi (s)]]=0,$$

$$[[{\mathbf{u}}_w\cdot \mathbf{n}]]=0,$$

$$[[v]]:=v_1-v_2,\text{ and }[[\mathbf{v}\cdot \mathbf{n}]]:={\mathbf{v}}_1\cdot {\mathbf{n}}_1+{\mathbf{v}}_2\cdot {\mathbf{n}}_2.$$

$$[[p-\beta (s)]]=0,$$

$$[[\mathbf{u}\cdot \mathbf{n}]]=0,$$

$${\mathbf{W}}_{h,i}:=\{{\mathbf{v}}_h\in {\mathbf{RTN}}_0({\mathcal{T}}_{h,i})\cap \mathbf{H}(\text{div},{\mathrm{\Omega }}_i);{\mathbf{v}}_h\cdot \mathbf{n}=0\text{ on }{\mathrm{\Gamma }}_i^{\text{N}}\},$$

$$M_{h,i}:={\mathbb{P}}_0({\mathcal{T}}_h).$$

$${\mathrm{\Lambda }}_h:=\{{\mu }_h\in L^2(\mathrm{\Gamma });{\mu }_h{|}_{\sigma }\in {\mathbb{P}}_0(\sigma ),\forall \sigma \in {\mathcal{E}}_h^{\mathrm{\Gamma }}\}.$$

$${\left({\mathcal{P}}_{ij}\varphi \right)}_{|I_i^{n,l}}:=\frac{1}{{\tau }_i^{n,l}}\sum _{m=1}^{{\mathcal{N}}_j^{\text{a}}}{\int }_{I_i^{n,l}\cap I_j^{n,m}}\varphi \text{d}t,\forall l\in \{1,\dots ,{\mathcal{N}}_i^{\text{a}}\}.$$

$$-(p_{h,i}^n,\nabla \cdot \mathbf{v}{)}_{{\mathrm{\Omega }}_i}+({\mathcal{A}}_i^n{\mathbf{u}}_{h,i}^n,\mathbf{v}{)}_{{\mathrm{\Omega }}_i}$$

$$=({\rho }_i(s_{h,i}^n){\mathbf{u}}_g,\mathbf{v}{)}_{{\mathrm{\Omega }}_i},$$

$$(\nabla \cdot {\mathbf{u}}_{h,i}^n,q{)}_{{\mathrm{\Omega }}_i}$$

$$\displaystyle\begin{array}[]{lll}\displaystyle\int_{\Gamma}[\![{\bf u}^{n}_{h}\cdot{\bf n}]\!]\mu=0,&\quad\forall\mu\in\mathrm{\Lambda}_{h}.\end{array}$$

$$p_{h,i}^n:={\overline{p}}_{h,i}^n(\lambda )+{\mathring{p}}_{h,i}^n\text{and}{\mathbf{u}}_{h,i}^n:={\overline{\mathbf{u}}}_{h,i}^n(\lambda )+{\mathring{\mathbf{u}}}_{h,i}^n,$$

$$-({\overline{p}}_{h,i}^n,\nabla \cdot \mathbf{v}{)}_{{\mathrm{\Omega }}_i}+({\mathcal{A}}_i^n{\overline{\mathbf{u}}}_{h,i}^n,\mathbf{v}{)}_{{\mathrm{\Omega }}_i}$$

$$(\nabla \cdot {\overline{\mathbf{u}}}_{h,i}^n,q{)}_{{\mathrm{\Omega }}_i}$$

$$-{\left({\mathring{p}}_{h,i}^n,\nabla \cdot \mathbf{v}\right)}_{{\mathrm{\Omega }}_i}+{\left({\mathcal{A}}_i^n{\mathring{\mathbf{u}}}_{h,i}^n,\mathbf{v}\right)}_{{\mathrm{\Omega }}_i}$$

$$={\left({\rho }_i(s_{h,i}^n){\mathbf{u}}_g,\mathbf{v}\right)}_{{\mathrm{\Omega }}_i},$$

$${\left(\nabla \cdot {\mathring{\mathbf{u}}}_{h,i}^n,q\right)}_{{\mathrm{\Omega }}_i}$$

$$s_{\Gamma ,i}^n(\lambda ,\mu )$$

Full text

Key words: Two-phase flow in porous media; continuous capillary pressure; non-overlapping domain decomposition; operator-splitting; mixed finite element method; finite volume method; multirate and non-conforming time grids.

1 Introduction

Numerical simulation of two-phase flow in porous formations has been the subject of investigation of many researchers owing to important applications in e.g., the management of petroleum reservoirs or environmental remediation. Two-phase flow modeling has received increased attention in connection with porous medium with different rock types, so that the permeability and the capillary pressure field are changing across the interfaces between the rocks.

Two-phase flow in porous media can be modeled by mass balance laws for each of the fluids [21, 23]. In particular, an equivalent formulation can be obtained for the system of equations governing two-phase immiscible, incompressible flow in porous media, by introducing an artificial variable called the global pressure. The dependent variables in such formulation are the global pressure and the wetting phase saturation [21]. Considering these variables, the governing equations consist of a nonlinear elliptic Darcy’s law for the global pressure and a nonlinear degenerate parabolic equation of advection-diffusion type for the wetting phase saturation. Flow simulation of such systems in porous media with different rock types is very difficult because of variations in physical parameters in the different rock types that require coupling [9, 22]. Particularly, due to changes in permeability and capillary force, the global pressure and saturation exhibit strong discontinuities across the interfaces between the rocks (cf. [27, 50]). Several numerical schemes have been developed for two-phase flow in heterogeneous media. Finite volume schemes have been proposed in [18, 26] for a simplified two-phase model involving only diffusion effects. One can see [19] for the extension to the full two-phase flow between different rock types. A discontinuous Galerkin (DG) method or a combination of DG method and mixed finite elements (MFE) method was employed in [13, 27, 39]. For domain decomposition methods applied to two-phase flow problems, one can cite in particular the work in [33], where a fully implicit-multiscale mortar mixed finite element method for two-phase flow in a heterogeneous porous medium is presented. In [7], a space–time domain decomposition method formulated using Robin and Ventcell type coupling conditions was applied to a simple two-phase flow model between different rock types (cf. [2, 6]). We also refer the reader to linear domain decomposition methods for two-phase flow problems in homogeneous porous media presented in [57, 58]. Therein, the problem of two-phase flow is decomposed into a set of subproblems in different subdomains, and solved at each time step with the $L$-scheme linearization method; a robust quasi-Newton method with a parameter $L>0$ mimicking the Jacobian from the Newton method [46, 49, 52, 53]. Other related works can be found in [3, 4, 35, 37, 51, 54].

Our contribution in this paper is to formulate robust domain decomposition methods for the two-phase flow between different rock types. The methods are motivated by the fact that the (global) pressure and saturation as well as advection and diffusion effects act on different time scales within the same rock type and also between different rock types, and it is easy to implement this in the domain decomposition context. However, classical domain decomposition methods are known to perform poorly if diffusion is dominated by advection as well as in purely advective problems (cf. [20] for more details). From these observations, we investigate new procedures based on a divide-and-conquer (operator-splitting) strategy for the numerical solution of the two-phase flow between different rock types [11, 16, 38, 41]. In the first approach, we split the global pressure, the saturation-advection and the saturation-diffusion problems and treat them individually using specialized numerical methods (see Scheme 1 for the workflow strategy). This workflow will then apply various domain decomposition approaches adapted to each problem, transforming the original problem into a sequential solves of reduced problems of pressure, saturation-diffusion and saturation-diffusion posed only on the physical interfaces between rock types, through the use of trace operators. The DD approaches are based either on the Steklov-Poincaré type methods [24, 35, 47] or the optimized Schwarz methods [30, 31, 37].

The resulting interface problems of the first approach can be solved by various iterative methods. For the discretization in the subdomains, the saturation-diffusion is approximated with the implicit Euler method in time, and in space with the hybrid mixed finite element method, as in the pressure problem. To handle efficiently the advection step, a numerical scheme is presented that uses the hybrid finite volume method based on the Godunov scheme (cf. [1, 12, 48]), with the explicit Euler method in time. The coupling strategy between the three processes is explicit and the time marching is without any iterations. It is known that the time step must be sufficiently small for the stability of such explicit-coupling schemes [25]. To overcome this difficulty and reveal the multi-physical process of the problem, the algorithm allows a multirate time stepping strategy, in which adequate time steps are used for each of the problems (see Scheme 1). Here, we consider applications, where advection effects are often much greater than capillary diffusion effects. To reduce the computational cost (i) we solve the (global) pressure and saturation-diffusion problems in the coarse time step and solve the advection problem in the finer time step (ii) we employ nonconforming time grids between different rock types for the advection step. For an overview and further insight on multirate schemes, we refer the reader to [11, 34, 56] where various coupled problems are considered (see also the adaptive strategy in [8]). For local/nonconforming time discretizations in the context of DD approaches, we refer to [35, 37] and the references therein.

In the second approach of this paper, we review the classical Implicit Pressure–Explicit Saturation (IMPES) method (cf. [42]) in the context of DD methods. This approach consists on decoupling only (global) pressure and saturation problems in order to solve sequentially each problem for its main variable. The pressure problem is discretized as in the first approach with MFE method, and the coupled advection-diffusion problem for the saturation, is discretized with a scheme explicit in time and with a cell-centered finite volume method in space. In this approach, we employ larger time steps for the pressure problem with respect to the saturation problem, and to further save computational effort, we use different time-steps for the saturation in the different rock types.

The two approaches of this paper are iteration-free and allow easy reuse of existing codes specialized to each component of the problem, which minimizes the computational effort. Furthermore, the approaches can easily integrate more complex problems; the first approach is extended to address reduced fracture model between different rock types presented in [5].

The remainder of this paper is organized as follows. In Section 2, we recall briefly the global pressure formulation for incompressible two-phase flow problem. We then rewrite this problem for the case of flow between different rock types and describe the physical conditions occurring at the interface between the rocks. In Section 3, we present the numerical approaches. We show numerical results in Section 4. Section 5 contains concluding remarks. In Appendix A, we present a brief application of the first approach on a reduced fracture model for two-phase flow between different rock types.

2 Two-phase Darcy

flow model between different rock types

Let $\Omega$ be an open bounded domain of ${\amsmathbb R}^{d}$ ($d=2$ or $3$) which is assumed to be polygonal if $d=2$ and polyhedral if $d=3$. We denote by $\partial\Omega$ its boundary (supposed to be Lipschitz-continuous) and by ${\bf n}$ the unit normal to $\partial\Omega$, outward to $\Omega$.

2.1 Global pressure formulation for incompressible two-phase flow

We consider the immiscible incompressible two-phase flow in a porous medium. Assuming that there are only two phases occupying the porous medium $\mathrm{\Omega}$, and that each phase is composed of a single component, the mathematical form of this problem is as follows [22]:

[TABLE]

where the unknowns are $s_{\ell}$, the phase saturation, and $p_{\ell}$, the phase pressure, $\ell\in\{w,n\}$. The subscripts $w$, $n$ stand for wetting and non-wetting, respectively. Typically, the non-wetting phase is oil and the wetting one is water. The function $\mathrm{\Phi}=\mathrm{\Phi}(\mathbf{x})$ denotes the porosity of the rock ($\mathrm{\Phi}\in(0,1)$ in the domain $\mathrm{\Omega}$) and ${\bf u}_{g}$ denotes the gravity field; the permeability of the porous medium $\mathbf{K}$ is assumed to be a positive scalar function, the mobility $k_{\ell}$ of the phase $\ell$ is an increasing function of the saturation $s_{\ell}$, satisfying $k_{\ell}(0)=0$ and $k_{\ell}(1)=1$, and $\rho_{\ell}$ denotes the density of phase $\ell$, $\ell\in\{w,n\}$. The function $\pi$ is the capillary pressure, defined as the difference between the phase pressures and assumed to be a function (positive and decreasing) of the saturation $s_{w}$. More details on capillary pressure and relative permeability can be found in [22, 59]. We assume that the initial phase saturation is known, say

[TABLE]

We also assume homogeneous Neumann boundary conditions on the phase fluxes, i.e.,

[TABLE]

In order to avoid some of the known difficulties related to the degeneracy of the problem (2.1)–(2.3), we follow the classical idea of introducing the so-called global pressure formulation [21]. For the simplicity of notation let $s:=s_{w}$, so that $s_{n}:=1-s$. Define the total mobility function $\mathbf{M}(s):=\mathbf{K}\left(k_{w}(s)+k_{n}(s)\right)$ and introduce the fractional flow function $f(s):={\displaystyle\frac{k_{w}(s)}{k_{w}(s)+k_{n}(s)}}$ and let $\beta(s):=\displaystyle{\int_{1}^{s}}(\frac{1}{2}-f(a))\pi^{\prime}(a)\textnormal{d}a,\,\forall s\in(0,1)$. We then introduce the global pressure function $p$ by

[TABLE]

and the Kirchoff transform

[TABLE]

Following [21], the problem (2.1)–(2.3) can be rewritten under the form

[TABLE]

where

[TABLE]

In view of the above formulation, in equation (2.5d) the saturation appears only through the coefficients $\mathbf{M}(s)$ and $\rho(s)$, and not through its gradient $\nabla s$. Hence, the equation (2.5d) looks very much like Darcy’s law for a fictitious fluid representing the global (water + oil) flow pattern. The global pressure is not a physical pressure but a good numerical and methematical unknown [21, 22]. Precisely, it is a smooth function defined in the whole domain, whether a phase vanishes or not, while this is not the case for the phase pressures. Still following [21], the global pressure satisfies $p_{w}\leq p\leq p_{n}$, and such that $p=p_{w}$ if $s_{w}=1$ and $p=p_{n}$ if $s_{w}=0$. Also, by plugging equation (2.5b) into (2.5a), the resulting saturation equation gives a nonlinear parabolic equation of advection-diffusion type given by the sum of a diffusion contribution due to capillary effects and an advection contribution due to gravity and total flow effects.

2.2 Flow between different rock types

In this part, we particularize the model problem (2.5) to a porous medium with different rock types, following [10]. Precisely, we suppose that $\mathrm{\Omega}$ is divided into two non-overlapping polyhedral subdomains $\mathrm{\Omega}_{i}$, $i\in\{1,2\}$. The exterior boundary of $\mathrm{\Omega}_{i}$ is denoted by $\mathrm{\Gamma}_{i}^{\textnormal{N}}:=\partial\mathrm{\Omega}_{i}\cap\partial\mathrm{\Omega}$, $i\in\{1,2\}$. Let $\mathrm{\Gamma}$ be the interface between the subdomains, i.e., $\mathrm{\Gamma}:=\partial\mathrm{\Omega}_{1}\cap\partial\mathrm{\Omega}_{2}$ (see Figure 1 (left)). Let ${\bf n}_{i}$ denotes the unit outward normal on $\partial\mathrm{\Omega}_{i}$. In the following, we consider that each subdomain $\mathrm{\Omega}_{i}$ corresponds to a (homogeneous) rock type. This means that across $\Gamma$, not only the porosity and the absolute permeability differ, but also the relative permeability and capillary pressure curves. For every function in (2.5) involving a space-dependent quantity (e.g., $f(s)$) that is valid separately in each of the two subdomains (see Figure 1 (left)), we will use subscripts (e.g., $f_{i}(s_{i})$).

Now that the decomposition of $\Omega$ is made, the two-phase flow equations (2.5) remain valid in each subdomain $\mathrm{\Omega}_{i}$, $i\in\{1,2\}$, i.e.,

[TABLE]

where the unknown functions for the saturation equations (2.7a)-(2.7b) are $s_{i}$, the subdomain saturation, and ${\bf u}_{wi}$, the subdomain velocity for the wetting phase, $i\in\{1,2\}$. That of the global pressure equations (2.7c)-(2.7d) are $p_{i}$, the subdomain (global) pressure, and ${\bf u}_{i}$, the subdomain total velocity, $i\in\{1,2\}$.

Now we come to the transmission conditions across the interface $\mathrm{\Gamma}$. The subdomain saturation problems (2.7a)-(2.7b), for $i\in\{1,2\},$ are coupled through the following interface conditions [10]:

[TABLE]

where the jumps on $\mathrm{\Gamma}$ are defined as

[TABLE]

The condition on the capillary pressure (2.8a) stems from the continuity of phase pressures since $\pi(s)=p_{n}-p_{w}$, resulting in a jump on the saturation as depicted in Figure 1 (right). The condition (2.8b) imposes the continuity of the normal component of the wetting flow ${\bf u}_{w}$, to preserve phase conservation. The subdomain pressure problems (2.7c)-(2.7d), for $i\in\{1,2\},$ are coupled through the following interface conditions [5]:

[TABLE]

Equation (2.9a) justifies the continuity of the phase pressures since $p-\beta(s)=\frac{1}{2}(p_{w}+p_{n})$. Equation (2.9b) represents conservation of the total fluid.

Remark 1** (Compatibility condition)**

We present the solution procedures for the case of continuous capillary pressure fields. That case can be realized under the compatibility condition that the capillary pressure curves have the same endpoints as in Figure 1 (right). The approaches can be applied to discontinuous capillary pressure by extending the curves of these functions into monotone graphs(cf. [12, 19]).

3 Non-overlapping DD methods with explicit coupling

We present in this section two numerical approaches for the two-phase flow problem (2.7)–(2.9) based on domain decomposition methods and operator splitting techniques.

3.1 A brief description of operator-splitting methods

The natural IMPES-based approach for solving (2.7)–(2.9) is to use a sequential splitting in which one fixes the saturation in the pressure equations (2.7c)-(2.7d) with the corresponding coupling conditions (2.9) and solves for the unknown pressure and velocity, which are subsequently held fixed when evolving the saturation a time step according to (2.7a)-(2.7b) with the coupling conditions (2.8). Here, we involve a further operator splitting of this problem and split the advection-diffusion operator (see Remark 2), and apply a multirate time stepping strategy exploiting the different time scales of the problem components (details on the flowchart of this method and its modules are given in Scheme 1). In grosso modo, once the pressure system is solved through the module PRESSURE_SOLVER at the coarse time step, say $t^{n}$, with time-lagging the saturation, we advance this later through the module ADVECTION_SOLVER for multiple finer time steps, say $t^{n,l}$, $1\leq l\leq\mathcal{N}^{\textnormal{a}}$ (with $t^{n,\mathcal{N}^{\textnormal{a}}}=t^{n+1}$), and finally advance saturation through the module DIFFUSION_SOLVER at the coarse time step $t^{n+1}$. Such approach allows for the use of various schemes and DD techniques taking the advantage of the specific properties of each problem, permits the use of different time steps for advection effects to adapt the different time scales in the different rocks, and generally reduce memory requirements.

Remark 2** (Advection-dominated effects)**

Recalling that plugging (2.7b) into (2.7a) results on a nonlinear advection-diffusion problem for the saturation. For many applications, advection effects are much greater than capillary-diffusion effects through porous formations. Therefore, it is efficient to compute fast (advection-effects) solutions on a fine mesh in time and slow (total velocity + capillarity-effects) solutions on a coarse mesh in time.

3.2 Meshes and function spaces

We introduce here the partitions of $\mathrm{\Omega}$, time discretization, notation, and function spaces; see [5, 7, 26] for the standard part of the notation.

3.2.1 Temporal meshes

For integer values $N\geq 0$, let $\left(\tau^{n}\right)_{1\leq n\leq N}$ denote a sequence of positive real numbers corresponding to the discrete time steps such that $T:=\sum_{n=1}^{N}\tau^{n}$. Let $t^{0}:=0$, and $t^{n}:=\sum_{j=1}^{n}\tau^{j},\ 1\leq n\leq N$ be the discrete times. Let $I^{n}:=(t^{n-1},t^{n}],\ 1\leq n\leq N$. For every time step $0\leq n\leq N$, we let $v_{h}^{n}:=v_{h}(\cdot,t^{n})$ for any sufficiently smooth function $v_{h}$.

3.2.2 Spatial meshes

We consider a spatial discretization $\mathcal{T}_{h,i}$ of the domain $\mathrm{\Omega}_{i}$ consisting of either simplicial or rectangular elements $K\in\mathcal{T}_{h,i}$. We assume that these meshes are such that $\mathcal{T}_{h}:=\displaystyle{\cup_{i=1}^{2}\mathcal{T}_{h,i}}$ forms a conforming finite element mesh on all of $\Omega$ (see Figure 2);

we assume that this partition is matching in the sense that for two distinct elements of $\mathcal{T}_{h}$ their intersection is either an empty set or their common vertex or edge. We denote by $h$ the maximal mesh size of $\mathcal{T}_{h,i}$. The interior mesh faces in $\mathcal{T}_{h,i}$ are collected into the set $\mathcal{E}_{h,i}^{\textnormal{int}}$, and we denote by $\mathcal{E}_{h,i}$ all the faces of $\mathcal{T}_{h,i}$ and we set $\mathcal{E}_{h}:=\displaystyle\cup_{i=1}^{2}\mathcal{E}_{h,i}$. We denote by $\mathcal{E}_{h,i}^{\textnormal{N}}$ the sides of $\mathcal{E}_{h}$ on $\mathrm{\Gamma}^{\textnormal{N}}_{i}$. Finally, let $\mathcal{E}^{\mathrm{\Gamma}}_{h}$ be a partition of $\mathrm{\Gamma}$ given by the sides of $\mathcal{T}_{h}$ on $\mathrm{\Gamma}$ and we denote by $\mathcal{E}_{K}$ the faces of the element $K\in\mathcal{T}_{h}$. The volume of an element $K$ is denoted by $|K|$ and that of a face $\sigma$ by $|\sigma|$. Finally, we use the notation $\mathbf{x}_{K}$ to denote the “center” of the cell $K\in\mathcal{T}_{h}$. If $\sigma=K|L\in\mathcal{E}_{h}$ separates the cells $K$ and $L$, $d_{K,L}$ denotes the Euclidean distance between $\mathbf{x}_{K}$ and $\mathbf{x}_{L}$, and $d_{K,\sigma}$ for $\sigma\in\mathcal{E}_{K}$ denotes the distance from $\mathbf{x}_{K}$ to $\sigma$. We assume that the composite mesh $\mathcal{T}_{h}$ satisfies the following orthogonality condition: for a face $\sigma=K|L\in\mathcal{E}_{h}$, the line segment $\mathbf{x}_{K}\mathbf{x}_{L}$ is orthogonal to $\sigma$ (see [54]).

3.2.3 Function spaces

We denote by $\mathbb{P}_{l}(S)$ the space of polynomials on a subdomain $S\subset\mathrm{\Omega}$ of total degree less than or equal to $l$, and by $\mathbb{P}_{l}(\mathcal{T}_{h}):=\{p_{h}\in L^{2}(\Omega);\,p_{h}|_{K}\in\mathbb{P}_{l}(K),\ \forall K\in\mathcal{T}_{h}\}$. Let $\mathbf{H}(\textnormal{div};\mathrm{\Omega})\;$ be the space of vector-valued functions from $\left[L^{2}(\mathrm{\Omega})\right]^{d}$ that admit a weak divergence in $L^{2}(\mathrm{\Omega})$. We additionally consider the piecewise lowest-order Raviart–Thomas–Nédélec space $\mathbf{RTN}_{0}(\mathcal{T}_{h})\subset L^{2}(\Omega)$ defined by $\mathbf{RTN}_{0}(\mathcal{T}_{h}):=\{\mathbf{v}_{h}\in L^{2}(\Omega);\,\mathbf{v}_{h}|_{K}\in\mathbf{RTN}_{0}(K),\;\forall K\in\mathcal{T}_{h}\}$, where $\mathbf{RTN}_{0}(K):=\left[\mathbb{P}_{0}(\mathcal{T}_{h})\right]^{d}+\mathbb{P}_{0}(\mathcal{T}_{h})\mathbf{x}$ (see Figure 3).

For the discretization of the subdomain variables, we introduce the discrete spaces:

[TABLE]

The DD algorithms utilize Lagrange multipliers on the interfaces between the rocks to impose weakly interface conditions. Thus, to impose the matching conditions (2.8)-(2.9), we introduce the mortar finite element space:

[TABLE]

3.2.4 Multirate and nonconforming time grids

For $1\leq n\leq N$, let $\mathcal{D}_{1}^{n}$ and $\mathcal{D}_{2}^{n}$ be two possibly different partitions of the time interval $I^{n}$. The finer time step for the advection problem in each interval $I^{n}$ is defined by a sequence $\left(\tau^{n,l}_{i}\right)_{1\leq l\leq\mathcal{N}^{\textnormal{a}}_{i}}$, where $\mathcal{N}^{\textnormal{a}}_{i}$ is a positive integer value such that $\tau^{n}:=\sum_{i=1}^{\mathcal{N}^{\textnormal{a}}_{i}}\tau^{n,l}_{i}$. We let $t^{n,0}:=t^{n}$, $t^{n,\mathcal{N}^{\textnormal{a}}_{i}}_{i}:=t^{n+1}$, and $t^{n,l}_{i}:=\sum_{j=1}^{l}\tau^{n,l}_{j}$, $1\leq l\leq\mathcal{N}^{\textnormal{a}}_{i}$, be the discrete times for the advection in $I^{n}$. We let $I^{n,l}_{i}=(t_{i}^{n,l-1},t^{n,l}_{i}],\,1\leq l\leq\mathcal{N}^{\textnormal{a}}_{i}$ (see Figure 4).

We denote by $\mathbb{P}_{0}(\mathcal{D}_{j}^{n},\mathrm{\Lambda}_{h})$ the space of piecewise constant functions in time on grid $\mathcal{D}_{j}^{n}$ with values in $\mathrm{\Lambda}_{h}$. To perform the exchange of the data between the fine time grids $\mathcal{D}_{1}^{n}$ and $\mathcal{D}_{2}^{n}$, $1\leq n\leq N$, (see Figure 4 (right)), we define the $L^{2}$-projection $\mathcal{P}_{ij}$ from $\mathbb{P}_{0}(\mathcal{D}_{j}^{n},\mathrm{\Lambda}_{h})$ to $\mathbb{P}_{0}(\mathcal{D}_{i}^{n},\mathrm{\Lambda}_{h})$, as follows [35, 6]: for all $0\leq n\leq N$, for $\phi\in\mathbb{P}_{0}(\mathcal{D}_{j}^{n},\mathrm{\Lambda}_{h})$,

[TABLE]

We refer the reader to [32] for details on the implementation of this projection using an optimal algorithm with no additional grid.

3.3 The module PRESSURE_SOLVER

Let $0\leq n\leq N$ be fixed. The module PRESSURE_SOLVER takes as input the saturation at the present time to linearize the nonlinear coefficients. So, for known saturations $s^{n}_{h,i}$, $i\in\{1,2\}$, we let $\mathcal{A}_{i}^{n}=\mathbf{M}_{i}(s^{n}_{h,i})^{-1}$, and introduce a new unknown $\lambda:=(p^{n}_{h,i}-\beta_{i}(s^{n}_{h,i}))|_{\mathrm{\Gamma}}$, such that $({\bf u}^{n}_{h,i},p^{n}_{h,i},\lambda)\in\mathbf{W}_{h,i}\times M_{h,i}\times\mathrm{\Lambda}_{h}$ solves

[TABLE]

This system is closed by the weak interface condition

[TABLE]

To efficiently solve the above multi-domain pressure problem, this system will be reduced to an interface problem posed only on $\Gamma$ [3]. The reduced problem is then solved by an iterative procedure, which requires solving subdomain pressure problems at each iteration. For $i\in\{1,2\}$, we let

[TABLE]

where for $\lambda\in\mathrm{\Lambda}_{h}$, $(\overline{{\bf u}}^{n}_{h,i},\overline{p}^{n}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$ solves

[TABLE]

and where $(\mathring{{\bf u}}^{n}_{h,i},\mathring{p}^{n}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$ solves

[TABLE]

Define the forms $s^{n}_{\Gamma,i}:\mathrm{\Lambda}_{h}\times\mathrm{\Lambda}_{h}\rightarrow{\amsmathbb R}$, $i\in\{1,2\}$, $s^{n}_{\Gamma}:\mathrm{\Lambda}_{h}\times\mathrm{\Lambda}_{h}\rightarrow{\amsmathbb R}$, and $g^{n}_{\Gamma}:\mathrm{\Lambda}_{h}\rightarrow{\amsmathbb R}$,

[TABLE]

where $\mathcal{S}^{\textnormal{DtN}}_{\Gamma,i,n}:\mathrm{\Lambda}_{h}\rightarrow\mathrm{\Lambda}_{h}$, $1\leq i\leq 2$, and $\mathcal{S}^{\textnormal{DtN}}_{\Gamma,n}:=\sum_{i=1}^{2}\mathcal{S}^{\textnormal{DtN}}_{\Gamma,i,n}$ are Dirichlet-to-Neumann or Steklov-Poincaré (SP) type operators. Obviously, the operator $\mathcal{S}^{\textnormal{DtN}}_{\Gamma,i,n}$ is linear. It is now easy to verify that the pressure problem (3.2)–(3.4) is reduced to an interface problem:

Definition 3** (Steklov-Poincaré problem)**

Find $\lambda\in\mathrm{\Lambda}_{h}$ such that

[TABLE]

One can show that the Dirichlet-to Neumann operator $\mathcal{S}^{\textnormal{DtN}}_{\Gamma,n}$ associated to (3.9) is a symmetric positive definite one so we can use a conjugate gradient method (CG) to calculate $\lambda$. The subdomain solutions are then retrieved via (3.5).

Remark 4** (Equivalence with the multidomain pressure problem)**

In Definition 3, the continuity condition $[\![p^{n}_{h}-\beta(s^{n}_{h})]\!]=0$ on $\Gamma$ is imposed in space by fixing $\lambda:=(p^{n}_{h,i}-\beta_{i}(s^{n}_{h,i}))|_{\mathrm{\Gamma}}$, $\forall i\in\{1,2\}$, while the continuity of the total flux $[\![{\bf u}^{n}_{h}\cdot{\bf n}]\!]=0$ is imposed weakly via (3.9).

In order to speed up the convergence of the CG iterations, we may use a Neumann–Neumann preconditioner with weighted matrices (cf. [35]): for the parameters $a_{i}\in{\mathbb{R}}$ such that $a_{1}+a_{2}=1$, we introduce the Neumann-to-Dirichlet operators $\mathcal{S}^{\textnormal{NtD}}_{\Gamma,i,n}$, $i\in\{1,2\}$, given by

[TABLE]

where $({\bf u}^{n}_{h,i}(\mathrm{\varphi}),p^{n}_{h,i}(\mathrm{\varphi}))$ is such that ${\bf u}^{n}_{h,i}\cdot\mathbf{n}_{i}=\mathrm{\varphi}$ on $\mathrm{\Gamma}$ and solves

[TABLE]

The preconditioned version (in strong form) of (3.9) is given by

Definition 5** (Preconditioned problem)**

Find $\lambda\in\mathrm{\Lambda}_{h}$

[TABLE]

Following [35], an alternative to the Steklov-Poincaré formulation of Definition 3, can be obtained through an equivalent formulation of the pressure model problem (2.7c)-(2.7d) with the interface conditions (2.9). Precisely, one can solve equations (2.7c)-(2.7d), for $i\in\{1,2\}$, together with Robin transmission conditions

[TABLE]

where $\mathrm{\gamma}_{1}$ and $\mathrm{\gamma}_{2}$ are strictly positive constants. Note that (3.12) are proper linear combinations of the original conditions within PRESSURE_SOLVER module given by $[\![{\bf u}^{n}_{h}\cdot{\bf n}]\!]=0$ and $[\![p^{n}_{h}-\beta(s^{n}_{h})]\!]=0$ on $\Gamma$. Therefore, adopting ideas from [2, 35, 37, 58], one can use an *Optimized Schwarz Waveform Relaxation *(OSWR) iterative method for solving this problem. This method can be written as follows:

Algorithm 1** (The linear OSWR)**

Choose an initial approximation $\mathrm{\eta}^{-1}_{i}\in\mathrm{\Lambda}_{h}$ of $\mathrm{\eta}_{i}:={\bf u}^{n}_{h,j}\cdot{\bf n}_{j}+\gamma_{i}(p^{n}_{h,j}-\beta(s^{n}_{h,j}))$, for $j=(3-i)$, and a tolerance $\epsilon>0$. Set $k:=-1$. 2. 2.

Do**

1. (a)

Increase $k:=k+1$. 2. (b)

Compute $({\bf u}^{n}_{h,i},p^{n}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$, $i\in\{1,2\}$, such that

[TABLE] 3. (c)

Set $\mathrm{\eta}^{k}_{i}:={\bf u}^{k,n}_{h,j}\cdot{\bf n}_{j}+\gamma_{i}(p^{k,n}_{h,j}-\beta(s^{n}_{h,j}))$, for $j=(3-i)$.

While* $\dfrac{\|(\mathrm{\eta}^{k}_{1},\mathrm{\eta}^{k}_{2})-(\mathrm{\eta}^{k-1}_{1},\mathrm{\eta}^{k-1}_{2})\|_{L^{2}(\Gamma)}}{\|(\mathrm{\eta}^{k-1}_{1},\mathrm{\eta}^{k-1}_{2})\|_{L^{2}(\Gamma)}}\geq\epsilon$.*

The OSWR method in the context of MFE methods has been studied in [37]. Note that the parameters $(\mathrm{\gamma}_{1},\mathrm{\gamma}_{2})$ can be optimized to give an improved convergence rate of this algorithm. This can be done by numerically minimizing the convergence factor corresponding to the Darcy problem (cf. [2, 37] for more details).

3.4 The module SPLIT_SATURATION_SOLVER

Here, we detail the ingredients of the module SPLIT_SATURATION_SOLVER. As pointed previously, the saturation-advection within ADVECTION_SOLVER module takes multiple finer time steps against one saturation-diffusion time step within DIFFUSION_SOLVER. The module ADVECTION_SOLVER employs different time grids between subdomains to adapt the different time scales in the different rocks [37].

3.4.1 STEP 1: ADVECTION_SOLVER

The scheme in ADVECTION_SOLVER, is based on conservation of mass element-by-element; the subdomain saturation $s^{n,\textnormal{a}}_{h,i}$ is piecewise constant and calculated using first order cell-centered finite volume method. The initial data $s_{i}^{0}$, $i\in\{1,2\}$, are discretized as follows:

[TABLE]

For the following time steps, we compute an intermediate saturation $s^{n,l}_{h,i}$ for all $l\in\{0,\dotsc,\mathcal{N}^{\textnormal{a}}_{i}-1\}$, by

[TABLE]

where $|\sigma|\varphi_{K,\sigma}^{n,l}$ is an approximation of the advection flux through the face $\sigma\in\mathcal{E}_{K}$, $\int_{\mathrm{\sigma}}f_{a,i}^{n}$, with $f_{a,i}^{n}={\bf f}_{i}(s)\cdot{\bf n}_{K,\sigma}=(f_{i}(s){\bf u}_{i}+f_{gi}(s){\bf u}_{g})\cdot{\bf n}_{K,\sigma}$ where ${\bf n}_{K,\sigma}$ denotes the outward normal to $\sigma$ with respect to $K$. The function $\varphi_{K,\sigma}^{n,l}$ will be defined as a function of the two values of the saturation on the two sides of $\sigma$. Let us first suppose the case of conforming time grids between the two rocks, i.e., $\mathcal{D}_{1}^{n}=\mathcal{D}_{2}^{n}$. The numerical flux $\varphi^{n,l}_{K,\sigma}$ in that case is calculated by

[TABLE]

where $\theta^{n,l}_{K,\sigma}$ and $\theta^{n,l}_{L,\sigma}$ are two additional unknowns on $\mathrm{\sigma}=K|L\in\mathcal{E}_{h}^{\mathrm{\Gamma}}$, such that $K\in\mathcal{T}_{h,1}$ and $L\in\mathcal{T}_{h,2}$, chosen to satisfy the following conditions

[TABLE]

We give now the definition of the flux function $\varphi_{i}$. Following [1, 48], and taking advantage that in our case $f_{ai}$ has a particular bell-shape, (see Figure 5), which has either one global maximum and no other local maxima or one global minimum and no other local minima, the function $\varphi_{i}(a,b)$, for $(a,b)\in[0,1]^{2}$ and $i\in\{1,2\}$, is given by:

When $f_{ai}$ has one maximum:

[TABLE]

When $f_{ai}$ has one minimum:

[TABLE]

In view of the above definition, the flux function $\varphi_{i}$ is exactly the Godunov numerical flux with respect to $f_{ai}$ [1, 12, 48].

Remark 6** (Phase-by-phase function)**

One can also use the phase-by-phase upstream flux function instead of the Godunov function, which is simpler for implementation, and whose expression relies on a particular structure of the flux function ${\bf f}_{i}$ (cf. [12, 14]).

To solve numerically the advection system (3.20) it is desirable to reduce the number of unknowns by eliminating the equation (3.20b). Following [12], finding $(\theta^{n,l}_{K,\sigma},\theta^{n,l}_{L,\sigma})$ solution of (3.20) can be reduced to the following problem.

Definition 7** (Conforming-in-time advection problem)**

At each time step $t^{n,l}$, find $q_{\sigma}^{n,l}\in{\mathbb{R}}$ such that

[TABLE]

This interface problem is with one implicit unknown per interface face (see Figure 6) and can be solved using Newton’s method or simply using a scalar root finder, e.g. Regula Falsi method (cf. [12]).

Now, it remains to extend the above setting to the nonconforming time grids. To this aim, we introduce the space-time variables $\theta^{n,a}_{K,\sigma}:=(\theta^{n,l}_{K,\sigma})_{0\leq l\leq\mathcal{N}^{\textnormal{a}}_{1}}$ and $\theta^{n,a}_{L,\sigma}:=(\theta^{n,l}_{L,\sigma})_{0\leq l\leq\mathcal{N}^{\textnormal{a}}_{2}}$, and similarly for $s^{n,a}_{K}$ and $s^{n,a}_{L}$. Therefore, based on the adoption of Robin interface conditions, i.e., proper linear combinations of the coupling conditions (3.20), we define the interface operator

[TABLE]

where $\mathrm{\kappa}_{1}$ and $\mathrm{\kappa}_{2}$ are strictly positive constants, and $\mathcal{P}_{12}$ and $\mathcal{P}_{21}$, are projection-in-time operators defined by (3.1) (see Subsection 3.2.4 for more details). The nonconforming-in-time counterpart of the flux continuity condition (3.21) is then replaced with a space-time formulation.

Definition 8** (Nonconforming-in-time advection problem)**

Find $(\theta^{n,a}_{K,\sigma},\theta^{n,a}_{L,\sigma})\in{\mathbb{R}}^{\mathcal{N}^{\textnormal{a}}_{1}}\times{\mathbb{R}}^{\mathcal{N}^{\textnormal{a}}_{2}}$ such that

[TABLE]

This interface system is solved iteratively with Newton’s method (cf. [2, 37] for a similar framework), where the interface unknowns $\theta^{n,a}_{K,\sigma}$ and $\theta^{n,a}_{L,\sigma}$ cooperate to solve the interface problem (3.25). Note that the operator within (3.21) is continuous and strictly monotone, therefore, this scheme admits a unique solution (cf. [12]). The space-time formulation (3.25) is not studied, but we can at least claim that the operator (3.24) inherits the monotonicity in function of each of the interface variables, as both the projection $\mathcal{P}_{ij}$ and intergration preserve that property. The stability of (3.25) is only verified numerically (see Section 4).

Remark 9** (Interface data)**

The couple $(\theta^{n,a}_{K,\sigma},\theta^{n,a}_{L,\sigma})$ on $\sigma=K|L\in\mathcal{E}_{h}^{\mathrm{\Gamma}}$ is the common saturation values, expressing the continuity of the capillary pressure across $\sigma$, while ensuring the flux transmission.

3.4.2 STEP 2: DIFFUSION_SOLVER

We inject the intermediate saturations $s^{n,*}_{h,i}:=s^{n,\mathcal{N}^{\textnormal{a}}_{i}}_{h,i}$ calculated through ADVECTION_SOLVER in the module DIFFUSION_SOLVER, to calculate the saturations at the coarse time step $n+1$. The scheme of this module is based on MFE scheme for the subdomain diffusion problems and the nonlinear Steklov-Poincaré and OSWR methods for the DD approaches, as used in the PRESSURE_SOLVER module. The reader is then invited to compare the two modules as the ideas and notations are analogous. First, we define for a specified capillary pressure $\lambda_{d}:=\pi_{i}(s^{n+1}_{h,i})|_{\Gamma}\in\mathrm{\Lambda}_{h}$, $i\in\{1,2\}$, the non-linear Dirichlet-to-Neumann operator:

[TABLE]

where $({\bf r}^{n+1}_{h,i},s^{n+1}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$, $i\in\{1,2\}$, is calculated by solving the diffusion problem inside each subdomain together with $\pi_{i}(s^{n+1}_{h,i})=\lambda_{d}$ on $\Gamma$, using MFE method (cf. [28, 55]);

[TABLE]

We then define the forms $d^{n+1}_{\Gamma,i}:\mathrm{\Lambda}_{h}\times\mathrm{\Lambda}_{h}\rightarrow{\amsmathbb R}$, $i\in\{1,2\}$, and $d^{n+1}_{\Gamma}:\mathrm{\Lambda}_{h}\times\mathrm{\Lambda}_{h}\rightarrow{\amsmathbb R}$:

[TABLE]

where $\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,i,n+1}:\mathrm{\Lambda}_{h}\rightarrow\mathrm{\Lambda}_{h}$, $1\leq i\leq 2$, and $\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,n+1}:=\sum_{i=1}^{2}\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,i,n+1}$ are non-linear Dirichlet-to-Neumann operators. Now, we are ready to rewrite the equations of this module as an interface problem posed only on $\Gamma$.

Definition 10** (Nonlinear Steklov-Poincaré problem)**

Find $\lambda_{d}\in\mathrm{\Lambda}_{h}$ such that

[TABLE]

Remark 11** (Equivalence with the multidomain diffusion problem)**

In Definition 10, the continuity of the capillary pressure across $\Gamma$, i.e., $[\![\pi(s^{n+1}_{h})]\!]=0$, is imposed in space by fixing $\lambda_{d}=\pi_{i}(s^{n+1}_{h,i})|_{\Gamma}$, $\forall i\in\{1,2\}$, while the continuity of the diffusive flux $[\![{\bf r}^{n+1}_{h}\cdot{\bf n}]\!]=0$ is retrieved weakly by solving (3.29).

Remark 12** (On Newton-Krylov method for (3.29))**

If we employ the Newton-GMRes method (cf. [40]) for solving (3.29), the Newton step is defined by $\lambda_{d}^{k+1}=\lambda_{d}^{k}+\eta^{k}$. Each step $\eta^{k}$ is computed by a forward difference GMRes iteration by solving the linear interface-Jacobian problem $\left(\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,n+1}\right)^{\prime}(\lambda_{d}^{k})\eta^{k}=-\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,n+1}(\lambda_{d}^{k})$. A well-known drawback of the GMRes algorithm for solving such interface problem is that the number of iterations depends essentially on the number of subdomain solves, therefore depends strongly on the subdomain discretization. A preconditioner is usually needed to reduce the number of iterations to a reasonable level. A left preconditioned GMRes strategy is based on solving $P^{-1}\left(\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,n+1}\right)^{\prime}(\lambda_{d}^{k})\eta^{k}=-P^{-1}\mathcal{Z}^{\textnormal{DtN}}_{\Gamma,n+1}(\lambda_{d}^{k}),$ where $P$ is an easily invertible approximation to the Jacobian. Physically, $P^{-1}$ can be interpreted as a non-linear Neumann-to-Dirichlet operator. Inverting this non-linear function would lead to a non-linear preconditioner and to resolve that one can construct an approximation of $P$ by solving a linear version of the problem as detailed in [60].

An alternative to the Steklov-Poincaré formulation (3.29), is to solve the diffusion problem in the module DIFFUSION_SOLVER with the OSWR scheme (see the linear case in PRESSURE_SOLVER) (cf. [2, 7]):

Algorithm 2** (The nonlinear OSWR)**

Choose an initial approximation $\mathrm{\eta}^{-1}_{i}\in\mathrm{\Lambda}_{h}$ of $\mathrm{\eta}_{i}:={\bf r}^{n+1}_{h,j}{\bf n}_{j}+\mathrm{\kappa}_{i}\pi_{j}(s^{n+1}_{h,j})$, with $j=(3-i)$, and a parameter $\mathrm{\kappa}_{i}\in{\amsmathbb R}^{+}\setminus\{0\}$. Define a tolerance $\epsilon>0$. Set $k:=-1$. 2. 2.

Do**

1. (a)

Increase $k:=k+1$. 2. (b)

Compute $({\bf r}^{n+1}_{h,i},s^{n+1}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$, $i\in\{1,2\}$, such that

[TABLE] 3. (c)

Set $\mathrm{\eta}^{k}_{i}:={\bf r}^{k,n+1}_{h,j}\cdot{\bf n}_{j}+\kappa_{i}\pi_{j}(s^{k,n+1}_{h,j})$, for $j=(3-i)$.

While* $\dfrac{\|(\mathrm{\eta}^{k}_{1},\mathrm{\eta}^{k}_{2})-(\mathrm{\eta}^{k-1}_{1},\mathrm{\eta}^{k-1}_{2})\|_{L^{2}(\Gamma)}}{\|(\mathrm{\eta}^{k-1}_{1},\mathrm{\eta}^{k-1}_{2})\|_{L^{2}(\Gamma)}}\geq\epsilon$.*

The well-posedness of the Robin subdomain problem as well as of the algorithm was shown in [7]. Note that the proof of convergence can be obtained only for the non-degenerate case, following the techniques of [35]. The free parameters $\mathrm{\kappa}_{i}$, $i\in\{1,2\},$ can be optimized to ensure faster convergence (cf. [6, 2, 35]).

Remark 13** (At convergence)**

At the convergence of Algorithm 2, the following Robin transmission conditions are retrieved

[TABLE]

These conditions are equivalent to the physical ones, i.e., $[\![{\bf r}^{n+1}_{h}\cdot{\bf n}]\!]=0$ and $[\![\pi(s^{n+1}_{h})]\!]=0$ on $\Gamma$.

3.5 The two-stage splitting algorithm

We provide the algorithm of Scheme 1, by assembling PRESSURE_SOLVER and saturation_SOLVER, which controls their execution and interaction.

Algorithm 3** (The overall algorithm of Scheme 1)**

Choose initial saturations $s^{0}_{h,i}\in M_{h,i}$, $i\in\{1,2\}$. Set $n:=-1$. 2. 2.

Do**

1. (a)

Increase $n:=n+1$. 2. (b)

Compute $({\bf u}^{n}_{h,i},p^{n}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$, $i\in\{1,2\}$, at the coarse pressure time step $t^{n}$, by performing a CG on (3.11), or by using Algorithm 1 (OSWR method). 3. (c)

Compute intermediate saturations $s^{n,l}_{h,i}\in M_{h,i}$, $i\in\{1,2\}$, through advection, at the finer advection time steps $t^{n,l}_{i}$, $1\leq l\leq\mathcal{N}^{\textnormal{a}}_{i}$, by performing Newton’s method on (3.25). 4. (d)

Set $s^{n,*}_{h,i}:=s^{n,\mathcal{N}^{\textnormal{a}}_{i}}_{h,i}$, $i\in\{1,2\}$. 5. (e)

Compute the saturation $s^{n+1}_{h,i}\in M_{h,i}$, $i\in\{1,2\},$ through diffusion, by performing a Newton-Krylov method on (3.29).

While* $n\leq N$.*

3.6 An improved IMPES method with multirate/nonconforming time grids

In this section, we review the classical IMPES method, where the problem is decoupled to pressure and saturation problems and then solved sequentially [25, 41]. Thus, we replace the module SPLIT_SATURATION_SOLVER in Scheme 1 with a specialized solver of the coupled advection-diffusion problem (ONE_SATURATION_SOLVER). This second approach is called Scheme 2. The flowchart of Scheme 2 can be drawn with simple modifications on Scheme 1. The drawback of this IMPES based scheme is that extremely small time steps should be used for overcoming the stability restriction. In order to accelerate the performance of the Scheme 2 and make it capable of solving larger problems (1) the pressure (fast solution) can cope with a much coarser time step compared to the saturation (slow solution) (2) the advection-diffusion problem for the saturation allows nonconforming time steps between different rock types.

3.6.1 The module ONE_SATURATION_SOLVER

The module ONE_SATURATION_SOLVER takes as input the pressures from the previous coarse time step $n$, and solves the saturation problem (Eqs. (2.7a)–(2.7b) with the IC (2.8)) for the finer time steps $t^{n,l}$, $l\in\{1,\dotsc,\mathcal{N}^{\textnormal{a}}_{i}\}$, $i\in\{1,2\}$. The scheme uses an explicit-in-time finite volume method for the discretization and the OSWR algorithm for the domain decomposition. To introduce the FV scheme, we first define the numerical flux $F^{n,l}_{K,\sigma}$ over a face $\sigma\in\mathcal{E}_{h}\cap\mathcal{E}_{K}$, $K\in\mathcal{T}_{h,i}$, by

[TABLE]

Therein, $\varphi_{K,\sigma}^{n,l}$ is the Godunov flux function given by (3.19) and where $\theta^{n,l}_{K,\sigma}$ is the unknown face saturation on $\sigma\in\mathcal{E}_{K}\cap\mathcal{E}_{h}^{\mathrm{\Gamma}}$. The numerical flux function $F_{K,\sigma}^{n,l}$ is then the sum of a diffusion contribution due to the capillary effects and an advection contribution due to the gravity effects and the total flow rate. The FV scheme for the multi-domain saturation is given by

[TABLE]

and for $l\in\{1,\dotsc,\mathcal{N}^{\textnormal{a}}_{i}\}$, the discrete saturations $s^{n,l}_{h,i}$, $i\in\{1,2\}$, is computed by

[TABLE]

and such that the saturation traces $\theta^{n,l}_{K,\sigma}$ and $\theta^{n,l}_{L,\sigma}$ on both sides of $\sigma=K|L$, $K\in\mathcal{T}_{h,1}$, and $L\in\mathcal{T}_{h,2}$, satisfies the following conditions

[TABLE]

where $\mathrm{\kappa}_{i}$, $i\in\{1,2\},$ are strictly positive constants. The reader is invited to compare the FV scheme (3.36)–(3.38) with that presented within ADVECTION_SOLVER, as the ideas are analogous and many of the explanations given there are carried over directly to the present case (see slo [7]).

Remark 14** (On equations (3.38))**

The equations (3.38) are simply a FV discretization over the time windows $\textnormal{on }\Gamma\times I^{n}$ of the Robin transmission conditions

[TABLE]

3.6.2 The OSWR Method

The solution of the FV scheme (3.36)–(3.38) is computed using the OSWR algorithm [2, 35]: At the iteration $k\geq 1$, for $i\in\{1,2\}$, let $\left(\theta_{h,i}^{0,n,l}\right)_{0\leq l\leq\mathcal{N}^{\textnormal{a}}_{i}}$ be given, we then compute the couple $(s^{k,n,l}_{h,i},\theta_{h,i}^{k,n,l})$ for all $l\in\{1,\dotsc,\mathcal{N}^{\textnormal{a}}_{i}\}$, by solving the following problem: for $i\in\{1,2\}$,

[TABLE]

with the transmission conditions

[TABLE]

for $j=(3-i)$, and $L\in\mathcal{T}_{h,j}$. For an overview and further details on the OSWR scheme for an explicit-in-time scheme, we refer the reader to [37].

3.6.3 The modified IMPES

The algorithm of Scheme 2 assembles PRESSURE_SOLVER and ONE_SATURATION_SOLVER in a coupling algorithm.

Algorithm 4** (The overall algorithm of Scheme 2)**

Choose initial saturations $s^{0}_{h,i}\in M_{h,i}$, $i\in\{1,2\}$. Set $n:=-1$. 2. 2.

Do**

1. (a)

Increase $n:=n+1$. 2. (b)

Compute $({\bf u}^{n}_{h,i},p^{n}_{h,i})\in\mathbf{W}_{h,i}\times M_{h,i}$, $i\in\{1,2\}$, at the coarse pressure time step $t^{n}$, by performing a CG method on (3.11), or by using Algorithm 1 (OSWR method). 3. (c)

Compute the saturation $s^{n,l}_{h,i}\in M_{h,i}$, $i\in\{1,2\}$, at the finer saturation time steps $t^{n,l}_{i}$, $l\in\{1,\cdots,\mathcal{N}^{\textnormal{a}}_{i}\}$, by performing the OSWR algorithm (3.40)-(3.41). 4. (d)

Set $s^{n+1}_{h,i}:=s^{n,\mathcal{N}^{\textnormal{a}}_{i}}_{h,i}$.

While* $n\leq N$.*

4 Numerical results

In this section, we discuss numerical solutions to the incompressible two-phase flow problem in three space dimensions using Algorithms 3 and 4. In many examples, we use the tolerances listed in Table 1. In all the experiments presented here the absolute permeability tensor $\mathbf{K}$ is actually a scalar value $K$, and the Mualem–Van Genuchten model is used for the relative permeability and capillary pressure curves [22]:

[TABLE]

with $n=2.8$, $m=1-1/n$ and $p_{\pi}=\mathrm{\alpha}\sqrt{\mathrm{\Phi}/K}$, where $\alpha$ is a proportionality constant or the Van Genuchten factor. Note that a small value of the Van Genuchten factor $\alpha$ indicates an advection-dominated problem. All of our tests include gravity effects. It is also worth noting that in all tests, the advective flux due to total flow rate dominates the one due to the gravity effects.

Remark 15** (Adaptive time stepping)**

In Algorithm 3 and 4, one can use an adaptive time stepping strategy in order to systematize the determination of adequate coarse and fine time steps for each problem. This strategy can be build based on a posteriori error estimates and distinguishing the various error components as used in the multirate scheme in [8].

Remark 16** (On the implementation)**

The two algorithms are implemented in the Matlab Reservoir Simulation Toolbox [43, 44]. The meshes are produced using the three-dimensional surface meshing software BLSURF interfaced with GHS3D [29] software for the three-dimensional volumetric meshes. The implementation uses the automatic differentiation feature of the toolbox to compute the Jacobian matrices for solving the nonlinear diffusion subdomain problems by the Newton method (cf. [44]).

4.1 Test case 1: saturation-diffusion problem between two rock types

In the first case, we particularize the two-phase flow model presented above for the sole saturation-diffusion equation (the total velocity being neglected) [7]. The goal of this test case is to assess and validate the module DIFFUSION_SOLVER. We study the convergence behavior of the Newton-GMRes method applied on (3.29) and the OSWR method given by Algorithm 2. We fix $T=250$ s and let $\mathrm{\Omega}=[0,1]^{3}$ be decomposed into two subdomains with two rock types (see Figure 7).

For the spatial discretization, we use a uniform tetrahedral mesh with 48000 elements. In (4.1), we let $\alpha=100$ for Rock 1 and so is $p_{\pi,1}=4.47$ psi, and then $\alpha=140$ so that $p_{\pi,2}=3.52$ psi for Rock 2. For the two rocks, we set $\mathrm{\Phi}=0.25$, and $K=1$. The saturation is set equal to $0.97$ on $\mathrm{\Gamma}^{\text{in}}:=\{x=0\}$. On the outflow boundary $\mathrm{\Gamma}^{\text{out}}:=\{x=1\}$, the saturation at time $t^{n+1}$ is set equal to that inside the closest cell at time $t^{n}$ (cf. [2, 5]). We assume homogeneous Neumann boundary conditions on the remaining part of the boundary. The initial condition is taken to be 0.95 in $\{(x,y,z)\in\mathrm{\Omega}\,\,|\,\,x<0.45\}$ and zero elsewhere which fits the continuity of the capillary pressure at the interface between the rocks.

The evolution of the saturation at two time steps is shown in Figure 8 (top). We remark a very sharp, discontinuous change in saturation at the interface between the two rocks. Clearly, the gas cannot penetrate to the domain $\mathrm{\Omega}_{2}$ with the same intensity due to the change in the capillary pressure function. Figure 8 (bottom) shows the capillary pressure at two time steps. We see that, unlike the saturation, the capillary pressure is continuous at the interface between the two rocks, highlighting numerically that the transmission conditions are satisfied.

We come now to the convergence analysis of the iterative methods involved by DIFFUSION_SOLVER (Newton-GMRes vs OSWR). First, we compare in Table 2 the results obtained with the Newton-GMRes method without and with preconditioner (see Remark 12). We then present in Table 3 the results obtained using the OSWR method.

We see that the Newton-GMRes method with preconditioner improves the convergence rate compared to the case with no preconditioner, i.e., the average number of Newton iterations is reduced from 6.36 to 4.11. That of the GMRes method, the average number of iterations is reduced from 16.72 to 11.2 iterations. In terms of number of subdomain solves, the results shows that the method without a preconditioner needs more subdomain solves. However, we note that each preconditioned GMRes iteration costs twice as much as one iteration of the method without preconditioning (to construct the preconditioner $P$).

For the results obtained with the OSWR method, the method is slower than the Newton-GMRes method and the average number of OSWR iterations is 8.15 iterations. This can be explained by the fact the free parameters in (3.30) are not the optimal ones. In terms of number of subdomain solves, the results are comparable with the method without preconditioning, even though OSWR method needs different subdomain solves (with nonlinear robin boundary condition).

4.2 Test case 2: fully two-phase flow between two rock types

We consider a numerical experiment from [9] given by the displacement of a non-wetting fluid by a wetting fluid in a domain $\mathrm{\Omega}=[0,10]^{3}$ made up of two different rock types and for $T=5\cdot 10^{6}$s. We let Rock 1 have an absolute permeability equal to $1$ millidarcies and five times larger than that on Rock 2. The porosity is fixed to $0.3$ in Rock 1 and $0.7$ in Rock 2. The injection boundary is taken orthogonal to the interface $\mathrm{\Gamma}$ between the two rock types. The initial saturation is set equal to $0.05$ in Rock 1 while in Rock 2 it is set to satisfy the equality of the capillary pressure on the interface.

4.2.1 Computational performance of Algorithm 3

We evaluate here the computational performance of Algorithm 3. The coarse time steps within the modules DIFFUSION_SOLVER and PRESSURE_SOLVER are fixed in size, i.e., $\tau^{n}=1\cdot 10^{4}$s. The finer time steps for the advection in the subdomains are also fixed in size and we choose $\tau^{n,l}_{1}=5\cdot 10^{2}$s in Rock 1 and $\tau^{n,l}_{2}=1\cdot 10^{3}$s in Rock 2. We choose the tolerances for the various algorithms involved in Algorithm 3 from Table 1.

In Figure 9, the saturation solution is depicted for four time steps, and we can show that the saturation is discontinuous across the interface so that it respects the continuity of the capillary pressure. We can see also that because of the contrast in the capillary pressure field, the saturation front moves faster around Rock 2. Thus, the capillary pressure has smoothed out a finger effect and a spike effect is observed at the interface. In Figure 10, we show the velocity profile between the rock types. We can see that the fluid inside Rock 2 snakes around the interface to travel again through Rock 1.

We come now to the convergence analysis of the iterative solvers involved in Algorithm 3. We start with the linear ones used in PRESSURE_SOLVER. In Figure 11, we plot the residual error for the CG solver with and without preconditioning and the OSWR method in one fixed time step $t^{n}=2\cdot 10^{6}$ (left) and the cumulative number of CG iterations as a function of time (right). One can clearly observe that the effect of the preconditioner on the convergence of the CG method is significant and that the desired residual tolerance is achieved with 6 iterations. The average number of CG iterations is reduced from 15.1 to 4.8. That of the OSWR method (Algorithm 1) is achieved after 21 iterations. For the solvers within the module DIFFUSION_SOLVER, the preconditioner of the GMRes iterations, used in the Newton-GMRes method, improves the linear convergence; we obtain the relative tolerance within 15 iterations (see Figure 12 (left)). Note that the reduction in the number of interface iterations is still small and the average is reduced from 17.81 to 13.87, but faster Newton convergence is observed during most of the iterations (see Figure 12 (right)). The OSWR method (Algorithm 2) for this step behaves as in the pressure problem so it converges slower than the Newton-GMRes method. For the Newton solver within ADVECTION_SOLVER, the cost is negligible as the solution of the interface advection problem (3.25) does not require subdomain solves, which is not the case in the other modules of the algorithm.

4.2.2 Comparison of the algorithms: accuracy-in-time and overall cost

Due to the nonconformity in time and the use of different splitting techniques, we want to see whether nonconforming time grids between subdomains as well as multirate time steps (within Algorithm 3 and 4) preserve the accuracy in time.

To study the accuracy of the nonconformity in time between subdomains, we compute with each algorithm a reference solution on a very fine time grid $(\tau^{n,l}_{i}=\tau^{n}/200=5\cdot 10^{1}$) and a fixed mesh. We then consider four initial time grids, refined then 4 times by a factor of 2:

• •

Time grid 1 (Conforming fine): $\tau^{n,l}_{i}=5\cdot 10^{2}$, $i\in\{1,2\}$,

• •

Time grid 2 (Nonconforming, fine in Rock 1): $\tau^{n,l}_{1}=5\cdot 10^{2}$, $\tau^{n,l}_{2}=1\cdot 10^{3}$,

• •

Time grid 3 (Conforming coarse): $\tau^{n,l}_{i}=2\cdot 10^{3}$, $i\in\{1,2\}$,

• •

Time grid 4 (Nonconforming, coarse in Rock 1): $\tau^{n,l}_{1}=2\cdot 10^{3}$, $\tau^{n,l}_{2}=2\cdot 10^{2}$.

In Figure 13, the error in the $L^{2}(0,T;L^{2}(\Omega))$-norm of the saturation versus the refinement level is depicted for the two algorithms. Clearly, both algorithms preserve the accuracy in time as the errors obtained in the nonconforming case with a fine time step in Rock 1 coincides with those obtained with the finer conforming case.

We come now to the comparison of Algorithm 3 and Algorithm 4 as well as the assessment of the accuracy in time of the multirate time strategies. In that case, we consider conforming time grids in the two algorithms. We compute using Algorithm 4 a reference solution on a very fine time grid $(\tau^{n,l}_{i}=\tau^{n}=2\cdot 10^{1})$ and a fixed mesh. We then test the two algorithms with $\tau^{n}=\mathcal{N}^{\textnormal{a}}\tau^{n,l}_{i}$, with a factor $\mathcal{N}^{\textnormal{a}}=48$ (the number of advection fine time steps within one coarse pressure and diffusion time step), divided then 4 times by a factor of 2.

In Figure 14, we show the error in the $L^{2}(0,T;L^{2}(\Omega))$-norm of the saturation versus the splitting level for the two algorithms. Clearly, an excellent quality of the solution is obtained from both algorithms even with $\mathcal{N}^{\textnormal{a}}=48$. Particularly, as stated previously Algorithm 3 decreases the number of inner time steps and Newton iterations required for the full saturation problem in Algorithm 4. As a result, the excellent quality of the solution from Algorithm 3 and its reduced cost compared to Algorithm 4 as well the reduced exchanged data (also smaller subdomain solves are required) between the different rocks make Algorithm 3 an efficient tool to deal with two-phase flow model in heterogeneous media.

4.3 Test case 3: extensions to multiple subdomains

Based on the above comparison, we choose here to apply Algorithm 3 on more complex configurations of porous media made-up with different rock types. For the next two cases, we let $\mathrm{\Omega}=[0,10]^{3}$ and $T=5\cdot 10^{6}$s.

4.3.1 Fractured porous media

In the first configuration, we consider the test case depicted in Figure 15 (left) where a three dimensional domain is divided into two equally sized subdomains by a fracture of the same size. Precisely, Rock 1 (rock matrix) appears in the left and right five layers of elements, separated by Rock 2 (fracture) in the center also five layers of elements.

In this example, the properties of Rock 1 are the same as the previous test case so that $K_{1}=1$ and $\phi_{1}=0.3$. In Rock 2, we fix the absolute permeability to be ten times larger than that of Rock 1. The coarser time steps are of fixed size and $\tau^{n}=1\cdot 10^{3}$s. The time step for the advection in the fracture is taken five times smaller than used in the surrounding rocks, i.e., $\tau^{n,l}_{1}=\tau^{n}/100$ and $\tau^{n,l}_{2}=\tau^{n}/20$. Water is injected uniformly through one side perpendicular to the fracture (see Figure 15). The production boundary is the opposite side. The other boundaries are impermeable.

Plots of the simulation results are shown in Figure 16. The results illustrate the discontinuous behavior of the saturation at rock interfaces due to the strong capillarity effects present in the fracture. The water saturation front snakes around the fracture to travel through the fracture and then moves from the injection boundary to the production boundary of the fracture. The convergence results (not shown here) are similar to what is observed for the previous test case, confirming the ability of Algorithm 3 to deal with more complex configuration of porous media. Particularly, the efficiency of the preconditioners to improve the convergence rate of the solvers compared to the ones without preconditioners is also confirmed in this test case.

4.3.2 Multiple rock types with strong heterogeneity

In this test case, we test the capability of Algorithm 3 to deal with multiple rock types with highly contrasting physical parameters. Precisely, we consider a domain made up with of four rock types (see Figure 17). The rock properties are listed in Table 4; we can see that Rock 1 is the most permeable and Rock 4 is the one with the lowest permeability. We also fix $\alpha$ in order to make an advection-dominated problem (with negligible capillarity effects) within Rock 1 and Rock 2, and in Rock 3 and Rock 4 we neglect totally the advection effects. In these rocks, we have only a saturation-diffusion problem as in Test case 1.

For the initial condition, we assume that the domain contains some quantity of water situated only within Rock 1 and Rock 2; we set $s=0.95$ in Rock 1 and elsewhere is set to satisfy the continuity of the capillary pressure. The time scales for the advection are taken different in the different rocks depending on the distribution of the absolute permeability, i.e., $\tau^{n,l}_{1}=\tau^{n}/10$, $\tau^{n,l}_{2}=\tau^{n}/5$, and $\tau^{n,l}_{3}=\tau^{n,l}_{4}=\tau^{n}/2$ with coarse and fixed time steps for diffusion and pressure, i.e., $\tau^{n}=\tau^{n}=5\cdot 10^{2}$s.

We show the results in Figure 18 for two time steps. One remark that a very sharp and discontinuous change in saturation at the rock type interfaces happens due to the higher contrast in the capillary pressure. The saturation of the wetting phase in the more permeable rocks (Rock 1 and 2) is increasing rapidly in these two rocks and a small quantity of the wetting phase penetrates into the Rocks 3 and 4 in which advection effects are neglected. We remark in Figure 18 that the wetting phase propagates in these subdomains with a finite speed due to their low permeabilities and due to the only diffusion effects that are present. In terms of computational effort, a good performance of the various inner algorithms is observed.

5 Conclusion

We propose in this paper a splitting-based domain decomposition methods to simulate two-phase flow in a porous medium composed of different rock types. The solution is resolved through a sequential approach that consists of splitting the original problem into three (or two) problems, where (global) pressure, saturation-advection and saturation-diffusion (or full saturation) problems are solved sequentially at each time step. The resulting schemes, that differ by how we treat the saturation problem (fully coupled or decoupled) provides us with flexible and efficient ways to treat the discontinuity of the saturation between rock types, as we can adapt the time scales for the advection and diffusion effects, as well as we can adapt the time scales for advection with respect to the rock type. Numerical experiments including those with several rock types and fractures showcase the computational efficiency of the methods and highlight their flexibility to handle complex transmission conditions between different rock types. Work underway addresses the stability of the derived schemes and provides adaptive multirate and local time steps based on a posteriori error estimates.

Acknowledgments

This research was partially funded by the Hydrinv Inria Euro Med 3+3: HYDRINV project. It has also received funding from EPIC project (within LIRIMA: http://lirima.inria.fr) and the Tunisian Ministry of Higher Education and Scientific Research. The author thanks Jérôme Jaffré and Jean E. Roberts for helpful discussions.

Appendix A Appendix: Application to a reduced fracture model between two rock types

A further important feature of the developed DD approaches is the ability to integrate models for reduced fractures for two-phase flow in a natural way. Precisely, we apply Algorithm 3 to the discrete fracture model presented in [5] in which a fracture is treated as an interface of dimension 1 in a 2-dimensional simulation, with fluid exchange between the 1-dimensional fracture flow and the 2-dimensional flow in the surrounding rock matrix (see Figure 19). Next, we give a short overview of the model and the used discrete scheme.

The model as presented in [5] is given by two-phase model problem (2.7) in each space–time domain $\mathrm{\Omega}_{i}\times(0,T)$ together with the following two phase flow in the fracture interface:

[TABLE]

where $\nabla_{{\tau}}$ denotes tangential component of the gradient operator, ${\bf u}_{g\tau}$ is the tangential component of ${\bf u}_{g}$ on $\mathrm{\Gamma}$, and where the functions $\mathbf{M}_{f}(\hat{s})$, $f_{gf}(\hat{s})$ and $\alpha_{f}(\hat{s})$ are defined using $\mathbf{K}_{f\tau}$ the tangential component of the permeability on the fracture $\mathbf{K}_{f}$ (see [5, 15, 45] for more details). We impose homogeneous Neumann boundary condition on $\partial\Gamma\times(0,T)$ and we assume that the initial saturation is known. For model closure, we introduce the following coupling conditions, for $i\in\{1,2\}$,

[TABLE]

The above systems are also coupled through the source terms appearing in the conservation equations in the fracture (see Eqs (A.1a) and (A.1d)), representing the difference between the fluid entering the fracture from one subdomain and that exiting through the other subdomain.

The scheme is given by extending Scheme 1 to the above setting, leading to solve sequentially reduced pressure, saturation-advection and saturation-diffusion interface problems posed only on the fracture. With the use of the definitions (3.26)–(3.28), we promptly arrive to rewrite the pressure problem (A.1d)-(A.1e) as an interface problem, in which simply we replace the pressure problem (3.9) by the following: at each coarse time step $n$, with known saturations $\hat{s}^{n}_{h}$, we solve for $\hat{p}^{n}_{h}\in\mathrm{\Lambda}_{h}$ such that

[TABLE]

The operator associated to this reduced problem is symmetric positive definite and can be solved using the CG method (cf. [4] for details). Now, we turn to the saturation problem. The flexibility in the time scales in the subdomains and in the fracture is a crucial asset in our numerical method, and it allows to significantly improve the accuracy of the scheme when highly permeable fractures are present between different rock types. The same notations to obtain (3.21) allows to take into account the advection effects on the fracture from the subdomains; we simply replace (3.21) by solving for $\theta^{n,l}_{f,h}:=\hat{s}^{n,l}_{h}\in\Lambda_{h}$, such that

[TABLE]

for all $l\in\{0,\cdots,\mathcal{N}^{\textnormal{a}}_{f}-1\}$, where $\mathcal{P}_{fi}$, $i\in\{1,2\}$, are projection-in-time type operators given by (3.1). Therein, $|e|\varphi^{n,l}_{\sigma,e}$ is an approximation of the advection flux through the edge $e$, $\int_{e}{\bf f}^{n}_{f}(\hat{s})\cdot\mathbf{n}_{e}$. Similarly to $\varphi^{n,l}_{K,\sigma}$ in (3.19), $\varphi^{n,l}_{\sigma,e}$ is a function of the two values of the saturation adjacent to $e=\sigma^{-}|\sigma^{+}$ and we calculate it using the introduced Godunov scheme with $i=f$. Now, it remains to extend the diffusion step. As for the reduced pressure problem, we replace (3.29) by solving $\theta^{n+1}_{f,h}\in\mathrm{\Lambda}_{h}$ such that

[TABLE]

This problem can be solved iteratively using fixed point iterations or the Newton method. The result of this system is then the saturation at the coarse time step $n+1$, i.e., $\hat{s}^{n+1}_{h}=\theta^{n+1}_{f,h}$.

Remark 17** (The interface preconditioners)**

To improve the convergence of the our DD method, the inverse of second order operator on the fracture $\nabla_{{\tau}}\cdot\left(-M_{f}(\hat{s}^{n}_{h})\nabla_{{\tau}}\cdot\right)$ is used as preconditioner for the interface pressure problem [17, 5, 36]. Similarly, a linearized version of the interface operator $\nabla_{\tau}\cdot\left(\mathbf{K}_{f\tau}\nabla_{\tau}(\alpha_{f}(\cdot)\right)$ is used as preconditioner for GMRes when Newton method is used to solve the interface diffusion problem (cf. [5, 4] for more details).

The test case we consider takes a unit square domain $\Omega$, divided into two equally sized subdomains by a fracture $\mathrm{\Gamma}$ located in the middle of the domain, perpendicular to the injection and production faces of $\mathrm{\Omega}$.

The permeability of the matrix is given by $\mathbf{K}=1$, very low compared to the permeability of the fracture $\mathbf{K}_{f}=1\cdot 10^{2}$. The porosity is equal to $\mathrm{\Phi}=0.1$ in the matrix and to $\mathrm{\Phi}=0.3$ in the fracture. We consider a triangular mesh with 800 grids. In time, we fix $T=5\cdot 10^{3}$, and use uniform coarse time steps in the fracture and subdomains, i.e., $\tau^{n}=T/100$. For the advection, we choose a finer and nonconforming time steps between the subdomains and the fracture; $\tau_{i}^{n,l}=\tau^{n}/5$, $i\in\{1,2\}$, and $\tau_{f}^{n,l}=\tau^{n}/20$ as depicted in Figure 19 (right). In Figure 20 (left), the saturation and the total velocity at $t^{n}=T/5$ is shown.

In this experiment, the flow is driven not only by the difference of the permeabilities and the capillary pressure fields, but also by the presence of the fracture. The wetting phase moves immediately through the fracture and the saturation front in the surrounding matrix snakes along the matrix-fracture interfaces. The wetting phase accumulates into the fracture near the production boundary until it eventually spreads out into the surrounding matrix. To verify the accuracy in time, we consider two initial time grids refined then 5 times by a factor of 2:

• •

Time grid 1 (Conforming fine): $\tau^{n,l}_{i}=\tau^{n}/5$, $i\in\{1,2\}$,

• •

Time grid 2 (Nonconforming, fine in fracture): $\tau^{n,l}_{i}=\tau^{n}/5$, $i\in\{1,2\},$ and $\tau^{n,l}_{f}=\tau^{n}/20$.

A reference solution is calculated with a very fine and single-rate time step and on fixed mesh with Algorithm 4. In Figure 20 (right) the error in the $L^{2}(0,T;L^{2}(\Gamma))$-norm of the saturation in the fracture versus the refinement level is depicted. As expected, first order convergence is ensured in the nonconforming case and the errors obtained in the nonconforming case with a finer time step in the fracture are nearly the same as in the finer conforming case. Thus, the use of nonconforming and multirate grids preserves the accuracy in time. Note that excellent results are obtained with a ratio 10 of the fine time step to the coarse time step.