Konstantin Brenner
email: konstantin.brenner@univ-cotedazur.fr

Nejmeddine Chorfi
email: nchorfi@ksu.edu.sa

R Masson

Sequential implicit Vertex Approximate Gradient discretization of incompressible two-phase Darcy flows with discontinuous capillary pressure

published or not. The documents may come

Sequential implicit Vertex Approximate Gradient discretization of incompressible two-phase Darcy flows

with discontinuous capillary pressure

Introduction

The numerical simulation of multi-phase Darcy flows plays an important role in many subsurface applications such as oil and gas recovery, basin modeling, geological storage, geothermal energy or hydrogeology [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF][START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media[END_REF][START_REF] Peaceman | Fundamentals of numerical reservoir simulation[END_REF][START_REF] Chen | Computational Methods for Multiphase Flows in Porous Media[END_REF]. This work focuses on incompressible and immiscible two-phase Darcy flows in heterogenous media with different rock types. The abrupt change of the pore sizes at different rock type interfaces induces spatially discontinuous capillary pressure curves at the Darcy scale triggering the discontinuity of the phase saturation and pressure solutions [START_REF] Van Duijn | The effect of capillary forces on immiscible two-phase flow in heterogeneous porous media[END_REF][START_REF] Enchéry | Numerical approximation of a two-phase flow problem in a porous medium with discontinuous capillary forces[END_REF][START_REF] Cancès | Finite volume scheme for two-phase flows in heterogeneous porous media involving capillary pressure discontinuities[END_REF][START_REF] Cancès | An existence result for multidimensional immiscible twophase flows with discontinuous capillary pressure field[END_REF][START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF]. These discontinuities play a major role and should be accurately captured in many important processes such as capillary driven imbibition in oil recovery or capillary barrier effects in oil migration and gas storage.

In the case of spatially homogeneous capillary pressures, the system of Partial Differential Equations accounting for incompressible two-phase Darcy flow can be typically expressed as an elliptic equation for the global pressure coupled to a degenerate parabolic equation for say the non-wetting phase saturation [START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media[END_REF]. The pressure equation is weakly coupled to the saturation via the total mobility function which varies smoothly in time. This naturally leads to sequential time integration schemes solving the linear pressure equation at fixed saturation followed by the solution of the saturation equation at given total velocity. The saturation equation can be integrated explicitly leading to the IMplicit in Pressure, Explicit in Saturation (IMPES) scheme or implicitly leading to the IMplicit in Pressure, IMplicit in Saturation (IMPIMS) or sequential implicit scheme [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF][START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media[END_REF][START_REF] Peaceman | Fundamentals of numerical reservoir simulation[END_REF][START_REF] Chen | Computational Methods for Multiphase Flows in Porous Media[END_REF]. Compared with Fully IMplicit (FIM) time integration schemes coupling the pressure and saturation unknowns with no time lag, sequential formulations have the advantage to allow for more adapted and efficient solvers for each sub-system with possibly different types of discretizations. Let us refer for example to [START_REF] Xie | Higher resolution total velocity vt and va finite-volume formulations on cell-centred structured and unstructured grids[END_REF][START_REF] Xie | Unstructured CVD-MPFA Reduced-Dimensional DFM Models for Two-Phase Flow, Coupled with Higher Resolution Hybrid Upwind Methods[END_REF] where a higher order discretization of the transport terms is used in the saturation equation or to [START_REF] Jenny | Adaptive multiscale finite-volume method for multiphase flow and transport in porous media[END_REF] where a multiscale method for the pressure equation is combined with a finite volume scheme on the fine scale for the saturation equation. In such formulations, let us remark that the use of the global pressure is not essential in practice. Natural (non-wetting and wetting) pressures can also be used at the price of an explicit approximation of the capillary pressure in the pressure equation, which is not an issue as long as homogeneous capillary pressures are considered. The key point of this type of splitting rather lies in the total velocity formulation combined with a conservative discretization of the total velocity fluxes. Let us also point out that this type of splitting strategy can be alternatively combined with an outer loop fixed point iteration leading to the so-called sequential fully implicit approach. If it converges, such approach can potentially keep the advantage of the sequential schemes, using specialized solvers, while converging to the fully coupled solution [START_REF] Moncorgé | Sequential fully implicit formulation for compositional simulation using natural variables[END_REF][START_REF] Moyner | A Mass-Conservative Sequential Implicit Multiscale Method for Isothermal Equation-of-State Compositional Problems[END_REF].

New difficulties arise to extend sequential integration schemes to the case of spatially discontinuous capillary pressure curves at heterogenous rock type interfaces. The transmission conditions at such interfaces strongly couple the pressure and saturation unknowns which raises the issue of the stability of the IMPES or the IMPIMS splitting schemes. As opposed to the case of homogeneous rock type, the global pressure formulation of the total velocity combined with an explicit approximation of the total mobility does not decouple the pressure equation from the saturation equation. This results from the saturation dependent jump of the global pressure at the heterogeneous rock type interfaces [START_REF] Cancès | An existence result for multidimensional immiscible twophase flows with discontinuous capillary pressure field[END_REF][START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF].

To our knowledge the stability of the IMPIMS scheme (or similarly of the IMPES scheme) has not been investigated in the case of heterogeneous rock types and is the main topic of this work. To fix ideas, our investigation is based on the Vertex Approximate Gradient (VAG) spatial discretization developed in [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF][START_REF] Brenner | Vertex approximate gradient scheme for hybrid dimensional two-phase Darcy flows in fractured porous media[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for two-phase Darcy flows in heterogeneous media on polyhedral meshes. Following [START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF], a total velocity formulation combined with an Hybrid Upwinding (HU) transport scheme [START_REF] Hamon | Implicit Hybrid Upwind scheme for coupled multiphase flow and transport with buoyancy[END_REF][START_REF] Hamon | Implicit hybrid upwinding for twophase flow in heterogeneous porous media with buoyancy and capillarity[END_REF][START_REF] Alali | Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions[END_REF][START_REF] Xie | Unstructured CVD-MPFA Reduced-Dimensional DFM Models for Two-Phase Flow, Coupled with Higher Resolution Hybrid Upwind Methods[END_REF][START_REF] Ahmed | A Cell-Centred CVD-MPFA Finite Volume Method for Two-Phase Fluid Flow Problems with Capillary Heterogeneity and Discontinuity[END_REF] will be used for the saturation equation. Thanks to the upwind positive discretization of the capillary and gravity flux terms, it is proved in [START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF] to guarantee the physical bounds on the saturations for any conservative discretization of the total velocity fluxes and whatever the time step. The FIM integration scheme combined with the VAG HU spatial discretization is described in [START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF] and will be used to provide reference solutions to be compared with the solutions of the sequential schemes investigated in this work. We will focus on sequential implicit integration schemes to avoid the strong limitations on the time steps induced by an explicit integration of the saturation equation in heterogeneous media.

All sequential implicit schemes investigated in this work satisfy the physical bounds on the saturations as well as the positivity of the capillary energy unconditionally on the time step. The schemes differ in the way the total velocity fluxes are computed. It will be shown numerically that the natural extension of the IMPIMS scheme, based on an explicit approximation of the total mobility and of the capillary pressure in the total velocity fluxes, suffers from instabilities for strongly contrasted capillary pressures. These instabilities increase for smaller time steps and consequently would also hold for an IMPES integration scheme. They are linked to the strong coupling in the transmission conditions at heterogeneous rock type interfaces between the pressure and saturation unknowns. A first tentative to circumvent this instability is an extension to the VAG discretization of a strategy commonly used for Two Point Flux Approximation (TPFA) [START_REF] Hamon | Implicit Hybrid Upwind scheme for coupled multiphase flow and transport with buoyancy[END_REF]. It is based on jumping over the interface using an harmonic averaging of the half transmissivities combined with an upwind approximation of the phase mobilities. The extension to the VAG discretization raises several issues addressed in this work. The second strategy investigated in this work extends an idea previously introduced in the case of an IMPES scheme in [START_REF] Andreianov | Approximating the vanishing capillarity limit of two-phase flow in multi-dimensional heterogeneous porous medium[END_REF] but to our knowledge not tested numerically. It amounts to maintain the coupling between the saturation and the pressure locally at heterogeneous rock type interfaces to compute the total velocity fluxes.

The remaining part of the paper is outlined as follows. The two-phase Darcy flow model with heterogeneous rock types together with its total velocity formulation is presented in Section 2. The VAG discretization in total velocity formulation following [START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF] is recalled in section 3 including the VAG HU implicit discretization of the saturation equation at given total velocity fluxes and the FIM discretization of the total velocity fluxes. Section 4 presents the usual IMPIMS discretization of the total velocity fluxes as well as the two strategies to fix its stability issues. In section 5, the three sequential implicit schemes are compared to the FIM VAG HU scheme both in terms of stability, accuracy and computational cost on three oil migration test cases with a large range of capillary pressure parameters in the drain and barrier rock types. The numerical investigation of the sequential schemes is extended to the cell-centered finite volume TPFA spatial discretization in subsection 5.4.

Two-phase Darcy flow model

Let Ω be a polyhedral bounded domain of R d (d ≥ 1) and (0, t f ) be the time interval. The wetting phase is indexed by w and the non-wetting phase by nw. To simplify the presentation, the porous medium domain Ω is constituted by two rock types rt ∈ RT = {b, d} partitioning the domain in two polyhedral open subsets (Ω rt ) rt∈RT such that Ω = rt∈RT Ω rt , where b typically denotes a barrier rock type and d a drain rock type. The interface ∂Ω b ∩ ∂Ω d between the two rock type subdomains is denoted by Γ. We consider the incompressible and immiscible two-phase Darcy flow model in heterogeneous porous medium defined by the following system of equations set on

Q t f = Ω rt × (0, t f )          φ rt (x)∂ t s α rt + divV α rt = 0, V α rt = -η α rt (s α rt )K rt (x)(∇p α rt -ρ α g), p c,rt = p nw rt -p w rt = P c,rt (s nw rt ), s nw rt + s w rt = 1, (1)
with α ∈ {nw, w}, rt ∈ RT . In (1) φ rt (x) denotes, for each rock type, the medium porosity, K rt (x) the medium permeability tensor, s α rt the phase saturation, p α rt the phase pressure, V α rt the phase Darcy velocity, and p c,rt the capillary pressure. The mass density ρ α of each phase α ∈ {nw, w} is assumed constant. The gravity acceleration vector is denoted by g and its norm by g. The phase mobility function, defined as the ratio of the relative permeability k α r,rt (s) to the constant phase dynamic viscosity µ α , is denoted by η α rt (s). Finally, P c,rt stands for the capillary pressure function, depending on the rock type (see the left Figure 1 for an example). In addition we introduce the monotone graph extension of P c,rt defined by P c,rt (s) = P c,rt (s), s > 0, (-∞, P c,rt (0)], s = 0.

The system (1) is complemented by transmission conditions at the interface Γ between the two rock types, involving the flux continuity equations

V α b • n| Γ = V α d • n| Γ , α ∈ {nw, w}, (2)
together with a formal phase pressure continuity condition. In view of [START_REF] Cancès | Finite volume scheme for two-phase flows in heterogeneous porous media involving capillary pressure discontinuities[END_REF], [START_REF] Cancès | An existence result for multidimensional immiscible twophase flows with discontinuous capillary pressure field[END_REF] and [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF], the formal continuity of the capillary pressure results in the following saturation jump relation

P c,b (s nw b | Γ ) ∩ P c,d (s nw d | Γ ) = ∅. (3)
Our numerical implementation of the condition (3) is based on the capillary pressure graph parametrization as introduced in [START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF]. Denoting by τ an auxiliary interface variable, we introduce a saturation function S nw rt (τ ) for each rock type rt as well as a capillary pressure function P c (τ ) which are such that The practical choice of this parametrization is based on a switch of variables technique between both rock type saturations and the capillary pressure (see Figure 1 for an example). This parametrization has the advantage to account for general capillary pressure curves including typically entry pressure in the matrix and vanishing capillary pressure in the drain zone. It also provides a more efficient nonlinear convergence than the more usual capillary pressure or saturation variable. Let us refer to [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for details.

P c (τ ) ∈ P c,rt (S nw rt (τ )), rt ∈ RT ,
To fix ideas, the system (1) is completed by an initial datum on the non-wetting phase saturation and by the boundary conditions

V α rt •n = 0 on Γ N ∩∂Ω rt ×(0, t f ), p α rt = p α Dir on Γ Dir ∩∂Ω rt ×(0, t f ) for α ∈ {nw, w},
where n is the unit normal vector to Γ N oriented outward to Ω, and {Γ N , Γ Dir } a partition of ∂Ω.

The discretization developed in the next section is essentially based on the total velocity formulation for which we recall the following definition of the total velocity:

V T rt = - α∈{nw,w} η α rt (s α rt )K rt (x)(∇p α rt -ρ α g). (4)
This allows to express the phase Darcy velocities at given total velocity V T as follows:

V nw rt = f nw rt (s nw rt )V T rt fractional flow term + η nw rt (s nw rt )η w rt (1 -s nw rt ) η T rt (s nw rt ) (-K rt (x)∇p c,rt ) capillary term + η nw rt (s nw rt )η w rt (1 -s nw rt ) η T rt (s nw rt ) (ρ nw -ρ w )K(x)g gravity term , (5)
and V w rt = V T rt -V nw rt with the non-wetting phase fractional flow function defined by

f nw rt (s) = η nw rt (s) η T rt (s) , (6)
and the total mobility by η T rt (s) = η nw rt (s) + η w rt (1 -s).

VAG discretization 3.1 VAG fluxes and porous volumes

The VAG discretization of two-phase Darcy flows introduced in [START_REF] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF] considers generalized polyhedral meshes of Ω in the spirit of [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF]. It is recalled in the following to fix the notations and for the convenience of the reader. Let M be the set of cells that are disjoint open polyhedral subsets of Ω such that K∈M k = Ω. For all k ∈ M, x k denotes the so-called center of the cell k under the assumption that k is star-shaped with respect to x k . We then denote by F k the set of faces of the cell k ∈ M and by F = k∈M F k the set of faces of the mesh. Remark that the faces are not assumed to be planar, hence the term "generalized polyhedral mesh". For σ ∈ F, let E σ be the set of edges of the face σ and V σ the set of vertices of σ. We denote by E = σ∈F E σ the set of all edges of the mesh. For each k ∈ M we define the set of nodes of the cell k by V k = σ∈F k V σ , and we also denote by V = K∈M V k the set of all vertices of the mesh, and by M s , the subset of cells sharing the node s ∈ V. It is then assumed that for each face σ ∈ F, there exists a so-called center of the face x σ ∈ σ \ e∈Eσ e such that x σ = s∈Vσ β σ,s x s , with s∈Vσ β σ,s = 1, and β σ,s ≥ 0 for all s ∈ V σ .

Moreover the face σ is assumed to be defined by the union of the triangles T σ,e defined by the face center x σ and each edge e ∈ E σ .

Note also that the mesh is supposed to be conforming w.r.t. the partition of Ω in subdomains Ω rt , rt ∈ RT , and w.r.t. the partition {Γ N , Γ Dir } of ∂Ω. We then denote by V Γ the set of nodes at the interface Γ and by V Dir the set of nodes located at the boundary Γ Dir .

The VAG discretization has been introduced in [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF] for diffusive problems on heterogeneous anisotropic media. It is based on the following vector space of degrees of freedom:

X D = {v ν ∈ R, ν ∈ D}, denoting by D the set D = M ∪ V.
A finite element discretization is built using a tetrahedral sub-mesh of M and a second order interpolation at the face centers x σ , σ ∈ F defined for all v D ∈ X D by

I σ (v D ) = s∈Vσ β σ,s v s .
The tetrahedral sub-mesh is defined by T = {T k,σ,e , e ∈ E σ , σ ∈ F k , k ∈ M} where T k,σ,e is the tetrahedron joining the cell center x k to the triangle T σ,e (see Figure 2). For a given v D ∈ X D , we define the function π T v D ∈ H 1 (Ω) as the continuous piecewise affine function on each tetrahedron of The VAG scheme is a control volume scheme in the sense that it amounts, for each d.o.f. not located at the Dirichlet boundary and each phase, to a volume balance equation. The two main ingredients are therefore the conservative fluxes and the porous volumes. For u D ∈ X D , the VAG fluxes F k,s (u D ) connect the cell k ∈ M to the nodes s ∈ V k (see Figure 2). They are derived from the finite element variational formulation (see [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF] for details) leading to the following generalized fluxes definition

T such that π T v D (x k ) = v k , π T v D (s) = v s , and π T v D (x σ ) = I σ (v D ) for all k ∈ M, s ∈ V, σ ∈ F.
F k,s (u D ) = k -Λ(x)∇π T u D • ∇ϕ s dx = s ∈V k T s,s k (u k -u s ), (7)
with

T s,s k = k Λ(x)∇ϕ s • ∇ϕ s dx.
Following [START_REF] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF], the porous volumes are obtained by distributing the porous volume of each cell k to its nodes s ∈ V k . For each k ∈ M, we define a set of non-negative volume fractions (a k,s ) s∈V k \V Dir satisfying s∈V k \V Dir a k,s < 1, and we set

φ k,s = a k,s k φ(x)dx.
Then, we consider for all k ∈ M

φ k = k φ(x)dx - s∈V k \V Dir φ k,s .
It is shown in [START_REF] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF] that the weights (a k,s ) s∈V k \V Dir can be chosen to avoid artificial drain enlargement at nodes sharing cells with highly contrasted permeabilities. Roughly speaking it suffices to select these weights proportional to the permeabilities of the cells k ∈ M s around a given node s.

Discrete physical unknowns

Let us define the set

D Γ = M ∪ (V \ V Γ ) = D \ V Γ .
The set of discrete physical unknowns (phase pressures, phase saturations, phase potentials and capillary pressure) are expressed in terms of the primary variables, which are chosen to be given by the set of non-wetting phase pressures p nw D = (p nw ν ) ν∈D , the set of non-wetting phase saturations s nw D Γ = (s nw ν ) ν∈D Γ and the set τ V Γ = (τ s ) s∈V Γ of parameters at each interface node. A rock type denoted by rt ν ∈ {b, d}, ν ∈ D Γ can be assigned to each cell and each node not belonging to the interface. Then, we define

                   p c,D = (p c,ν ) ν∈D , with p c,ν = P c,rtν (s nw ν ), ν ∈ D Γ , p c,s = P c (τ s ), s ∈ V Γ , p w D = (p w ν ) ν∈D , with p w ν = p nw ν -p c,ν , Φ α D = p α D + ρ α gZ D , with Z D = (z ν ) ν∈D , s w k = 1 -s nw k , k ∈ M, s nw k,s = s nw s , s ∈ V \ V Γ , k ∈ M s , s w k,s = 1 -s nw s , s ∈ V \ V Γ , k ∈ M s , s α k,s = S α rt k (τ s ), s ∈ V Γ , k ∈ M s . (8)
It results from these definitions that the pressures p α s at a given interface node s ∈ V Γ are single sided while the saturations s α k,s depend on the rock type of each cell k ∈ M s , capturing the saturation jump (see Figure 3).

General form of the discrete equations

For N ∈ N * , let us consider the time discretization

t 0 = 0 < t 1 < • • • < t n-1 < t n • • • < t N = t f of the time interval [0, t f ]. We denote the time steps by ∆t n = t n -t n-1 for all n = 1, • • • , N .
Before giving a detailed description of the discretization methods that we aim to compare in this article, let us present a general layout of the discrete equations. For all k ∈ M and s ∈ V k we denote by V T k,s the total velocity flux between the control volume associated with the cell k and the node s, similarly we denote by V nw k,s the non-wetting phase flux.

The total velocity satisfies the following system formed by the conservation equations and the boundary conditions

             s∈V k V T,n k,s = 0, k ∈ M, k∈Ms V T,n k,s = 0, s ∈ V \ V Dir , p nw,n s = p nw s,Dir , s ∈ V Dir , (9)
where p nw s,Dir are specified non-wetting phase pressures at the Dirichlet nodes. The conservation of the non-wetting phase is written as

                   φ k ∆t n s nw,n k -s nw,n-1 k + s∈V k V nw,n k,s = 0, k ∈ M, k∈Ms\V Dir φ k,s ∆t n s nw,n k,s -s nw,n-1 k,s - k∈Ms V nw,n k,s = 0, s ∈ V \ V Dir , s nw,n s = s nw s,Dir , s ∈ V Dir \ V Γ , τ n s = τ s,Dir , s ∈ V Dir ∩ V Γ , (10)
where s nw s,Dir and τ s,Dir are Dirichlet boundary conditions at the Dirichlet nodes. The discrete total velocity fluxes V T,n k,s are going to be constructed in terms of the local mobilities and the values of the non-wetting phase and capillary pressures; in particular we will have

V T,n k,s = V T,n k,s s nw, D , s w, D , p nw,n D , p c,D (11)
where the symbol may stand for n or n -1 independently in each term. Using V T,n k,s the non wetting fluxes will be defined in the form

V nw,n k,s = V nw,n k,s s nw,n k , s nw,n k,s , p n c,D , V T,n k,s . (12)
In the following we will consider several decoupling methods for the system ( 9) - [START_REF] Brenner | Convergence of a Vertex Centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF]. The construction of the discrete total velocity V T,n k,s will differ from one method to another. In contrast the definition of the discrete non-wetting phase velocity, detailed in the next section, will not be affected by the choice of the decoupling approach.

Saturation equation

In the context of the VAG discretization, as opposed to TPFA discretizations, the capillary gradient flux is not monotone [START_REF] Brenner | Convergence of a Vertex Centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF][START_REF] Eymard | TP or not TP, that is the question[END_REF]. Following [START_REF] Brenner | Vertex approximate gradient discretization preserving positivity for two-phase darcy flows in heterogeneous porous media[END_REF], this difficulty is overcome thanks to a first order upwind approximation of the mobility term. The same strategy is applied to the gravity flux. For a parameter F ∈ R, let us define

D k,s (s nw k , s nw k,s , F ) =              η nw rt k (s nw k )η w rt k (1 -s nw k,s ) η T rt k (s nw k ) if F ≥ 0, η nw rt k (s nw k,s )η w rt k (1 -s w k ) η T rt k (s nw k ) if F < 0. ( 13
)
Let us set G k,s = (ρ nw -ρ w )gF k,s (Z D ), and

D cap k,s = D k,s (s nw k , s nw k,s , F k,s (p c,D )), D g k,s = D k,s (s nw k , s nw k,s , G k,s
), let us assume that discrete total velocity VAG fluxes V T,n k,s satisfy (9), then the non-wetting phase VAG fluxes are defined by

V nw k,s (s nw k , s nw k,s , p c,D , V T k,s ) = f nw rt k (s nw k )(V T k,s ) + -f nw rt k (s nw k,s )(V T k,s ) - +D cap k,s F k,s (p c,D ) + D g k,s G k,s . (14)
Let us recall the results proved in [START_REF] Brenner | Vertex approximate gradient discretization preserving positivity for two-phase darcy flows in heterogeneous porous media[END_REF] at given total velocity fluxes satisfying the conservation equations ( 9) : (i) any solution of the saturation equation [START_REF] Brenner | Convergence of a Vertex Centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF] satisfies the physical bounds unconditionally on the time step, in the sense that the saturations are in [0, 1] provided that it is the case at initial time. Moreover, it satisfies the positivity of the capillary energy flux term in the sense that for all time step n one has

k∈M s∈V k D cap,n k,s F k,s (p n c,D )(p n c,k -p n c,s ) ≥ 0.
(ii) There exists a solution to (10).

Fully implicit scheme

The usual choice of the total velocity is based on a phase potential upwind approximation of the phase mobilities, leading to define

V T,up k,s s nw D , s w D , p nw D , p c,D , F nw κ,s , F w κ,s = η nw κ,s F k,s (p nw D + ρ nw gZ D ) + η w κ,s F k,s (p nw D -p c,D + ρ w gZ D ), (15)
with the upwind approximation

η α κ,s = η α rt k (s α κ ) if F α κ,s ≥ 0, η α rt k (s α κ,s ) else.
The Fully IMplicit (FIM) VAG HU scheme (see [START_REF] Brenner | Vertex approximate gradient discretization preserving positivity for two-phase darcy flows in heterogeneous porous media[END_REF]) writes as follows:

Find p nw,n D , s nw,n D Γ and τ n V Γ satisfying ( 9) and ( 10) with

V T,n k,s = V T,up k,s s nw,n D , s w,n D , p nw,n D , p n c,D , F k,s (Φ nw,n D ), F k,s (Φ w,n D ) , = α∈{nw,w} η α rt k (s α,n k ) (F k,s (Φ α,n D )) + -η α rt k (s α,n k,s )(F k,s (Φ α,n D )) -.
The FIM VAG HU scheme is investigated in [START_REF] Brenner | Vertex approximate gradient discretization preserving positivity for two-phase darcy flows in heterogeneous porous media[END_REF] where it is shown to be more robust than the more usual phase potential upwind discretization [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF]. It has the advantage to capture accurately the transmission conditions at the different rock types interfaces for general capillary pressure laws, to allow for large time steps and be easy to extend to more complex physics than immiscible incompressible flows. Its main drawback is to lead to coupled nonlinear and linear systems which can be difficult and costly to solve. It does not fully exploit the different nature of the pressure and saturation equations, not allowing the use of different time stepping or different spatial discretizations. This motivates the use of alternative time discretizations based on sequential implicit schemes decoupling the pressure and saturation equations as discussed in the following section.

Sequential implicit schemes

This type of splitting schemes combines sequentially at each time step, the computation of the total velocity fluxes solving a pressure system combined with the solution of the saturation system (10) at given total velocity fluxes. Note that thanks to the total velocity formulation combined with an upwind approximation of the mobility terms, the physical bounds on the saturations are garanteed whatever the time step.

In the case of a single rocktype, it is well known that an explicit approximation of the phase mobilities in the total velocity combined with a global pressure formulation [START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media[END_REF] of the total velocity leads to a rigorous and stable decoupling between the saturation and the pressure equations. It is usually prefered to keep the phase pressures rather than the cumbersome global pressure, which requires to approximate explicitly the capillary pressure in addition to the phase mobilities in the pressure equation. This does not seem to lead to instability issues and can be made rigorous in the case of a Two-Point Flux Approximation using a specific approximation of the total mobility following [START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF]. This type of decoupling leads to the so-called IMplicit in Pressure, IMplicit in Saturation (IMPIMS) scheme also termed sequential implicit scheme which is frequently used as an alternative to the fully implicit scheme for its potential lower computational cost and the flexibility to use different discretizations for each equation [START_REF] Chen | Computational Methods for Multiphase Flows in Porous Media[END_REF].

The case of different rock types raises an additional issue due to the strong coupling between the pressures and saturations in the transmission conditions at the interface Γ. As opposed to the case of a single rock type, the global pressure is discontinuous at the interface Γ and strongly coupled to the saturation [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF]. It will be shown in the numerical experiments that the usual IMPIMS scheme leads to severe instabilities in the case of strong capillary pressure effects at heterogeneous rock type interfaces. Two strategies are investigated to recover the stability of the splitting scheme based either on a jump over the interface strategy or on an additional coupling localized at the interface Γ.

IMplicit in Pressure, Implicit in Saturation (IMPIMS) scheme

The decoupling is obtained using explicit approximation of the phase mobilities and capillary pressures in the total velocity equation. In addition, compared to FIM, the upwind approximation of the phase mobilities in the total velocity fluxes is done w.r.t. the sign of the phase Darcy fluxes computed from the previous time step pressures in order to preserve the linearity of the pressure equation.

The IMPIMS sequential implicit scheme is defined as follows:

Step 1: Find p nw,n D , satisfying (9) with

V T,n k,s = V T,up k,s s nw,n-1 D , s w,n-1 D , p nw,n D , p n-1 c,D , F k,s (Φ nw,n-1 D ), F k,s (Φ w,n-1 D ) . (16)
Step 2: Find s nw,n D Γ and τ n V Γ satisfying [START_REF] Brenner | Convergence of a Vertex Centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF].

The main asset of this simple and usual decoupling strategy is to lead to a linear pressure equation which can be solved efficiently. On the other hand, it will be shown to lead to severe instabilities in the case of heterogeneous capillary pressures.

Partially decoupled scheme

In order to fix the stability issue of the IMPIMS splitting scheme we propose to account for the strong coupling of the pressures and saturations localized at the interface Γ. Following [START_REF] Andreianov | Approximating the vanishing capillarity limit of two-phase flow in multi-dimensional heterogeneous porous medium[END_REF], where this idea is introduced in the framework of an IMPES scheme (but to our knowledge not tested), the saturation equation at interface nodes is kept fully coupled with the total velocity conservation equations. It results that the capillary pressure is kept implicit at the interface using the auxiliary unknowns τ n s for s ∈ V Γ .

Let us decompose the set of the discrete capillary pressure variables as p c,D = p c,D Γ , p c,V Γ the modified sequential scheme will keep an implicit approximation of p c,V Γ . More precisely the modified sequential implicit scheme denoted by SIM-Γ writes as follows:

Step 1: Find p nw,n D and τ n V Γ = ( τ n s ) s∈V Γ , satisfying [START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] and

             k∈Ms\V Dir φ k,s ∆t n S nw rt k ( τ n s ) -S nw rt k ( τ n-1 s ) - k∈Ms V nw k,s s nw,n-1 k , S nw rt k ( τ n s ), p n-1 c,D Γ , P c τ n V Γ , V T,n k,s = 0, s ∈ V Γ \ V Dir τ n s = τ s,Dir , s ∈ V Dir ∩ V Γ . with V T,n k,s = V T,up k,s s nw,n-1 D , s w,n-1 D , p nw,n D , p n-1 c,D Γ , P c τ n V Γ , F nw κ,s , F w κ,s
. and

F nw κ,s = F k,s p nw,n D + ρ nw gZ D , F w κ,s = F k,s p nw,n D -p n-1 c,D Γ , P c τ n V Γ + ρ w gZ D .
Step 2: Find s nw,n D Γ and τ n V Γ satisfying [START_REF] Brenner | Convergence of a Vertex Centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF].

As opposed to the IMPES version introduced in [START_REF] Andreianov | Approximating the vanishing capillarity limit of two-phase flow in multi-dimensional heterogeneous porous medium[END_REF], the parameter τ s at the interface nodes s ∈ V Γ is recomputed at the saturation equation step using the total velocity fluxes computed at the pressure equation step (Step 1). Note that, as opposed to the previous IMPIMS scheme, the SIM-Γ scheme leads to a nonlinear pressure system (Step 1) which will reduce the potential gain in terms of computational cost compared with the fully implicit scheme. A potential gain still remains thanks to the reduced sizes of the Jacobian linear systems compared with the fully coupled Jacobian systems of the FIM algorithm. The availability of efficient preconditioners both for the pressure and saturation equations will consequently play a key role. For the saturation equation, a zero fill-in incomplete factorization (ILU0) provides a very efficient and cheap preconditioner which is a key asset of the sequential approach. For the pressure equation which concentrates the stiffness of the system, the CPR-AMG preconditioner [START_REF] Lacroix | Decoupling preconditioners in the implicit parallel accurate reservoir simulator (IPARS)[END_REF] considered as the state of the art preconditioner for the fully coupled system is a possible choice but does not seem yet optimal. A possible cheaper and still efficient choice could be based on an approximate Schur complement defined by the approximate diagonal elimination of the interface unknowns τ s , s ∈ V Γ . It could then be combined with an AMG preconditioner of the approximate Schur complement. Note that, in the case of the TPFA scheme instead of the VAG scheme, this Schur complement is exact. The investigation of such a customized preconditioner is left for future work.

Sequential implicit scheme based on Harmonic averaging at the interface

A common way to decouple the pressure and saturation equations in the case of Two-Point Flux Approximation (TPFA) is based on an harmonic averaging of the half transmissivities at the interface combined with an upwind approximation of the phase mobilities jumping over the interface [START_REF] Hamon | Implicit Hybrid Upwind scheme for coupled multiphase flow and transport with buoyancy[END_REF]. The extension of these ideas to the case of the VAG discretization first requires to introduce two independent pressure unknowns at each interface node s ∈ V Γ . These unknowns are used to express the linear VAG Darcy fluxes continuity equations at each interface node and for each phase; extending the TPFA harmonic averaging to MultiPoint Flux Approximations. Second, in order to ensure the conservativity of the total velocity fluxes, we need to define for each interface node, a single phase mobility obtained by averaging the mobilities over the upwind cells around the given node.

The set of discrete unknowns of the pressure system is now defined by the set of nonwetting phase pressures p nw,n D together with the set of wetting phase pressures at interface nodes (p w,n s ) s∈V Γ . Let us set

p w,n,n-1 D = (p w,n,n-1 ν ) ν∈D with p w,n,n-1 ν = p nw,n ν -p n-1 c,ν , ν ∈ D Γ , p w,n,n-1 s = p w,n s , s ∈ V Γ , and V nw,n k,s = F k,s (p nw,n D + ρ nw gZ D ), V w,n k,s = F k,s (p w,n,n-1 D + ρ w gZ D ).
A single explicit phase mobility must be defined at interface nodes. It is obtained using the following mean upwind approximation:

ηnw s = k∈Ms|V nw,n-1 k,s ≥0 η nw rt k (s nw,n-1 k ) #{k ∈ M s |V nw,n-1 k,s ≥ 0} , (17)
ηw s = k∈Ms|V w,n-1 k,s ≥0 η w rt k (s w,n-1 k ) #{k ∈ M s |V w,n-1 k,s ≥ 0} . ( 18
)
Then, the total fluxes are defined by

V T,n k,s = η nw rt k (s nw,n-1 k )V nw,n k,s + η w rt k (s w,n-1 k )V w,n k,s , for s ∈ V \ V Γ , k ∈ M s , and by V T,n k,s = ηnw s V nw,n k,s + ηw s V w,n k,s , for s ∈ V Γ , k ∈ M s .
The discrete system to compute the total velocity fluxes combines the total velocity flux conservation equations with each phase Darcy flux conservation at interface nodes:

                       s∈V k V T,n k,s = 0, k ∈ M, k∈Ms V T,n k,s = 0, s ∈ V \ (V Dir ∪ V Γ ), k∈Ms V α,n k,s = 0, s ∈ V Γ \ V Dir , α ∈ {nw, w}, p nw,n s = p nw s,Dir , s ∈ V Dir , p w,n s = p w s,Dir , s ∈ V Dir ∩ V Γ . (19)
Note that the total velocity flux conservation is ensured at interface nodes thanks to the definition of a single (depending only on s) phase mobility at each interface node. This scheme is denoted by SIM-H-1 in the following. As will be shown in the numerical experiments, this harmonic averaging of the phase pressures at interface nodes is not consistent in the case of a vanishing phase mobility, hence it should not be used in other equations than interface node flux conservation equations (note that this is always true for TPFA). A possible remedy is to modify the definition of the fluxes as follows for

s ∈ V \ V Γ , k ∈ M s : V nw,n k,s = F k,s (p nw,n D + ρ nw gZ D ), V w,n k,s = F k,s (p nw,n D -p n-1 c,D + ρ w gZ D ). ( 20
)
This modified scheme is denoted by SIM-H in what follows. Both schemes can be described as follows:

Step 1: Find p nw,n D and (p w,n s ) s∈V Γ satisfying [START_REF] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF].

Step 2: Find s nw,n D Γ and τ n V Γ satisfying [START_REF] Brenner | Convergence of a Vertex Centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF].

Let us remark that, as for the IMPIMS scheme, the pressure system ( 19) is linear, but on the other hand it incorporates additional pressure unknowns at interface nodes.

Comparison of the FIM, IMPIMS, SIM-Γ and SIM-H VAG schemes

The set of two rock types is defined by {b, d} with b the barrier rock type and d the drain rock type. The relative permeabilities are set to k α r,rt (s α ) = (s α ) 2 for both rock types rt ∈ {b, d}. The phase mass densities are set to ρ nw = 700 Kg.m -3 and ρ w = 1000 Kg.m -3 , and the phase dynamic viscosities to µ nw = 5 10 -3 Pa.s and µ w = 10 -3 Pa.s.

The test case considers the capillary pressure relations

P c,b (s nw ) = p e -b b log(1 -s nw ), P c,d (s nw ) = -b d log(1 -s nw ), (21)
parametrized by the entry pressure p e and b b for the barrier rock type and by b d for the drain rock type.

In the following simulations, the time stepping is defined by ∆t 1 = ∆t init and for all n ≥ 1 by ∆t n+1 = max(∆t max , 1.2∆t n ), in case of a successful time step ∆t n , and ∆t n+1 = ∆t n 2 , in case of non convergence of the Newton algorithm in 25 iterations.

The nonlinear systems are solved using a Newton algorithm. To obtain a more robust convergence of the nonlinear solver for the saturation system, a damping of the Newton step forces a maximum variation of the weighted discrete unknowns. The cell unknowns are eliminated without any fill-in before solving the linear systems using a GMRES iterative algorithm preconditioned by a zero fill-in incomplete factorization preconditioner. The GMRES stopping criterion on the relative residual is fixed to 10 -7 . The Newton solver is convergent if the relative residual is lower than 10 -4 or if the weighted maximum norm of the Newton increment is lower than 10 -4 . We denote by N t the number of successful time steps, by N Chop the number of time step chops, by N Newton the total number of Newton iterations per successful time step, and by N GMRES the average number of GMRES iterations per Newton iteration. For the sequential schemes, N Newton and N GMRES are defined separately for both the pressure and saturation solves. Finally, CPU (s) stands for the CPU time in seconds using an Intel Core i7 2.6 Ghz processor.

Oil migration in a 1D basin with one barrier layer

The basin domain is defined by Ω = (0, 80 m) × (0, 80 m) × (0, 800 m) and includes a barrier rock type at the center of the domain of width 100 m for z ∈ (350, 450) m. It is discretized using a uniform mesh of size n x × n y × n z with n x = n y = 4 and n z = 80 (see Figure 4). The absolute permeabilities are set to K b = 10 -15 m 2 and K d = 10 -13 m 2 , and the porosity is fixed to φ b = φ d = 0.2. The basin is initially water saturated and the non-wetting phase is injected from the bottom boundary with the Dirichlet boundary condition s nw = 0.5 during 400 years. Then, the bottom non-wetting phase saturation is set to zero for t > 400 years. The bottom non-wetting phase pressure is set to 8.1 MPa. Dirichlet boundary conditions are set at the top boundary with p w = 10 5 Pa and s nw = 0. Homogeneous Neumann boundary conditions are imposed at the lateral sides of the domain. The simulation time is t f = 4000 years and the time stepping is defined by the initial time step ∆t init = 30 days, and the maximum time step ∆t max = 1 year for t ≤ 200 years and ∆t max = 10 years for t > 200 years.

Figures 5 and7 compare respectively the non-wetting phase saturation at final time and the variable τ at the bottom interface as a function of time obtained with the FIM reference scheme and the IMPIMS and SIM-Γ sequential implicit schemes for different values of the capillary pressure parameters. The non-wetting phase saturation accumulates below the barrier during the injection phase and the leakage through the barrier at long times after the injection phase strongly depends on the capillary pressure contrast between the barrier and drain rock types. For small capillary pressure contrast (left plots in Figures 567), the solution is the same for all three schemes. A small oscillation of the parameter τ at the bottom interface is visible for the IMPIMS scheme with no significant impact on the saturation solution. These oscillations increase with higher capillary contrast (middle and right plots of Figure 7) despite the L ∞ stability of the saturation and parameter τ solutions. Note that these oscillations resulting from the explicit coupling of the pressure and saturation unknowns at the interface Γ actually increase for smaller time steps. As shown in the middle and right plots of Figure 5, it results that the IMPIMS scheme provides a wrong solution over-estimating the leakage across the barrier.

On the other hand the SIM-Γ scheme provides basically the same saturation solutions than the FIM scheme for all capillary pressure parameters. It shows that it suffices to maintain the pressure saturation coupling at heterogeneous rock type interfaces in the computation of the total velocity fluxes to prevent the oscillations and recover the right solution. We note still a small oscillation for the case with entry pressure (right plot in Figure 7) but with no significant impact on the saturation solution. Such oscillation is somewhat expected in the range (1, 2) of the parameter τ which corresponds to the vertical part of the barrier capillary pressure where the pressure and saturation unknowns are strongly coupled by the transmission conditions.

The comparison between the reference FIM scheme and the SIM-H-1 and SIM-H schemes is exhibited in Figures 6 and8. It is clearly seen that the SIM-H-1 scheme provides a wrong solution for the largest capillary pressure constrasts (middle and right plots in Figures 678) which is due to the MultiPoint nature of the VAG scheme combined with the inconsistent approximation of the phase pressures at interface nodes for vanishing mobilities. This issue is solved by the SIM-H version implementing the correction (20) of the total velocity fluxes. As for the SIM-Γ scheme, the SIM-H scheme provides a rather bad approximation of the bottom parameter τ in the range (1, 2) for the case with entry pressure due to the stiff pressure and saturation coupling.

Remark 5.1 Let us point out that an upwind approximation of the phase mobilities w.r.t. the current phase Darcy fluxes rather than the lagged ones in the definition of the total velocity [START_REF] Enchéry | Numerical approximation of a two-phase flow problem in a porous medium with discontinuous capillary forces[END_REF], leading to a nonlinear pressure equation, does not improve the stability of the IMPIMS algorithm. For the test case with nonzero entry pressure, we have also investigated the convergence of the IMPIMS splitting combined with an outer loop fixed point iteration using

V T,n,q k,s = V T,up k,s s nw,n,q-1 D , s w,n,q-1 D , p nw,n,q D , p n,q-1 c,D , F k,s (Φ nw,q-1 D ), F k,s (Φ w,q-1 D ) ,
with q the fixed point iteration count. We have noticed that the convergence of the outer loop fails at times when instabilities of the IMPIMS sequential scheme occur and that the solution is not improved compared with the sequential algorithm.

Oil migration in a 2D basin with 2 barriers

The basin domain is defined by Ω = (0, 800 m) × (0, 10 m) × (0, 800 m) discretized using a uniform mesh of size n x × n y × n z with n x = 40, n y = 1 and n z = 80. It includes two horizontal barrier layers exhibited in blue in Figure 9 10, all schemes provide basically the same solution with no significant differences. For the case with highly contrasted parameters b b = 6 10 5 , b d = 10 2 and no entry pressure in Figure 11, the IMPIMS scheme over-estimates the leakage across the barriers. This is also slightly the case of the SIM-H scheme while the solutions of the FIM and SIM-Γ schemes cannot be distinguished. For the case with entry pressure in Figure 12, the FIM and SIM-Γ solutions perfectly match while the SIM-H and IMPIMS schemes over-estimate the leakage across the barriers.

Tables 2-1-3 exhibit the numerical efficiency of the four schemes for each test case. For midly contrasted parameters b b = 6 10 5 , b d = 10 5 and no entry pressure, the IMPIMS and SIM-H schemes are the cheapest schemes thanks to the linearity of the pressure equations with a higher cost for the SIM-H scheme due to the additional pressure unknowns at heterogeneous rock type interfaces. The SIM-Γ scheme is somewhat more expensive than the SIM-H scheme since the pressure equation becomes nonlinear. It still provides a gain in CPU time of a factor larger than 2 compared with the FIM scheme thanks to the reduced size of the linear systems and the much lower condition number of the saturation system compared with the fully coupled system or the pressure system. The same conclusions hold for the two other test cases, except that the IMPIMS scheme exhibit a large number of time step failures and Newton iterations as a result of the oscillations in time at the different rock type interfaces. Table 4 exhibits the numerical behavior obtained when using CPR-AMG as preconditioner for the FIM systems and for the SIM-Γ pressure systems (ILU0 is still used for the SIM-Γ saturation systems). It shows that the gain in CPU time is as expected reduced but still reaches 50 pourcents. As discussed in Subsection 4.2, there is potentially room for improvement using a customized preconditioner of SIM-Γ pressure systems.

Oil migration in a 3D basin with 3 barriers

The 3D basin domain Ω ⊂ (0, 100 m) 2 is exhibited in Figure 14 with its hexahedral mesh.

The mesh includes degeneracies of some hexahedra due to erosion leading to the collapse of up to three vertical edges of the cell which are accounted for by the VAG discretization.

The three barriers are exhibited in red and meshed with only one layer of cells. Note that the barrier located in the middle of the domain extends on the full domain horizontally while the two others extend only on the right horizontal side of the domain. The absolute permeabilities are set to K b = 10 -14 m 2 for the barriers and to K d = 10 -12 m 2 for the drain zone, and the porosity is fixed to φ b = φ d = 0.2 for both rock types. The capillary pressures are again given by ( 21) with different values of the parameters used in the following test cases.

The basin is initially saturated by the wetting phase and the Dirichlet conditions s nw = 0.9, p w = 8 MPa are fixed at the portion of the bottom boundary located below the first barrier and defined by {(x, y) ∈ (0, 100) × (0, 100) | x 2 + y 2 ≤ 25} × {0}. Dirichlet boundary conditions are also imposed at the top boundary with s nw = 0 and p w = 8. 10 6 -ρ w gz Pa. Impervious boundary conditions are set at the remaining boundaries. The simulation time is t f = 10 years with time stepping defined by the initial time step ∆t init = 1 days, and the maximum time step ∆t max = 20 days.

From Figures 14-15-16 exhibiting the non-wetting phase saturation at final time, it is clear that the solutions obtained with the FIM and SIM-Γ schemes perfectly match while both the SIM-H and IMPIMS schemes over-estimate the leakage across the barriers. For the IMPIMS scheme, these discrepancies are explained by the instability induced by the decoupling of the pressure and saturation in the transmission conditions at the interfaces between different rock types. This instability generates oscillations in time and space which can be observed in the bottom right Figure 16. The SIM-H scheme tries to circumvent the decoupling issue by jumping over the interfaces. This strategy stabilizes the coupling but suffers from an inconsistent approximation of the interface pressures especially on coarse meshes which is particularly the case for this test case with a single layer of cells in the width of the barriers. On the over hand, the SIM-Γ scheme, which maintains the pressure saturation coupling localized at the heterogeneous interfaces in the computation of the total velocity fluxes, provides a stable and consistent solution.

Tables 5-6-7 exhibit the numerical behavior of the four schemes on the three capillary pressures data sets. It shows, as for the previous test case, that the IMPIMS scheme exhibits a quite high number of nonlinear iterations due to the oscillations in time at heterogeneous rock type interfaces. The SIM-H scheme is the cheapest thanks to its linear pressure approximation. The SIM-Γ scheme is somewhat penalized by the nonlinearity of the pressure equation localized at different rock type interfaces, but still outperforms the FIM scheme in CPU time by a factor 2 or larger.

Figure 17 and Table 8 exhibit the final non-wetting phase saturations and the numerical behavior of the FIM and SIM-Γ schemes obtained with a time stepping using larger time steps defined by ∆t init = 5 days and ∆t max = 100 days. The final saturations are still in good agreement with a slightly increased discrepancy as could be expected compared with the previous time stepping using smaller time steps. The gain in CPU time is still roughly of a factor 2. We note however that the Newton solver for the saturation equation of the SIM-Γ scheme has roughly twice as many iterations as the fully coupled FIM scheme. More advanced nonlinear solvers for this saturation equation need to be investigated for further improvement. Figure 18 shows the CFL numbers obtained with the two time steppings and the SIM-Γ scheme. This shows that the SIM-Γ scheme can accommodate large CFL numbers with still a good accuracy and a good numerical behavior. Note that these CFL numbers are computed based on the total velocity only.

Figure 13: Two views on the 3D basin domain with its hexahedral mesh (including degeneracies of some hexahedra due to erosion) with 31422 cells. The barrier layers are exhibited in red and the surrounding drain rock type in blue. Next, we present in Figure 20 and Table 9 the results obtained for the 2D test case of Subsection 5.2 with nonzero entry capillary pressure. As for the VAG spatial discretization, the TPFA IMPIMS scheme exhibits the same over-estimation of the leakage accross the barriers triggered by the instabilities at the highly contrasted rock type interfaces. The saturation solutions of the FIM and SIM-Γ schemes cannot be distinguished. The saturation solution obtained with the SIM-H scheme is much better than in the VAG case. It exhibits only a slight over-estimation of the leakage accross the barriers. Note that the saturation front exhibits a reduced numerical diffusion on the right upper part compared with the VAG scheme solutions. This is a consequence of the first order upwind approximation for the saturation equation combined with the connectivity of the TPFA fluxes which is aligned with the vertical velocity in this part of the domain, while the diagonal connectivity of the VAG fluxes is not. The 3D test case of Subsection 5.3 has not been simulated with the TPFA scheme since the mesh does not satisfy the orthogonal consistency property. Alternatively, we mimick the challenging thin barrier configuration of this 3D test case using the 2D test case configuration with a single cell in the thicknesses of the two barriers (see Figure 21). 10. It shows the perfect match between the FIM and SIM-Γ schemes while the IMPIMS scheme exhibits severe instabilities and the SIM-H scheme significantly over-estimates the leakage accross the bottom barrier.

Conclusions

This work investigates sequential implicit schemes for incompressible two-phase Darcy flow with discontinuous capillary pressure. The spatial discretization is based on the Vertex Approximate Gradient (VAG) scheme accounting for unstructured polyhedral meshes combined with an Hybrid Upwinding (HU) of the transport term and an upwind positive approximation of the capillary and gravity fluxes. All splitting schemes are built from the total velocity formulation of the two-phase flow model and only differ in the way the conservative VAG total velocity fluxes are approximated. It results that the physical bounds on the saturations and the conservation of both phases are guaranteed whatever the time steps. It is shown in this work that the classical IMPIMS scheme, based on an explicit approximation of both the total mobility and of the capillary pressure in the total velocity fluxes, exhibits severe instabilities for contrasted capillary pressures. These instabilities increase with reduced time steps and are due to the strong coupling of the pressure and saturation unknowns localized at the different rock type interfaces. Two strategies to solve this issue are investigated, the first is based on a jump over the interface idea classicaly used for TPFA. It has the advantage to maintain the linearity of the pressure equation but suffers from an inconsistency of the pressure approximation at interfaces. The numerical experiments show that it allows to recover the stability of the coupling but leads to an over-estimation of the leakage through the barrier rock type. The second idea developed in this work is to maintain the coupling of the pressure and saturation unknowns but localized at heterogeneous interfaces in the definition of the pressure system which consequently becomes nonlinear. This allows to retreive the same solutions than the ones obtained with the fully implicit scheme on all test cases with a significant gain in CPU time depending on the choice of the preconditioner. These conclusions seem to be in a large extent independent of the choice of the spatial discretization since they are shown to extend to the cell-centered finite volume TPFA spatial discretization.

  τ . It results that the saturation jump condition (3) at the interface Γ is replaced by s nw rt | Γ = S nw rt (τ ), rt ∈ RT . We will also use in the following, the notation S w rt (τ ) = 1 -S nw rt (τ ).

Figure 1 :

 1 Figure 1: Example of parametrization of the capillary pressure graphs based on a switch of variables between the drain and barrier non-wetting phase saturations and the capillary pressure. (Left): drain and barrier capillary pressures scaled by 5b b with P c,rt (s nw ) = p e,rt -b rt log(1 -s nw ) using b d = 10 4 Pa, b b = 10 5 Pa, p e,b = 5 10 4 Pa, p e,d = 0. (Right): Plots of S nw rt (τ ), rt ∈ {d, b} and of Pc(τ ) 5b b . The variable τ matches with the drain saturation for τ ∈ (0, τ 1 ), with the scaled capillary pressure for τ ∈ (τ 1 , τ 2 ), and with the barrier saturation for τ ∈ (τ 2 , τ 3 ), with τ 1 = 1 -b d p e,b , τ 2 = τ 1 + 1 + b d p e,b log(1 -τ 1 ), τ 3 = τ 2 + 1.

  The nodal finite element basis functions associated to this interpolation operator are denoted by ϕ ν , ν ∈ D.

Figure 2 :

 2 Figure 2: (Left): example of one cell k with a tetrahedron T k,σ,e of the sub-mesh T . (Right): example of VAG fluxes in a cell k.

Figure 3 :

 3 Figure 3: Example of a 2D mesh with the interface Γ in bold. (Left): primary unknowns at a given cell k ∈ M, node s ∈ V \ V Γ and interface node s ∈ V Γ . (Right): phase pressures and saturations at a given cell k ∈ M, node s ∈ V \ V Γ and interface node s ∈ V Γ .

Figure 4 :

 4 Figure 4: Cartesian Mesh with 4 × 4 × 80 cells of the domain Ω = (0, 80 m) × (0, 80 m) × (0, 800 m), the barrier (in blue) and the drain zone (in red).

Figure 5 :

 5 Figure 5: Non-wetting phase saturation at final time obtained with the FIM, SIM-Γ and IMPIMS VAG HU schemes with (left): p e = 0, b b = 6 10 5 and b d = 10 5 Pa, (middle): p e = 0, b b = 6 10 5 and b d = 10 2 Pa, (right): p e = 6 10 5 , b b = 10 5 and b d = 10 3 Pa.

Figure 6 :

 6 Figure 6: Non-wetting phase saturation at final time obtained with the FIM, SIM-H and SIM-H-1 VAG HU schemes with (left): p e = 0, b b = 6 10 5 and b d = 10 5 Pa, (middle): p e = 0, b b = 6 10 5 and b d = 10 2 Pa, (right): p e = 6 10 5 , b b = 10 5 and b d = 10 3 Pa.

Figure 7 :

 7 Figure 7: Variable τ at the bottom rock type interface as a function of time obtained with the FIM, SIM-Γ and IMPIMS VAG HU schemes with (left): p e = 0, b b = 6 10 5 and b d = 10 5 Pa, (middle): p e = 0, b b = 6 10 5 and b d = 10 2 Pa, (right): p e = 6 10 5 , b b = 10 5 and b d = 10 3 Pa.

Figure 8 :

 8 Figure 8: Variable τ at the bottom rock type interface as a function of time obtained with the FIM, SIM-H and SIM-H-1 VAG HU schemes with (left): p e = 0, b b = 6 10 5 and b d = 10 5 Pa, (middle): p e = 0, b b = 6 10 5 and b d = 10 2 Pa, (right): p e = 6 10 5 , b b = 10 5 and b d = 10 3 Pa.

  Figures 10-11-12 exhibit the non-wetting phase saturation obtained at final time using the FIM, SIM-Γ, SIM-H and IMPIMS schemes with different values of the capillary

Figure 9 :

 9 Figure 9: Cartesian Mesh with 40 × 1 × 80 cells of the domain Ω = (0, 800 m) × (0, 10 m) × (0, 800 m), the two barriers (in blue) and the drain zone (in red).

Figure 10 :Table 1 :

 101 Figure 10: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS VAG HU schemes with p e = 0, b b = 6 10 5 and b d = 10 5 Pa.

Figure 11 :Table 2 :

 112 Figure 11: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS VAG HU schemes with p e = 0, b b = 6 10 5 and b d = 10 2 Pa.

Figure 12 :Table 3 :Table 4 :

 1234 Figure 12: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS VAG HU schemes with p e = 3 10 5 , b b = 10 5 and b d = 10 3 Pa.

Figure 14 :Table 5 :

 145 Figure 14: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS VAG HU schemes with b b = 10 5 , b d = 10 4 Pa and p e = 0.

Figure 15 :Table 6 :

 156 Figure 15: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS VAG HU schemes with b b = 10 5 , b d = 10 3 Pa and p e = 0.

Figure 16 :Table 7 :

 167 Figure 16: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS VAG HU schemes with p e = 5 10 4 Pa, b b = 10 5 Pa and b d = 10 3 Pa.

Figure 17 :

 17 Figure 17: From left to right, non-wetting phase saturation at final time for the FIM and SIM-Γ schemes with the time stepping defined by ∆t init = 5 days and ∆t max = 100 days. The capillary parameters are set to p e = 5 10 4 Pa, b b = 10 5 Pa and b d = 10 3 Pa.

Figure 18 :Table 8 :

 188 Figure 18: CFL numbers obtained with the SIM-Γ scheme using the two time steppings defined either by ∆t init = 1 day and ∆t max = 20 days or by ∆t init = 5 days and ∆t max = 100 days. The capillary parameters are set to p e = 5 10 4 Pa, b b = 10 5 Pa and b d = 10 3 Pa.

Figure 19 :

 19 Figure 19: Non-wetting phase saturation at final time (left) and variable τ at the bottom rock type interface as a function of time (right) obtained with the FIM, SIM-Γ, SIM-H and IMPIMS TPFA HU schemes.

Figure 20 :Table 9 :

 209 Figure 20: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS TPFA HU schemes with p e = 3 10 5 , b b = 10 5 and b d = 10 3 Pa.

  The simulation time is t f = 400 years with time stepping defined by the initial time step ∆t init = 30 days, and the maximum time step ∆t max = 1 year. The results obtained with the capillary pressure parameters set to p e = 4 10 5 , b b = 10 5 and b d = 10 3 Pa are exhibited in Figure 22 and Table

Figure 21 :

 21 Figure 21: Cartesian Mesh with 40×1×80 cells of the domain Ω = (0, 800 m)×(0, 10 m)× (0, 800 m), the two thin barriers (in red) and the drain zone (in blue).

Figure 22 :

 22 Figure 22: From left to right and top to bottom, non-wetting phase saturation at final time for the FIM, SIM-Γ, SIM-H and IMPIMS TPFA HU schemes for the thin barrier 2D test case with p e = 4 10 5 , b b = 10 5 and b d = 10 3 Pa.

Table 10 :

 10 Numerical behavior of the four TPFA-HU schemes for the 2D basin test case with 2 thin barriers.

Acknowledgments

This project was funded by the National Plan for Science, Technology and Innovation (MAARIFAH), King Abdulaziz City for Science and Technology, Kingdom of Saudi Arabia, award number 14-MAT739-02.

Data availability Statement

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Comparison of the sequential schemes combined with the TPFA spatial discretization

The objective of this subsection is to investigate the potential extension of the previous analysis of the sequential schemes to another spatial discretization widely used for the simulation of two-phase Darcy flows in heterogeneous media, namely the cell-centered finite volume Two-Point Flux Approximation (TPFA). To be consistent with the previous framework, we consider here the TPFA discretization including the primary unknowns p nw σ , τ σ at each face σ ∈ F Γ shared by two cells k and l with different rock types, in addition to each cell k primary unknowns

Let us refer to [START_REF] Brenner | Total-velocity-based finite volume discretization of two-phase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF] for a detailed description of the TPFA HU discretization combined with the FIM integration scheme. The definitions of the corresponding TPFA sequential IMPIMS, SIM-Γ and SIM-H schemes are straightforward adaptions of their VAG counterpart, replacing the fluxes F k,s by either F k,σ if σ ∈ F Γ ∩ F k , and by F k,l for cells k and l sharing a face and having the same rock type. Note that SIM-H and SIM-H-1 fully match when combined with the TPFA discretization and that the upwinding of the mobilities ( 17)-( 18) in the total velocity flux corresponds to the usual upwinding between cells k and l sharing the face σ ∈ F Γ (see [START_REF] Hamon | Implicit Hybrid Upwind scheme for coupled multiphase flow and transport with buoyancy[END_REF] for this type of approach). For the TPFA discretization, our quick implementation of the sequential schemes is not optimized in terms of CPU time in the sense that it keeps the 2 by 2 Jacobian data structure for all sub-systems. It results that we only provide in the following tables a comparison of the nonlinear convergence behavior.

We first consider the 1D basin test case of Subsection 5.1 with nonzero entry capillary pressure and a modified final simulation time of 8000 years to better enhance the instability of the IMPIMS scheme. As exhibited in Figure 19, the results are basically in line with the ones obtained for the VAG schemes. The IMPIMS scheme exhibits severe oscillations of the parameter τ σ at the bottom interface and provides a wrong final saturation solution. On the other hand, the SIM-Γ and SIM-H schemes exhibit only a quite small oscillation of the parameter τ σ in the range [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Ahmed | A Cell-Centred CVD-MPFA Finite Volume Method for Two-Phase Fluid Flow Problems with Capillary Heterogeneity and Discontinuity[END_REF] with no significant differences on the saturation final solution compared with the FIM scheme.