K Brenner

J Hennicker
email: julian.hennicker@unige.ch

R Masson

Nodal Discretization of Two-Phase Discrete Fracture Matrix models

Introduction

Many real life applications in the geosciences like oil and gas recovery, basin modelling, energy storage, geothermal energy or hydrogeology involve two-phase Darcy flows in heterogeneous porous media. Such models are governed by nonlinear partial differential equations typically coupling elliptic and degenerate parabolic equations. Next to the inherent difficulties posed by such equations, further challenges are due to the heterogeneity of the medium and the presence of discontinuities like fractures. This has a strong impact on the complexity of the models, challenging the development of efficient simulation tools.

This work focuses on the numerical modelling of two-phase Darcy flows in fractured porous media, for which the fracture network is represented as a manifold of co-dimension one with respect to the matrix domain. These reduced models are obtained by averaging the physical unknowns as well as the conservation equations along the fracture width. They are termed hybrid-dimensional or also Discrete Fracture Matrix (DFM) Darcy flow models. Given the high geometrical complexity of real life fracture networks, the main advantages of these hybrid-dimensional compared with equi-dimensional models are both to facilitate the mesh generation and the discretisation of the model, and to reduce the computational cost of the resulting schemes. This type of hybrid-dimensional models is the object of intensive researches since the last 15 years due to the ubiquity of fractures in geology and their considerable impact on the flow and transport in the porous medium.

DFM models are closed with appropriate transmission conditions at matrix fracture (mf) interfaces which differ for fractures acting as drains or as barriers. For single-phase flows there are two major approaches. The first, designed for modelling highly conductive fractures and referred to as continuous pressure model [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF][START_REF] Brenner | Gradient discretization of Hybrid-Dimensional Darcy Flows in Fractured Porous Media[END_REF], assumes the continuity of the fluid pressure at the mf interfaces. The second approach, referred to as discontinuous pressure model [START_REF] Flauraud | Domain Decomposition for an asymptotic geological fault modeling[END_REF][START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF][START_REF] Jaffré | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Angot | Asymptotic and numerical modeling of flows in fractured porous media[END_REF][START_REF] Formaggia | A reduced model for Darcy´s problem in networks of fractures[END_REF][START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces[END_REF][START_REF] Chave | A Hybrid High-Order Method for Darcy Flows in Fractured Porous Media[END_REF], allows to represent fractures acting as permeability barriers by imposing Robin-type transmission conditions at mf interfaces.

When the modelling of two-phase flow is concerned, three major types of models can be distinguished. The first and most common type is based on a straightforward adaptation of the single-phase continuous pressure model to the two-phase setting (see [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for multiphase flow in fractured porous media[END_REF][START_REF] Monteagudu | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF][START_REF] Hoteit | An efficient numerical model for incompressible twophase flow in fracture media[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 twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF]), it assumes the continuity of each phase pressure at mf interfaces which allows to capture the saturation jump for fractures acting as drains and matrix as barrier. As for single-phase flow, this approach cannot account for fractures acting as barriers. In contrast to the single-phase context, let us stress that, due to heterogeneous capillary pressures, fractures having a large absolute permeability may still act as barriers for a given phase, typically for the wetting phase for fractures filled by the non-wetting phase (see [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF]). Another existing type of models, accounting for both drains or permeability barriers, is based on the linear (without mobility but including gravity) single-phase Darcy flux conservation equation imposed at mf interfaces for each phase. It is usually combined with Two-Point [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF][START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF] or Multi-Point [START_REF] Tunc | A model for conductive faults with non matching grids[END_REF][START_REF] Sandve | An efficient multi-point flux approximation method for Discrete Fracture-Matrix simulations[END_REF][START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF][START_REF] Gläser | A discrete fracture model for two-phase flow in fractured porous media[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][START_REF] Xie | Unstructured CVD-MPFA Reduced-Dimensional DFM Models for Two-Phase Flow, Coupled with Higher Resolution Hybrid Upwind Methods[END_REF] cell-centred finite volume schemes for which the interfacial discontinuous pressures are eliminated when building the single phase Darcy flux transmissibilities. These models account for the discontinuity of the pressures but not of the mobilities at mf interfaces. Both previous types of models are based on linear mf transmission conditions. The last type of models considers nonlinear mf transmission conditions which are based on the nonlinear (including mobility) two-phase normal flux continuity equations at mf interfaces. This type of models is considered in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF][START_REF] Aghili | Hybrid Finite Volume discretization of two-phase Discrete Fracture Matrix models with nonlinear interface solver[END_REF][START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Alali | Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions[END_REF] using a two-point flux approximation in the fracture width with upwinding of the mobilities, and in [START_REF] Jaffré | A discrete fracture model for two-phase flow with matrix-fracture interaction[END_REF][START_REF] Ahmed | A reduced fracture model for twophase flow with different rock types[END_REF] using a global pressure formulation. Such nonlinear transmission conditions account for the discontinuity of both the phase pressures and the mobilities at mf interfaces. A comparison of these three types of models using reference equi-dimensional solutions can be found in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF][START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF].

Having in mind that tetrahedral meshes are commonly used to cope with the geometrical complexity of fracture networks, nodal discretizations of DFM two-phase Darcy flow models have a clear advantage over cell-centred or face based discretizations thanks to their much lower number of degrees of freedom (d.o.f.). This is in particular the case when considering fully coupled implicit time integration which are necessary to avoid severe time step restrictions in high velocity regions such as fractures and to account for the strong coupling between the pressure and saturation unknowns at mf interfaces [START_REF] Andreianov | Approximating the vanishing capillarity limit of two-phase flow in multi-dimensional heterogeneous porous medium[END_REF]. Alternatively, cell centred discretizations have been considered for DFM two-phase flow models using the Two-Point Flux Approximation (TPFA) as in [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF][START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Alali | Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions[END_REF] or Multi-Point Flux Approximations (MPFA) as in [START_REF] Gläser | A discrete fracture model for two-phase flow in fractured porous media[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][START_REF] Xie | Unstructured CVD-MPFA Reduced-Dimensional DFM Models for Two-Phase Flow, Coupled with Higher Resolution Hybrid Upwind Methods[END_REF]. Face based discretizations have been considered in [START_REF] Hoteit | An efficient numerical model for incompressible twophase flow in fracture media[END_REF][START_REF] Ahmed | A reduced fracture model for twophase flow with different rock types[END_REF] using the Mixed Hybrid Finite Element (MHFE) method and in [START_REF] Groza | Modelization and discretization of two-phase flows in porous media with discrete fracture networks[END_REF][START_REF] Aghili | Hybrid Finite Volume discretization of two-phase Discrete Fracture Matrix models with nonlinear interface solver[END_REF] using the Hybrid Finite Volume (HFV) scheme. Non conforming discretizations have also been developed for this type of models using XFEM discretizations as in [START_REF] Fumagalli | A numerical method for two-phase flow in fractured porous media with non-matching grids[END_REF] or Embedded Discrete Fracture Models as in [START_REF] Tene | Projection-based Embedded Discrete Fracture Model (pEDFM)[END_REF].

Nodal discretizations, such as the Control Volume Finite Element (CVFE) method, have been first introduced in [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for multiphase flow in fractured porous media[END_REF][START_REF] Monteagudu | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF][START_REF] Geiger | Black-oil simulations for threecomponent, three-phase flow in fractured porous media[END_REF] for DFM two-phase Darcy flow models with continuous pressures at mf interfaces accounting for fractures acting as drains. In this work, we review the Vertex Approximate Gradient (VAG) discretization introduced in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Brenner | Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for continuous pressure models and in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF] for discontinuous pressure models. The VAG scheme is based on nodal d.o.f. like CVFE methods but it also includes the cell d.o.f. which are eliminated at the linear algebra level at each Newton iteration without any fill-in. These cell d.o.f. provide an additional flexibility in the design of the discretization allowing to cope with traditional issues raised at mf interfaces by nodal discretizations of the transport equation. On practical meshes, for which the cell sizes at mf interfaces are much larger than the fracture width, these issues are induced by the use of dual control volumes combined with heterogeneous petrophysical and hydrodynamical properties defined on the primal mesh.

The outline of the remaining of this article is as follows. Section 2 describes the DFM continuous and discontinuous pressure two-phase Darcy flow models as introduced in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF]. Section 3 presents the VAG discretizations of DFM continuous pressure twophase Darcy flow models. Several techniques to cope with the issues raised by nodal discretizations at mf interfaces are discussed, including the adaptation of the control volumes at mf interfaces, a new Multi-Point upwind approximation of the mobilities in Subsection 3.3, and taking into account the saturation jump for general capillary pressure curves in Subsection 3.5. Section 4 reviews the VAG discretizations of the three types of DFM discontinuous pressure two-phase Darcy flow models as presented in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF]. For each type of model and its VAG discretization, numerical experiments are exhibited on 2D and 3D DFM models including comparisons of the VAG discretizations to a face based scheme, as well as the comparison between the hybrid-dimensional DFM models and the reference equi-dimensional model.

Two-phase DFM discontinuous and continuous pressure models

Let Ω be a bounded domain of R d , d = 2, 3 assumed to be polyhedral for d = 3 and polygonal for d = 2. To fix ideas, the dimension will be fixed to d = 3 when it needs to be specified, for instance in the naming of the geometrical objects or for the space discretization. The adaptations to the case d = 2 are straightforward. Let Γ = i∈I Γ i denotes the network of fractures Γ i ⊂ Ω, i ∈ I, such that each Γ i is a planar polygonal simply connected open domain included in some plane of R d (see Figure 1). In the matrix domain Ω, we denote by φ m (x) the porosity and by Λ m (x) the permeabil-ity tensor. Along the fracture network x ∈ Γ, we denote by φ f (x) the porosity averaged on the fracture width and by d f (x) the fracture aperture. The permeability tensor is assumed constant along the width of the fracture and the normal vector to the fracture is assumed to be a principal direction. It results that we can define along the fracture network x ∈ Γ, the tangential permeability tensor Λ f (x) and the normal permeability λ n,f (x). It is assumed, for the sake of simplicity, that the matrix (resp. the fracture network) has a single rock type. Hence, for each phase α ∈ {nw, w} (where nw stands for the nonwetting phase and w for the wetting phase) we denote by M α m (s α ) (resp. M α f (s α )), the matrix (resp. fracture network) phase mobility, and by P c,m (s nw ) (resp. P c,f (s nw )), the matrix (resp. fracture network) capillary pressure function. The inverse of the monotone graph extension of the matrix (resp. fracture network) capillary pressure is denoted by S nw m (p) (resp. S nw f (p)). We will also denote by ρ α the phase density which for the sake of simplicity is assumed constant for both phases α ∈ {nw, w}.

Let α ∈ {nw, w}, we denote by u α m (resp. u α f ) the phase pressure and by s α m (resp. s α f ) the phase saturation in the matrix (resp. the fracture network) domain. The Darcy velocity of phase α ∈ {nw, w} in the matrix domain is defined by

q α m = -M α m (s α m )Λ m (∇u α m -ρ α g),
where g = -g∇z stands for the gravity vector with g the gravitational acceleration constant. The flow in the matrix domain is described by the volume balance equation

φ m ∂ t s α m + div(q α m ) = 0, (1)
for α ∈ {nw, w}, and the closure laws defined by the macroscopic capillary pressure law together with the sum to one of the phase saturations

s nw m = S nw m (p c,m ), p c,m = u nw m -u w m , s w m = 1 -s nw m . (2)
On the fracture network Γ, we denote by ∇ τ the tangential gradient and by div τ the tangential divergence. In addition, we can define the two sides ± of the fracture network Γ in Ω \ Γ and the corresponding unit normal vectors n ± at Γ inward to the sides ±. Let γ n ± (resp. γ ± ) formally denote the normal trace (resp. trace) operators at both sides of the fracture network Γ for vector fields in H div (Ω \ Γ) (resp. scalar fields in H 1 (Ω \ Γ). The Darcy tangential velocity of phase α ∈ {nw, w} in the fracture network Γ integrated over the width of the fracture is defined by

q α f = -d f M α f (s α f )Λ f (∇ τ u α f -ρ α g τ ),
with g τ = g -(g • n + )n + . The flow in the fracture network Γ is described, for each phase α ∈ {nw, w}, by the volume balance equation

d f φ f ∂ t s α f + div τ (q α f ) + γ n + q α m + γ n -q α m = 0, (3)
and by the closure laws

s nw f = S nw f (p c,f ), p c,f = u nw f -u w f , s w f = 1 -s nw f . (4)

Two-phase DFM discontinuous pressure model

We consider the transmission conditions introduced in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF]. They are based on a twopoint approximation of each phase normal flux within the fracture combined with a phase potential upwinding of the phase mobility taking into account the phase saturation jump at the mf interface. Let us first define, for both phases α ∈ {nw, w}, the "single" phase normal flux in the fracture network

V α,± f,n = λ f,n γ ± u α m -u α f d f /2 -ρ α g • n ± , (5)
which does not include the phase mobility. For any a ∈ R, let us set a + = max{0, a} and a -= min{0, a}. The conditions coupling the matrix and fracture unknowns then read, for α ∈ {nw, w} (see the right Figure 2):

γ n ± q α m = q α f,n ± , q α f,n ± = M α f (S α f (γ ± p c,m ))(V α,± f,n ) + + M α f (s α f )(V α,± f,n ) -. (6)
The hybrid dimensional two-phase flow discontinuous pressure model looks for u α m , u α f , s α m , s α f , α ∈ {nw, w}, satisfying (1)-( 2) and ( 3)-( 4) together with the transmission conditions (6).

Figure 2: (Left): example of a 2D DFM discontinuous pressure model with the normal vectors n ± at both sides of a fracture, the matrix phase pressure and saturation u α m , s α m , the fracture phase pressure and saturation u α f , s α f , the matrix Darcy phase velocity q α m and the fracture network tangential Darcy phase velocity q α f . (Right): illustration of the coupling condition q α f,n + = γ n + q α m for the hybrid-dimensional discontinuous pressure model.

Two-phase DFM continuous pressure model

In the case of pervious fractures, for which the ratio of the transversal permeability of the fracture to the width of the fracture is large compared with the ratio of the permeability of the matrix to the size of the domain, it is classical to assume that the phase pressures are continuous at the interfaces between the fractures and the matrix domain. Let us also mention that in the context of two-phase flows the continuous pressure DFM models have to be used with caution. It has been shown in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF] and [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF] that even highly pervious fractures may still act as barriers. This is due to the potential degeneracy of the mobilities in the transmission condition [START_REF] Alali | Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions[END_REF] and to the saturation jumps resulting from the high contrast of the capillary pressure curves across mf interface. Typically a fracture filled with the non-wetting phase would act as a barrier for the wetting phase, and therefore would induce a discontinuity of the wetting phase's pressure. We refer to [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF] and [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF] for a detailed comparison of continuous and discontinuous pressure models in case of very pervious fractures.

The continuous pressure model replaces the transmission condition (6) by the following phase pressure continuity conditions at mf interfaces:

γ + u α m = γ -u α m = u α f on Γ, α ∈ {nw, w}. (7)
It results that we can denote by u α the matrix pressure of phase α ∈ {nw, w} and by γu α the fracture pressure of phase α ∈ {nw, w}, where γ is the trace operator on Γ for functions in H 1 (Ω).

The hybrid dimensional two-phase flow continuous pressure model looks for s α m , s α f , and u α , α = nw, w satisfying (1)-( 2) and ( 3)-(4).

Figure 3: Example of a 2D DFM continuous pressure model with the normal vectors n ± at both sides of a fracture, the phase pressure u α and its trace γu α on the fracture network Γ, the matrix phase saturation s α m , the fracture phase saturation s α f , the matrix Darcy phase velocity q α m and the fracture network tangential Darcy phase velocity q α f .

For both continuous and discontinuous pressure models, a no-flux boundary conditions is prescribed at the tips of the immersed fractures, that is to say on ∂Γ\∂Ω, and the volume conservation and pressure continuity conditions are imposed at the fracture intersections. We refer to [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] and [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF] for more details on those conditions.

Finally, one should provide some appropriate initial and boundary data. To fix ideas, we consider in a non homogeneous Dirichlet boundary conditions on the matrix boundary ∂Ω Dir ⊂ ∂Ω and on the fracture boundary Σ Dir ⊂ ∂Γ ∩ ∂Ω. Homogeneous Neumann boundary conditions are set on ∂Ω N = ∂Ω \ ∂Ω Dir and on Σ N = (∂Γ ∩ ∂Ω) \ Σ Dir .

3 Vertex Approximate Gradient (VAG) discretization of two-phase DFM continuous pressure models

The VAG discretization of hybrid dimensional two-phase Darcy flows introduced in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] considers generalised polyhedral meshes of Ω in the spirit of [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF]. Let us briefly recall some notations related to the space discretization. We denote by M the set of disjoint open polyhedral cells, by F the set of faces and by V the set of nodes of the mesh. For each cell K ∈ M we denote by F K ⊂ F the set of its faces and by V K the set of its nodes.

Similarly, we will denote by V σ the set of nodes of σ ∈ F. The set M σ denotes the two cells sharing an interior face σ or the single cell to which the boundary face σ belongs. The set M s (resp. F s ) is the subset of cells (resp. faces) sharing the node s ∈ V.

Let E σ denote the set of edges of the face σ ∈ F. It is then assumed that for each face σ ∈ F, there exists a so-called "centre" 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 σ . The face σ is not necessarily planar, hence the term generalised polyhedral mesh. More precisely, each face σ is assumed to be defined by the union of the triangles T σ,e defined by the face centre x σ and each edge e ∈ E σ .

The mesh is supposed to be conforming w.r.t. the fracture network Γ in the sense that there exists a subset F Γ of F such that Γ = σ∈F Γ σ. We set

V Γ = σ∈F Γ V σ ,
and, for s ∈ V Γ , we define F Γ,s = F s ∩ F Γ as the subset of faces in F Γ sharing the node s.

The VAG discretization proposed in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] is based upon the following set of degrees of freedom (d.o.f.)

D = M ∪ V ∪ F Γ
and the corresponding vector space:

X D = {v ν ∈ R, ν ∈ D}.
The d.o.f. are exhibited in Figure 4 for a given cell K with one fracture face σ in bold.

Let us denote by

V Dir = {s ∈ V | x s ∈ ∂Ω Dir ∪ Σ Dir },
the subset of Dirichlet nodes.

A finite element discretization is built from the vector space of d.o.f. X D using a tetrahedral sub-mesh of M and a second order interpolation at the face centres x σ , σ ∈ F \ F Γ defined by the operator

I σ : X D → R such that I σ (v) = s∈Vσ β σ,s v s .
The tetrahedral sub-mesh is defined by

T = {T K,σ,e , e ∈ E σ , σ ∈ F K , K ∈ M}, (8)
where T K,σ,e is the tetrahedron joining the cell centre x K to the triangle T σ,e . For a given v D ∈ X D , we define the function π T v D as the continuous piecewise affine function on each tetrahedron of T such that π

T v D (x K ) = v K , π T v D (x s ) = v s , π T v D (x σ ) = v σ , and π T v D (x σ ) = I σ (v) for all K ∈ M, s ∈ V, σ ∈ F Γ , and σ ∈ F \ F Γ .
The nodal basis of this finite element discretization will be denoted by η

K , η s , η σ , for K ∈ M, s ∈ V, σ ∈ F Γ .
The VAG scheme is a control volume scheme in the sense that it results, for each d.o.f. not located at the Dirichlet boundary and each phase, in a volume balance equation. The two main ingredients are therefore the conservative fluxes and the porous volumes. The VAG matrix and fracture fluxes are exhibited in Figure 4. They are derived from the variational formulation on the finite element subspace. For u D ∈ X D , the matrix fluxes

F K,ν (u D ) connect the cell K ∈ M to all the d.o.f. located at the boundary of K, namely ν ∈ Ξ K = V K ∪ (F K ∩ F Γ ).
They are defined by

F K,ν (u D ) = K -Λ m (x)∇π T u D (x) • ∇η ν (x)dx = ν ∈Ξ K T ν,ν K (u K -u ν ),
with the cell transmissibilities

T ν,ν K = K Λ m (x)∇η ν (x) • ∇η ν (x)dx.
The fracture fluxes F σ,s (u D ) connect each fracture face σ ∈ F Γ to its nodes s ∈ V σ and are defined by

F σ,s (u D ) = σ -d f Λ f ∇ τ γπ T u D (x) • ∇ τ γη s (x)dσ(x) = s ∈Vσ T s,s σ (u σ -u s ),
with the fracture face transmissibilities

T s,s σ = σ d f (x)Λ f (x)∇ τ γη s (x) • ∇ τ γη s (x)dσ(x),
where dσ(x) denotes the Lebesgue d -1 dimensional measure on Γ. The porous volumes are obtained by distributing the porous volumes of each cell K ∈ M and fracture face σ ∈ F Γ to the d.o.f. located on their respective boundaries. For each K ∈ M we define a set of non-negative volume fractions (α K,ν ) ν∈Ξ K \V Dir satisfying ν∈Ξ K \V Dir α K,ν ≤ 1, and we set

φ K,ν = α K,ν K φ m (x)dx.
Similarly, for all σ ∈ F Γ we set

φ σ,s = α σ,s σ φ f (x)d f (x)dσ(x),
with the non-negative volume fractions (α σ,s ) s∈Vσ\V Dir satisfying s∈Vσ\V Dir α σ,s ≤ 1. Then, we set for all K ∈ M and σ ∈ F Γ :

φ K = K φ m (x)dx - ν∈Ξ K \V Dir φ K,ν , φ σ = σ φ f (x)d f (x)dσ(x) - s∈Vσ\V Dir φ σ,s .
On practical meshes with cell sizes at mf interfaces much larger than the fracture width, the flexibility in the choice of the weights α K,s and α σ,s is shown in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] to be a crucial asset compared with usual CVFE approaches, allowing to improve significantly the accuracy of the scheme. As exhibited in Figure 5, and in contrast with the usual CVFE approaches, the fracture porous volumes can be defined with no contribution of the matrix porous volume, thus avoiding to enlarge artificially the flow path in the fractures and to slow down the front speed. This is achieved by choosing the volume fractions such that

α K,σ = 0 for all σ ∈ F Γ , K ∈ M σ , α K,s = 0 for all s ∈ V Γ , K ∈ M s .

VAG phase potential Two-Point (TP) upwind formulation

We consider in the following of Section 3, the usual approach (termed f-upwind model) for which a single rock type is assigned to each d.o.f. Quite naturally, the fracture rock type is associated with d.o.f. located on Γ, while the matrix rock type is associated to the remaining d.o.f., that is we set

P c,ν (s) = P c,m (s) if ν ∈ (V Γ ∪ F Γ ), P c,f (s) if ν ∈ (V Γ ∪ F Γ ),
and

M α ν (s) = M α m (s) if ν ∈ (V Γ ∪ F Γ ), M α f (s) if ν ∈ (V Γ ∪ F Γ ), α ∈ {nw, w}.
The set of discrete unknowns is defined by the set of phase pressure u α D ∈ X D and phase saturation s α D ∈ X D for each phase α ∈ {nw, w}.

The "single" phase VAG Darcy fluxes, not including the phase mobility, are defined, for each phase α ∈ {nw, w}, by

F α K,ν (u α D ) = F K,ν (u α D ) + ρ α gF K,ν (z D ), F α σ,s (u α D ) = F σ,s (u α D ) + ρ α gF σ,s (z D ), with z D = (x ν ) ν∈D , and for K ∈ M, σ ∈ F Γ , ν ∈ Ξ K , s ∈ V σ .
They are combined with the usual Two-Point (TP) phase potential upwinding of the mobilities, leading to the following two-phase Darcy VAG fluxes

q α K,ν = M α K (s α K )(F α K,ν (u α D )) + + M α ν (s α ν )(F α K,ν (u α D )) -, q α σ,s = M α σ (s α σ )(F α σ,s (u α D )) + + M α s (s α s )(F α σ,s (u α D )) -.

Let us define the accumulation terms by

           A α K = φ K s α K , K ∈ M, A α σ = φ σ + K∈Mσ φ K,σ s α σ , σ ∈ F Γ , A α s = K∈Ms φ K,s + σ∈F Γ,s φ σ,s s α s , s ∈ V \ V Dir .
Note that neither the accumulation terms A α σ and A α s nor the mf fluxes take into account the discontinuity of the saturations across mf interface. In other terms, the discrete problem does not involve quantities such as P c,m (s f ). An alternative approach is described in Section 3.5.

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

t 0 = 0 < t 1 < • • • < t n-1 < t n • • • < t N = T of the time interval [0, T ].
We denote the time steps by ∆t n = t nt n-1 for all n = 1, • • • , N . The superscript n will be used to denote the unknowns at time t n . To reduce the amount of notation, only the previous time step superscript n -1 will be specified in the following, while the superscript n will not be specified by default.

The set of discrete equations couples the volume balance equations at each d.o.f. excluding the Dirichlet nodes

                   A α K -A α,n-1 K ∆t n + ν∈Ξ K q α K,ν = 0, K ∈ M, α = nw, w, A α σ -A α,n-1 σ ∆t n + s∈Vσ q α σ,s - K∈Mσ q α K,σ = 0, σ ∈ F Γ , α = nw, w, A α s -A α,n-1 s ∆t n + K∈Ms -q α K,s + σ∈F Γ,s -q α σ,s = 0, s ∈ V \ V Dir , α = nw, w, (9)
combined with the closure laws

   s nw ν + s w ν = 1, ν ∈ D, u nw ν -u w ν = P c,ν (s nw ν ), ν ∈ D, (10)
and the Dirichlet boundary conditions

s nw s = s nw Dir,s u nw s = u nw Dir,s , s ∈ V Dir , (11)
for given s nw Dir,s ∈ [0, 1], u nw Dir,s , s ∈ V Dir .

To solve the discrete nonlinear system (9), one first uses the closure equations (10) to eliminate the unknowns s w ν and u w ν for ν ∈ D reducing the system to the primary unknowns u nw ν , s nw ν , ν ∈ D coupled by the set of equations ( 9) and the Dirichlet boundary conditions [START_REF] Berre | Call for participation: Verification benchmarks for single-phase flow in three-dimensional fractured porous media[END_REF]. A Newton's method is used to solve this nonlinear system at each time step of the simulation. At each Newton step, the Jacobian matrix is assembled and the cell unknowns u nw K , s nw K , K ∈ M are eliminated without any fill-in using the linearized cell volume balance equations reducing the system to the node and fracture face primary unknowns only. This elimination results in a huge gain in terms of system size in particular for tetrahedral meshes. The reduced linear system is solved using a Krylov subspace solver preconditioned by a CPR-AMG preconditioner. This preconditioner combines multiplicatively an AMG preconditioner on a pressure block (elliptic part of the system) with a zero fill-in incomplete factorization of the full system. Let us refer to [START_REF] Lacroix | Decoupling preconditioners in the Implicit Parallel Accurate Reservoir Simulator (IPARS)[END_REF][START_REF] Scheichl | Decoupling and block preconditioning for sedimentary basin simulations[END_REF] for its detailed description. In the following numerical experiments, the pressure block is simply obtain as the sum over both phases of the volume balance equations on each fracture face and non Dirichlet node.

What is wrong with Two-Point upwinding at mf interfaces

In this Section, we discuss one particular difficulty that the nodal discretizations have in regard of the discrete fluxes reconstruction. As shown below, due to the dual control volumes at mf interfaces, nodal schemes may result in fluxes having an opposite sign compared to the fluxes computed at the physical mf interfaces. Using Two-Point upwinding, this results in an artificial diffusion of the saturation toward an upstream direction. To avoid this drawback we propose below an alternative Multi-Point upwinding technique.

For a given constant velocity q, let us choose u D ∈ X D such that u ν = -Λ -1 m q • x ν for all ν ∈ D. From -Λ m ∇π T u D = q, we obtain

F K,s (u D ) = q • K ∇η s (x)dx,
at a given fracture node s ∈ V Γ , and cell K ∈ M s .

Figure 6: Example of a 2D mesh with three isosceles triangular cells at the interface with a fracture in bold. It is assumed that the unit normal vectors are such that n 2 = -n, n 1 = -n 3 and that n 1 • n = 0. The cell centers are chosen as the isobarycenters of their 3 nodes.

For the sake of simplicity, let us assume the geometrical configuration illustrated in Figure 6. It results that

F K,s (u D ) = -1 2 |s 2 σ 1 |q • n 1 , F J,s (u D ) = -1 2 |s 3 σ 2 |q • n 3 = -F K,s (u D ), F I,s (u D ) = + 1 2 |s 2 s 3 |q • n, F K,σ 1 (u D ) = + 1 2 |ss 1 |q • n, F J,σ 2 (u D ) = + 1 2 |ss 4 |q • n = F I,σ 1 (u D ).
We remark that, whatever the velocity q, either the flux F K,s (u D ) or the flux F J,s (u D ) have the opposite sign as the one of q • n. Assuming, to fix ideas that q • n > 0, it results from the Two-Point upwinding used for the transport scheme of a given phase with Darcy velocity q, that the phase propagates from the fracture either to the upstream cell K or the upstream cell J.

On the other hand, let us remark that the ill-orientated discrete fluxes cancel out when summing over the cells connected to the node s and located on the same side with respect to the planar fracture, that is we have

F K,s (u D ) + F I,s (u D ) + F J,s (u D ) = 1 2 |σ 1 σ 2 |q • n. (12)
This property actually holds for an arbitrary number of polygonal cells sharing the node s and whatever the choice of the cell centers. In the three-dimensional case, this property also holds for tetrahedral meshes.

In the following Subsection, this property on the sum of the fluxes is exploited to avoid the artificial diffusion of the phase toward an upstream direction.

Multi-Point (MP) upwind fluxes at mf interfaces

We first define an equivalence relation on each subset M s of cells, for any fixed node s ∈ V, by

K ≡ Ms L ⇐⇒ there exists n ∈ N and a sequence (σ i ) i=1,...,n in F s \F Γ , such that K ∈ M σ 1 , L ∈ M σn and M σ i+1 ∩ M σ i = ∅ for i = 1, . . . , n -1.
Let us then denote by M s the set of all classes of equivalence of M s and by K s the element of M s containing K ∈ M. Obviously M s might have more than one element only if s ∈ V Γ . Note that K s is both considered as a subset of cells of M s as well as an additional d.o.f. located at the same point than the node s i.e. we set x Ks = x s . Let us define the phase mobilities

M α Ks , s ∈ V Γ , K s ∈ M s , α ∈ {nw, w}, (13)
as matrix fracture additional unknowns. Then, using phase potential upwinding of the mobilities, let us define for all L ∈ K s the half Darcy fluxes between L and K s by

q α L,Ks = M α L (s α L )(F α L,s (u α D )) + + M α Ks (F α L,s (u α D )) -,
as well as the half Darcy flux between K s and s by

q α Ks,s = M α s (s α s )(F α Ks,s ) -+ M α Ks (F α Ks,s ) + ,
where we set by flux conservation

F α Ks,s = L∈Ks F α L,s (u α D ).
The flux continuity equation

q α Ks,s = L∈Ks q α L,Ks , (14)
is used to eliminate the mobility unknown M α Ks leading to the following convex linear combination of the cells L ∈ K s and node s mobilities:

M α Ks = L∈Ks (F α L,s (u α D )) + M α L (s α L ) -M α s (s α s )( L∈Ks F α L,s (u α D )) - L∈Ks (F α L,s (u α D )) + -( L∈Ks F α L,s (u α D )) - . ( 15
)
We deduce the definition of the new Multi-Point upwind flux

q α K,s = q α K,Ks = M α K (s α K )(F α K,s (u α D )) + + M α Ks (F α K,s (u α D )) -,
denoted by VAG MP in the following and to be used in the conservation equations [START_REF] Andreianov | Approximating the vanishing capillarity limit of two-phase flow in multi-dimensional heterogeneous porous medium[END_REF]. Compared with the Two-Point upwind flux

q α K,s = M α K (s α K )(F α K,s (u α D )) + + M α s (s α s )(F α K,s (u α D )) -,
denoted by VAG TP, the VAG MP flux uses the fracture node saturation s α s only if L∈Ks F α L,s (u α D ) < 0, which, in view of [START_REF] Berre | Verification benchmarks for single-phase flow in three-dimensional fractured porous media[END_REF], ensures that the phase will not go out from the fracture on the wrong side in the case of a linear phase pressure field.

Note also that, if K s contains only one cell, both the VAG TP and VAG MP fluxes match, this is why the fluxes q α K,σ , σ ∈ F K ∩ F Γ do not need to be modified.

The matrix fracture mobility unknowns [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] and flux continuity equations ( 14) can be kept in the nonlinear system and solve simultaneously with the other unknowns and equations. Let us recall that the CPR-AMG preconditioner combines multiplicatively an AMG preconditioner on a pressure block (elliptic part of the system) with a zero fill-in incomplete factorization of the full system. The matrix fracture mobility unknowns M α Ks and the flux continuity equations ( 14), s ∈ V Γ , K s ∈ M s , are not included in the definition of the pressure block due to their hyperbolic nature. It results that the pressure block has the same number of unknowns and sparsity pattern as the one of the usual VAG TP scheme. Since the AMG step is the most expensive part of the CPR-AMG two stage preconditioner, this explains why keeping the matrix fracture mobility unknowns is quite efficient.

On the other hand, the elimination of the matrix fracture mobility unknowns together with the flux continuity equations in [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces[END_REF] leads to a rather large fill-in of the Jacobian (depending on the density of the fracture network) and also prevents the elimination of the cell unknowns connected to the fractures. The following numerical experiments confirm that it is much more efficient in terms of CPU time to keep the matrix fracture mobility unknowns in the linear system.

Numerical experiments

The Objectives of this Subsection is to compare the solutions obtained with the following schemes:

• the CVFE like VAG scheme with rock type mixture and Two-Point upwinding of the mobilities at mf interfaces (VAG CVFE),

• the VAG scheme with no rock type mixture and Two-Point upwinding of the mobilities at mf interfaces (VAG TP),

• the VAG scheme with no rock type mixture and Multi-Point upwinding of the mobilities at mf interfaces (VAG MP). The VAG MP scheme is implemented either with elimination of the interface mobility unknowns (VAG MP) or without elimination of these unknowns (VAG MP no elim).

• the Hybrid Finite Volume (HFV) scheme with cell, face and fracture edge unknowns as described in [START_REF] Groza | Modelization and discretization of two-phase flows in porous media with discrete fracture networks[END_REF] (HFV).

All these schemes are implemented in the same code using the Fortran 90 programming language combined with the gfortran compiler. The linear systems are solved using the Slatec library [START_REF]SLATEC Common Mathematical Library[END_REF] for the GMRes iterative solver and the ILU0 preconditioner as well as the AMG1R5 library for the Algebraic MultiGrid preconditioner [START_REF] Ruge | Algebraic multigrid (AMG)[END_REF].

Tables 1 and 2 exhibit the following entries:

• mesh: number of cells,

• dof : number of degrees of freedom of each scheme (with 2 physical primary unknowns per d.o.f.),

• dof lin : number of degrees of freedom in the linear system after reduction. Let us recall that the cell unknowns are eliminated for VAG CVFE, VAG TP, HFV, and VAG MP no elim, while the interface mobilities together with the cell unknowns not connected to the fractures are eliminated for VAG MP,

• N z : number of nonzero elements in the reduced Jacobian (with 2 × 2 matrix elements).

Note that for the VAG MP no elim implementation, the pressure block is stored seperately after reduction with a lower number of d.o.f. and nonzero elements than the remaining part of the Jacobian. This is a key point to lower the CPU time when the CPR-AMG preconditioner is used. Then, the first dof lin (resp. N z ) entry corresponds to the pressure block, and the second entry to the remaining part.

These tables also include the following entries:

• N ∆t : number of successful time steps,

• N chop : number of time step chops,

• N N ewton : average number of Newton iterations per successful time step,

• N GM Res : average number of GMRes iteration per Newton step ,

• CPU (s): CPU time in seconds.

The CPU time takes into account the full time loop including the outputs in ensight format files at each time step but excluding the preprocessing computations (mesh reading, mesh connectivity, VAG transmissibilities, CSR format of the Jacobian) which are negligeable in terms of CPU time compared with the time loop.

Tracer DFM model with a single fracture

Let us denote by (x, y) the Cartesian coordinates of x and let us set Ω = (0, 1 m) 2 , x 1 = (0, 1 4 ), x 2 = (1, 0.875). We consider a single fracture defined by Γ = (x 1 , x 2 ) with tangential permeability Λ f = 200 m 2 and width d f = 10 -3 m. The matrix permeability is isotropic and set to Λ m = 1 m 2 . The matrix and fracture porosities are set to φ m = φ f = 1. Let us set t = 1 0.625 , q = 1 1 3

.

We consider the hybrid-dimensional tracer model obtained from the two-phase DFM model by setting M α m (s) = M α f (s) = s for α ∈ {nw, w}, P c,m (s) = P c,f (s) = 0, g = 0. The pressure analytical solution is defined for α ∈ {nw, w} by

u α (x, t) = 1 -Λ -1 m x • q,
leading to the matrix Darcy velocity

q α m = q,
and the tangential fracture velocity integrated over the width

q α f = d f Λ f (t • Λ -1 m q) |t| 2 t,
This pressure solution is exactly solved by the VAG scheme using Dirichlet condition at the boundary of the domain. An input Dirichlet boundary condition is imposed for the non-wetting phase saturation (tracer) with zero value at the matrix boundary and a value of 1 at the fracture boundary x 1 . The initial condition is defined by a zero non-wetting phase saturation both in the fracture and matrix domains. Figure 8 illustrates that the tracer VAG TP solution goes out on the wrong side of the fracture on a few layers of cells, while it is not the case for the HFV and VAG MP solutions as expected. The VAG CVFE stationary tracer solution is not plotted since it is the same than the VAG TP stationary tracer solution. Figure 9 exhibits the stationary solutions along the fracture showing that the HFV and VAG MP solutions match on both meshes while the VAG TP solution is not fully converged even on the fine mesh. Figure 10 exhibits the tracer volume in the fracture as a function of time. Again, the HFV and VAG MP solutions match on both meshes, while the tracer front in the fracture is clearly slown down for the VAG TP solution on both meshes. This is much worse for the VAG CVFE solution due to the fracture enlargement resulting from the rock type mixture at mf interfaces.

Large 2D DFM model

This test case considers the DFM model with the matrix domain Ω = (0, 100 m) × (0, 186.5 m) and a fracture network including 581 connected components both exhibited in Figure 11. The fracture width is d f = 1 cm and the fracture network is homogeneous and isotropic with Λ f = 10 -11 m 2 , φ f = 0.2. The matrix is homogeneous and isotropic with Λ m = 10 -14 m 2 , φ m = 0.4. The relative permeabilities are given by k α r,f (s α ) = s α and k α r,m (s α ) = (s α ) 2 , α ∈ {nw, w} and the capillary pressure is fixed to P c,m (s nw ) = -10 4 ln(1s nw ) Pa in the matrix and to P c,f (s nw ) = 0 Pa in the fracture network. The fluid properties are defined by their dynamic viscosities µ nw = 5. 10 -3 , µ w = 10 -3 Pa.s and their mass densities ρ w = 1000 and ρ nw = 700 Kg.m -3 .

The reservoir is initially saturated with the wetting phase. Dirichlet boundary conditions are imposed at the top boundary with a wetting phase pressure of 1 MPa and s w m = 1, as well as at the bottom boundary with s nw m = 0.9 and u w = 4 MPa. The remaining boundaries are assumed impervious and the final simulation time is fixed to t f = 1800 days.

The time stepping is defined by ∆t 1 = ∆t init = 10 days, and for all n ≥ 1 by ∆t n+1 = max(∆t max , 1.2∆t n ) with ∆t max = 10days, [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF] 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 N ewton max = 30 iterations. This last value is chosen not to small to avoid too many time step failures even on the finest mesh but also not to large to avoid increased CPU time in case of time step failures induced by residual oscillations.

The criterion of convergence for the Newton algorithm is based on a relative residual in l 1 norm smaller than Res max or on a Newton step in l ∞ norm (scaled by 10 -6 for the primary pressure unknown and by 1 for the other primary unknown) smaller than dx max with

Res max = 10 -5 , dx max = 10 -4 .

(

) 17
Note also that the Newton step is relaxed such that its l ∞ norm (scaled by 10 -6 for the primary pressure unknown and by 1 for the other primary unknown) is smaller than dx obj with dx obj = 1.

The non-wetting phase saturation is exhibited at final simulation time in Figure 12 in the matrix and in the fracture network, and the volume of the non-wetting phase as a function of time is presented in Figure 13. We clearly see in Figures 12 and13 that the VAG CVFE discretization considerably slows down the non-wetting phase front in the fracture network due to the drain enlargement induced by the mixing of matrix and fracture porous volumes at mf interfaces. The VAG TP discretization does a better job but still underestimates the front speed in the fracture network. As clearly exhibited by Figure 12, this is due to the fact that the VAG TP scheme propagates the non-wetting phase on the wrong side of the fractures as explained in Subsection 3.2. From Figure 13, the VAG TP solution gets very close to the VAG MP solution after two level of refinement of the coarse mesh, while the VAG CVFE solution has not yet converged on the finest mesh. The comparison between the VAG MP and HFV solutions shows that they are in good agreement for all meshes. It appears in Figure 13 that the HFV scheme converges more slowly then the VAG MP scheme.

The numerical behavior of the four schemes is reported in Table 1 with CPU time is in seconds on Intel E5-2670 2.6 GHz. We remark that the average number of Newton iterations is in all cases quite smaller than N ewton max due to significant variations in the number of Newton iterations during the simulation. This can be explained typically by a higher number of Newton iterations when the non-wetting phase reaches the tips of the fracture network.

For this large 2D network, the VAG MP implementation with elimination of the matrix fracture mobilities leads to a twice large CPU time than the VAG MP implementation with no elimination. Regarding the comparison between VAG MP and VAG TP, we notice a twice larger CPU time, which is a rather good result for such a large network. The comparison between HFV and VAG MP shows for this 2D test case that HFV is competitive on the coarse mesh due to the additional matrix fracture unknowns for VAG MP, but becomes more expensive on the two refined meshes. We will see in the next test case that the situation is much more in favor of the VAG schemes on tetrahedral 3D meshes. for the description of the entries.

3D DFM model

The DFM model of matrix domain Ω = (0, 100 m) 3 and its coarsest tetrahedral mesh conforming to the fracture network are illustrated in Figure 14. The fracture network is assumed to be of constant aperture d f = 1 cm. The matrix and fracture porosities, permeabilities, relative permeabilities and capillary pressures are the same as in the previous test case. The fluid properties are also the same than in the previous test case. At initial time, the reservoir is fully saturated with the wetting phase. Then, nonwetting phase is injected from below, which is managed by imposing Dirichlet conditions at the bottom and at the top of the reservoir. We impose at the bottom boundary either an overpressure ∆p = 2 MPa or no overpressure ∆p = 0 MPa w.r.t. the hydrostatic distribution of the water pressure. The remaining boundaries are assumed impervious and the final simulation time is fixed to t f = 360 days for ∆p = 2 MPa and to t f = 3600 days for ∆p = 0 MPa. The time stepping is defined as in ( 16) using ∆t init = 0.1 days, N ewton max = 30, and either ∆t max = 10 days for ∆p = 2 MPa or ∆t max = 100 days for ∆p = 0 MPa. The criterion of convergence for the Newton algorithm is defined as in [START_REF] Brenner | Gradient discretization of Hybrid-Dimensional Darcy Flows in Fractured Porous Media[END_REF] with Res max = 10 -6 and dx max = 10 -5 , and the relaxation of the Newton step is controlled as in ( 18) by the parameter dx obj = 1.

From Figures 15 and16, we observe that the VAG TP and VAG CVFE schemes are far from convergence even on the finest mesh with 450k cells while the solution provided by the VAG MP scheme is quite close to the one of the HFV scheme. The discrepancy between, on the one hand, the VAG TP and VAG CVFE, and, on the other hand, the VAG MP and HFV schemes is even more striking on the coarse mesh for the no-overpressure gravity dominant test case exhibited in Figures 17 and18. In terms of CPU time, as exhibited in Table 2, the VAG MP scheme implemented with no elimination of the matrix fracture mobility unknowns is competitive compared with the VAG TP scheme. It is also much cheaper than the HFV scheme which leads to a much larger number of d.o.f. and requires both more Newton and GMRes iterations than the VAG schemes. Note that the HFV scheme cannot be run in a reasonable CPU time for the finest mesh of size 1600k cells. Table 2: Numerical behavior of the simulation for the 3D DFM test case with the overpressure ∆p = 2 MPa on the three meshes of sizes 47k, 450k and 1600k cells. The VAG MP scheme is implemented either with elimination (VAG MP) or without elimination (VAG MP no elim) of the mf interface mobility unknowns. We refer to the beginning of Subsection 3.4 for the description of the entries.

Capturing the Saturation jumps at mf interfaces

Given cellwise and fracture facewise constant rock types, the idea introduced in [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for multiphase flow in fractured porous media[END_REF][START_REF] Monteagudu | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF] for CVFE methods and in [START_REF] Eymard | Gradient schemes for twophase flow in heterogeneous porous media and Richards equation[END_REF][START_REF] Brenner | Gradient discretization of Hybrid-Dimensional Darcy Flows in Fractured Porous Media[END_REF][START_REF] Brenner | Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for the VAG scheme is to define as many saturations as rock types shared at a given node or fracture face. This allows to capture the saturation jumps at rock type interfaces resulting from the continuity of the capillary pressure in the graphical sense [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].

The choice of the primary unknowns may greatly affect the convergence of Newton's method used to solve the nonlinear system at each time step of the simulation. For the cells and the nodal d.o.f. associated with a single rock the choice of the primary unknowns does not change compared to Section 3.1. That is we use the non-wetting phase's pressure and saturation as pair of primary unknowns. In contrast the d.o.f. located at rock type interfaces require a special treatment. For such d.o.f. ν ∈ V Γ ∪F Γ we set again the pressure of the non-wetting phase as the first primary unknown, while the second primary unknown is chosen based on the variable switching strategy introduced in [START_REF] Brenner | Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF]. For a given rock type rt ∈ RT = {m, f } let P c,rt denote the monotone graph extension of P c,rt as introduced in [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]. For each subset χ ∈ {{m}, {m, f }} of RT , non-decreasing continuous functions

P c,χ (τ ), S nw χ,rt (τ ), for all rt ∈ χ, (19)
are built such that P c,χ (τ ) ∈ P c,rt (S nw χ,rt (τ )), for all τ and rt ∈ χ, and such that P c,χ (τ ) + rt∈χ S nw χ,rt (τ ) is strictly increasing. Then, we set S w χ,rt (τ ) = 1 -S nw χ,rt (τ ).

The variable τ is going to be used as the second primary unknown.

The main advantage of this framework, which applies to an arbitrary number of rock types, is to incorporate in the construction of the functions [START_REF] Brenner | Positivity-preserving Vertex Approximate Gradient Discretization of two-phase Darcy flows in heterogeneous porous media[END_REF] the saturation jump condition at different rock type interfaces and to apply to general capillary pressure functions. In practice, we use τ = s nw for χ = {m} and the parametrization defined in [START_REF] Brenner | Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for χ = {m, f }. This parametrization is based on a generalization of variable switch approaches (see also [START_REF] Monteagudu | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF]) between s nw f , s nw m , p c and applies to general, including non invertible, capillary functions (see numerical section for an example).

Let us set

rt K = m, K ∈ M, rt σ = f, σ ∈ F Γ , χ ν = {m}, ν ∈ M ∪ (V \ V Γ ), χ ν = {m, f }, ν ∈ V Γ ∪ F Γ .
Using the above framework, given the primary unknowns u nw D = (u nw ν ) ν∈D and τ D = (τ ν ) ν∈D , we set u w D = (u w ν ) ν∈D with u w ν = u nw ν -P c,χν (τ ν ) for all d.o.f. ν ∈ D, and we define the discrete values of the saturation as follows. For all cells K ∈ M and the nodes s ∈ V \ V Γ associated with the single matrix rock type, we set

s α K = S α χ K ,rt K (τ K ) = S α {m},m (τ K ) = τ K s α K,s = S α χs,rt K (τ s ) = S α {m},m (τ s ) = τ s , K ∈ M s .
For the fracture faces σ ∈ F Γ , we set

s α σ = S α χσ,rtσ (τ σ ) = S α {m,f },f (τ σ ) s α K,σ = S α χσ,rt K (τ σ ) = S α {m,f },m (τ σ ), K ∈ M σ .
For the nodes s ∈ V Γ , located at the mf interface, we set As exhibited in Figure 19, the above definition of the saturations at the mf interfaces takes into account the jump of the saturations induced by the different rock types. Let us remark that, in our specific example, since the matrix domain is homogeneous in terms of capillary pressure-saturation relation, the variables s α K,s , K ∈ M s (resp.

Numerical experiments

In this subsection, we compare the m-upwind version of the VAG TP and VAG MP schemes using the same code implementation as described in the beginning of Subsection 3.4. The test case considers the large DFM model exhibited in Figure 21 with domain Ω = (0, 85)×(0, 60)×(0, 140) m kindly provided by the authors of the Benchmark [START_REF] Berre | Call for participation: Verification benchmarks for single-phase flow in three-dimensional fractured porous media[END_REF][START_REF] Berre | Verification benchmarks for single-phase flow in three-dimensional fractured porous media[END_REF]. The fracture width is d f = 1 cm and the fracture network is homogeneous and isotropic with Λ f = 10 -11 m 2 , φ f = 0.2. The matrix is homogeneous and isotropic with Λ m = 10 -14 m 2 , φ m = 0.4.

The relative permeabilities are given by k The parametrization τ at mf interfaces introduced in [START_REF] Brenner | Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] is recalled below and illustrated in Figure 20 for the convenience of the reader.

S nw {m,f },f (τ ) = τ, τ ∈ [0, τ 1 ), 1 -(τ 1 + (1 -τ 1 ) b f bm -τ ) bm b f , τ ∈ [τ 1 , τ 2 ), (21)
S nw {m,f },m (τ ) = 1 -(1 -τ ) b f bm , τ ∈ [0, τ 1 ), τ -τ 1 + 1 -(1 -τ 1 ) b f bm , τ ∈ [τ 1 , τ 2 ), (22)
and

P c,{m,f } (τ ) = -b f ln(1 -τ ), τ ∈ [0, τ 1 ), -b m ln τ 1 + (1 -τ 1 ) b f bm -τ , τ ∈ [τ 1 , τ 2 ), (23)
where The reservoir is initially saturated with the wetting phase. Output Dirichlet boundary conditions are imposed at the boundary {0, 85} × (0, 20) × (110, 140) with a wetting phase pressure of 1 MPa and s w m = 1, and input Dirichlet boundary conditions are set at the boundary {0} × (40, 60) × (0, 30) ∪ (0, 30) × (40, 85) × {0} with s nw m = 0.9 and u w = 4 MPa. The remaining boundaries are assumed impervious and the final simulation time is fixed to t f = 3600 days. The time stepping is defined as in ( 16) using ∆t init = 0.01 days, ∆t max = 100 days and N ewton max = 25. The criterion of convergence for the Newton algorithm is defined as in [START_REF] Brenner | Gradient discretization of Hybrid-Dimensional Darcy Flows in Fractured Porous Media[END_REF] with Res max = 10 -5 and dx max = 10 -4 , and the relaxation of the Newton step is controlled as in [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF] by the parameter dx obj = 1. Table 3: Numerical behavior of the simulation Field test case for the VAG TP and MP m-upwind schemes. We refer to the beginning of Subsection 3.4 for the description of the entries.

τ 1 = 1 -( b f bm ) bm bm-b f and τ 2 = τ 1 + (1 -τ 1 ) b f bm .
The same issue at mf interfaces as for the VAG TP f-upwind approximation can be noticed in Figure 23 for the VAG TP m-upwind discretization in the sense that the non-wetting phase can go out from the fractures on the wrong side for the VAG TP approximation. Nevertheless, thanks to the rather large saturation jump captured by the m-upwind model in this test case, it involves small amounts of the non-wetting phase and does not have visible effects on overall quantities (see Figure 22) nor on the non-wetting phase saturation front (see Figure 23). In terms of CPU time, as exhibited in Table 3, a factor of roughly 1.7 is observed in favor of the TP discretization due to the additional mf interface unknowns on this rather large fracture network and to the slightly larger number of Newton iterations for the MP scheme.

Let us refer to [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid-Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] for a numerical comparison between the m-upwind VAG TP scheme and the m-upwind VAG CVFE scheme (i.e. without adaptive distribution of the porous volumes at mf interfaces). It shows that the m-upwind VAG CVFE scheme still slows down the transport in the fractures in particular for a high matrix fracture permeability ratio.

4 Vertex Approximate Gradient (VAG) discretization of two-phase DFM discontinuous pressure models

Discontinuous pressure models are required to account for fractures acting as barriers. Such barriers are usually induced by a low fracture normal permeability combined with a capillary barrier effect. Note that even in the case of a high normal fracture permeability, a barrier behavior can still be observed for a given phase due to the degeneracy of the phase mobility when the fracture is filled by the other phase (see [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF][START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF]). Compared to the single phase flow models the possibility of such capillary barriers constitutes an additional motivation for the use of discontinuous pressure models. VAG discrete unknowns: as exhibited in Figure 24, the discrete unknowns are defined by the matrix d.o.f.

D m = M ∪ {K s | K s ∈ M s , s ∈ V \ V Γ } ∪ D mf
and by the fracture d.o.f.

D f = F Γ ∪ V Γ ,
where D mf ⊂ D m are the mf interface d.o.f.

D mf = M Γ ∪ F Γ , with M Γ = {K s | K s ∈ M s , s ∈ V Γ }, F Γ = {K σ | K ∈ M σ , σ ∈ F Γ }.
Let us set D = D m ∪ D f and let us remark that for s ∈ V \ V Γ , M s is reduced to the set of cells around s and the d.o.f. K s ∈ M s is considered to match with the node s.

For each cell K ∈ M, let us also define the following subset of d.o.f. located at the boundary of the cell:

Ξ K = {K s , s ∈ V K , K σ , σ ∈ F Γ ∩ F K }.
The subset of Dirichlet d.o.f. is denoted by D Dir ⊂ D.

As in Section 3.5, the definition of the primary and secondary unknowns at the d.o.f. located at the rock type interfaces is based on the parametrization of the capillary pressure graphs [START_REF] Brenner | Positivity-preserving Vertex Approximate Gradient Discretization of two-phase Darcy flows in heterogeneous porous media[END_REF]. To fix ideas, we assume the presence of 3 rock types RT = {m, f d , f b } where f d is a fracture drain rock type and f b is a fracture barrier rock type while m denote again the matrix rock type Let the fracture network Γ be partitioned into the networks Γ d of fractures acting as drains and the network Γ b of fractures acting as barriers. In order to simplify the presentation of the numerical scheme, we will assume that Γ d ∩ Γ b = ∅. Then, the collection χ of rock types associated with any given d.o.f. take values in

{{m}, {f d }, {f b }, {m, f d }, {m, f b }}.
corresponding to assume no intersections between fractures acting as drain and barrier. In practice, we use the parametrization τ = s nw for χ = {m}, {f d }, {f b } and the parametrizations defined in [START_REF] Brenner | Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] for χ = {m, f d } or {m, f b }. More precisely, let us set for = b, d

rt K = m, K ∈ M, rt σ = f , {σ ∈ F Γ | x σ ∈ Γ },    χ ν = {m}, ν ∈ D m \ D mf , χ ν = {f }, {ν ∈ D f | x ν ∈ Γ } := D f , χ ν = {m, f }, {ν ∈ D mf | x ν ∈ Γ } := D mf ,
Using the above framework, given the primary unknowns u nw D = (u nw ν ) ν∈D and τ D = (τ ν ) ν∈D , we set u w D = (u w ν ) ν∈D with u w ν = u nw ν -P c,χν (τ ν ) for all d.o.f. ν ∈ D, and we define the discrete values of the saturation as follows. For all d.o.f. associated with a single rock type, that is K ∈ M and σ ∈ F Γ we set

s α K = S α χ K ,rt K (τ K ) = τ K , s α σ = S α χσ,rtσ (τ σ ) = τ K , for all ν ∈ Ξ K ∩ D m \ D mf , K ∈ M, we set s α K,ν = S α χν ,rt K (τ ν ) = S α {m},m (τ ν ) = τ ν ,
for all s ∈ V s , σ ∈ F Γ , we set s α σ,s = S α χs,rtσ (τ s ) = τ s , while for all mf interface d.o.f. from D mf , and with = b, d, we impose

     s α K,ν = S α χν ,rt K (τ ν ) = S α {m,f },m (τ ν ), ν ∈ Ξ K ∩ D mf , K ∈ M, s α σ,Kσ = S α χ Kσ ,rtσ (τ Kσ ) = S α {m,f },f d (τ Kσ ), K σ ∈ F Γ ∩ D mf , s α σ,Ks = S α χ Ks,rtσ (τ Ks ) = S α χ {m,f },f d (τ Ks ), K s ∈ M Γ ∩ D mf , σ ∈ F Γ,Ks , where F Γ,Ks = F Γ,s ∩ ( K∈Ks F K ).
Discrete fluxes: the VAG fluxes connect each cell K (resp. each fracture face σ) to its boundary d.o.f. ν ∈ Ξ K (resp. s ∈ V σ ) using the same transmissibility coefficients as for the continuous pressure model

F K,ν (u Dm ) = ν ∈Ξ K T ν,ν K (u K -u ν ), F σ,s (u D f ) = s ∈Vσ T s,s σ (u σ -u s ).
Additionally, two-point matrix fracture fluxes are defined by The tests are driven on triangular meshes, extended to 3D prismatic meshes by adding a second layer of nodes as a translation of the original nodes in normal direction to the plane of the original 2D domain (cf. Figure 4.1). The equi-dimensional mesh contains two layers of cells in the fractures. In order to focus on modelling errors, the meshes are chosen to be fine with cell sizes of the same order as the fracture aperture. The final simulation time is fixed to t f = 54 days. The time stepping is defined as in (16) using ∆t init = 0.01 days and ∆t max = 0.1 days for the equi-dimensional and hybrid dimensional mf nonlinear models, and ∆t init = 0.002 days and ∆t max = 0.27 days for the hybrid-dimensional mupwind and f-upwind models. The maximum number of Newton iterations per time step is fixed as N ewton max = 35. The criterion of convergence for the Newton algorithm is defined as in [START_REF] Brenner | Gradient discretization of Hybrid-Dimensional Darcy Flows in Fractured Porous Media[END_REF] with Res max = 10 -6 and dx max = 10 -4 , and the relaxation of the Newton step is controlled as in [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF] by the parameter dx obj = 0.5.

F Ks,s (u Ks , u s ) = T Ks,s (u Ks -u s ), F Kσ,σ (u Kσ , u σ ) = T Kσ,σ (u Kσ -u σ ),
The hybrid dimensional mf nonlinear and m-upwind models make use of the parametrization ( 21)-( 22)-( 23) at the mf interfaces. Table 4: Numerical behavior of the simulation obtained by the VAG scheme for the equi-dimensional model and for the mf nonlinear, m-upwind and f-upwind discontinuous pressure DFM models, as presented in section 4. We refer to the beginning of Subsection 3.4 for the description of the entries with dof lin and N z accounting for the elimination of the cell unknowns in the linear systems.

In this test case, we study the presence of a fracture, which acts as a barrier, both by its low permeability and by its high capillarity compared to the rock matrix. As a result of the higher capillarity, the sign of the matrix-fracture non-wetting phase saturation jump S nw m (γ ± p c,m ) -S nw f (γ ± p c,m ) at the mf interfaces is non negative. Figures 27 and 28 compare the above mf nonlinear, m-upwind and f-upwind discontinuous pressure models to a reference equi-dimensional model. For the f-upwind and m-upwind models, mass transfer of the non-wetting phase from the matrix to the barrier is overestimated, since in this direction, saturation jumps are not accounted for. The assumption of constant saturation across the fracture for these models consequently leads to an overestimation of the non-wetting phase leaving the barrier. This overestimation is most severe for the m-upwind model, which takes into account saturation jumps for fluxes directed from the fracture to the matrix. Again, the mf nonlinear model does not suffer from the difficulties described above, since it provides mass transport that passes by the mf interfaces and takes into account the saturation jumps.

Conclusions and Perspectives

This article reviews the nodal VAG discretization of DFM two-phase Darcy flow models. For linear transmission conditions, the adaptation of the control volumes combined with a Multi-Point upwind approximation of the mobilities for f-upwind models or taking into account the saturation jump for m-upwind models, allows to obtain a similar accuracy as face based discretizations with a much lower CPU time on tetrahedral meshes. Nonlinear mf transmission conditions provide a more accurate DFM model than linear transmission conditions. As discussed in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF][START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF], they can account for a large range of physical processes at mf interfaces which cannot be captured by linear mf transmission conditions even in the case of fractures acting as drains. It is typically the case for fractures acting as capillary barriers, or for highly permeable fractures filled with a given phase acting as a barrier for the other phase. The VAG discretization of DFM models with nonlinear mf transmission conditions still raises the issue of numerical efficiency regarding the nonlinear convergence due to the combination of highly nonlinear transmission conditions with tiny volumes at mf interfaces. Improving the numerical efficiency for this type of DFM models is the object of ongoing researches in two directions. The first is to go back to face based discretizations allowing the elimination of the mf interface unknowns with a local nonlinear interface solver as in [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF] using TPFA discretization on orthogonal meshes and in [START_REF] Aghili | Hybrid Finite Volume discretization of two-phase Discrete Fracture Matrix models with nonlinear interface solver[END_REF] using an HFV discretization. The second perspective is to use the more robust Hybrid Upwinding approximation of the mobilities to define the two-phase Darcy fluxes at mf interfaces as proposed in [START_REF] Alali | Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions[END_REF] for TPFA schemes and in [START_REF] Brenner | Positivity-preserving Vertex Approximate Gradient Discretization of two-phase Darcy flows in heterogeneous porous media[END_REF] for the VAG discretization.

Figure 1 :

 1 Figure 1: Example of a 2D DFM with the matrix domain Ω and 3 intersecting fractures Γ i , i = 1, 2, 3.

Figure 4 :

 4 Figure 4: For a cell K and a fracture face σ (in bold), examples of VAG d.o.f. u K , u s , u σ , u s and VAG fluxes F K,σ , F K,s , F K,s , F σ,s .

Figure 5 :

 5 Figure 5: Example of control volumes at cells, fracture face, and nodes, in the case of two cells K and L splitted by one fracture face σ (the width of the fracture has been enlarged in this Figure). (left): VAG choice of the porous volumes avoiding mixing between fracture and matrix porous volumes. (right): CVFE like choice of the porous volumes mixing fracture and matrix porous volumes leading to a considerable enlargement of the fracture drain on practical meshes.

Figure 7 :

 7 Figure 7: (left) 2D mesh with 3 fracture faces in bold and the 3 d.o.f. in M s at the node s ∈ V Γ , (right) Darcy fluxes joining each cell L ∈ K s to the new d.o.f. K s , and joining the new d.o.f. K s to the node s (the node s and K s are located at the same point s but they have been separated for the sake of clarity of the Figure).

Figure 8 :

 8 Figure 8: Stationary solution for the non-wetting phase saturation (tracer) in the matrix and in the fracture obtained by, from left to right, the HVF, VAG MP and VAG TP schemes, and, from top to bottom, on the 16 × 16 and 128 × 128 topologically Cartesian meshes.

Figure 9 :

 9 Figure 9: Stationary non-wetting phase saturation along the fracture as a function of x obtained by the HVF, VAG MP, VAG TP schemes on the 16 × 16 and 128 × 128 topologically Cartesian meshes.

Figure 10 :

 10 Figure 10: Volume of the non-wetting phase in the fracture as a function of time for the HVF, VAG MP, VAG TP, VAG CVFE scheme solutions on the 16 × 16 and 128 × 128 topologically Cartesian meshes.

Figure 11 :

 11 Figure 11: Triangular mesh of the DFM model with 32340 (32k) cells and 5344 fracture faces (Courtesy of M. Karimi-Fard, Stanford, and A. Lapène, Total). This mesh is refined uniformly to obtain the 129k and 517k cells meshes.

Figure 12 :

 12 Figure 12: Non-wetting phase saturation in the matrix and fracture network at time t f = 1800 days for the HVF, VAG MP, VAG TP, VAG CVFE schemes from left to right, and the 32k, 129k, 517k cells meshes from top to bottom.

Figure 13 :

 13 Figure 13: Non-wetting phase volume in the fracture network as a function of time for the VAG MP, HFV, VAG TP and VAG CFVE schemes on the 3 meshes of sizes 32k, 129k and 517k cells.

Figure 14 :

 14 Figure 14: Geometry of the domain Ω = 100m×100m×100 m with the fracture network in red (left), coarsest tetrahedral mesh with 47670 cells (right).

Figure 15 :

 15 Figure 15: Non-wetting phase saturation solutions obtained with the HVF, VAG MP, VAG TP, VAG CVFE schemes from left to right, at time t f = 360 days (top), and at time t = 100 days (bottom),with overpressure ∆p = 2 MPa, and the mesh of size 450k cells. The threshold in the matrix is S nw m > 0.1 (bottom).

Figure 16 :

 16 Figure 16: Non-wetting phase volume in the fracture network as a function of time for the 3D DFM test case with the overpressure ∆p = 2 MPa using the VAG MP, HFV, VAG TP and VAG CFVE schemes on the 2 meshes of sizes 47k and 450k cells.

Figure 17 :

 17 Figure 17: Non-wetting phase saturation solutions obtained with the HVF, VAG MP, VAG TP, VAG CVFE schemes from left to right, at time t f = 3600 days with no overpressure ∆p = 0 MPa, and the mesh of size 47k cells. The threshold in the matrix is S nw m > 0.1 (bottom).

Figure 18 :

 18 Figure 18: Non-wetting phase volume in the fracture network as a function of time for the 3D DFM test case with no overpressure ∆p = 0 MPa using the VAG MP, HFV, VAG TP and VAG CFVE schemes on the 47k cells mesh.

  s α K,s = S α χs,rt K (τ s ) = S α {m,f },m (τ s ), s α σ,s = S α χs,rtσ (τ s ) = S α {m,f },f (τ s ), σ ∈ F Γ,s .

Figure 19 :

 19 Figure 19: Saturations inside the cells K and L, the fracture face σ and at the mf interfaces taking into account the saturation jumps induced by the different rock types.

  α r,f (s α f ) = s α f and k α r,m (s α m ) = (s α m ) 2 , α ∈ {nw, w} and the capillary pressure is fixed to P c,m (s nw m ) = -b m ln(1-s nw m ) Pa in the matrix and to P c,f (s nw f ) = -b f ln(1s nw f ) Pa in the fracture network, with b f = 10 3 Pa, and b m = 10 4 Pa. The fluid properties are defined by their dynamic viscosities µ nw = 5. 10 -3 , µ w = 10 -3 Pa.s and their mass densities ρ w = 1000 and ρ nw = 700 Kg.m -3 .

Figure 20 :

 20 Figure 20: (Top): capillary pressure as a function of the non-wetting phase saturation for both the fracture (f) and matrix (m) rock types with b m = 10 4 and b f = 10 3 Pa. (Bottom): capillary pressure and fracture and matrix non-wetting phase saturations as functions of the parameter τ ∈ [0, τ 2 ).

Figure 21 :

 21 Figure 21: Large DFM model with its mesh of size 495233 tetrahedral cells and 66908 fracture faces provided by the authors of the Benchmark [11].

Figure 22 :

 22 Figure 22: Non-wetting phase saturation volumes in the matrix (left) and in the fracture network (right) as a function of time obtained for the VAG TP and the VAG MP m-upwind schemes.

Figure 23 :

 23 Figure 23: Non-wetting phase saturation in the matrix (top) and in the fracture network (bottom) obtained for the VAG TP (left) and VAG MP (right) m-upwind schemes at time t = 350 days.

Figure 25 :

 25 Figure 25: Coarse mesh over the domain under consideration, which contains two intersecting fractures with high permeability and low capillarity and one upper fracture with low permeability and high capillarity. The size of the domain is 4m×8m and the fractures have an aperture of 4 cm.

Figure 26 :

 26 Figure 26: Comparison of the equi-dimensional model and of the mf nonlinear, m-upwind and f-upwind discontinuous pressure DFM models (from left to right) numerical solutions for non-wetting phase saturation at final time t = 54 days.

Figure 27 :Figure 28 :

 2728 Figure 27: Comparison of the equi-dimensional model and of the mf nonlinear, m-upwind and f-upwind discontinuous pressure DFM models (from left to right) numerical solutions for water overpressure at final time t = 54 days.

Table 1 :

 1 Numerical behavior of the simulation for the 2D DFM test case on the 3 meshes and for the HFV, VAG CVFE, VAG TP and VAG MP schemes. The VAG MP scheme is implemented either with elimination (VAG MP) or without elimination (VAG MP no elim) of the mf interface mobility unknowns. We refer to the beginning ofSubsection 3.4

	Scheme	mesh	dof	dof lin	N z	N ∆t	N chop	N N ewton	N GM Res	CPU (s)
	HFV	32k	87k	54k	280k	180	0	5.9	16.0	465
	VAG CVFE	32k	54k	21k	161k	180	0	3.0	8.1	165
	VAG TP	32k	54k	21k	161k	180	0	3.7	14.3	224
	VAG MP	32k	65k	43k	459k	188	4	5.2	19.5	794
	VAG MP no elim.	32k	65k	21/32k	161/356k	180	0	4.8	20.4	491
	HFV	129k	335k	205k	1045k	180	0	8.5	56	7747
	VAG CVFE	129k	205k	75k	549k	180	0	4.8	26	1301
	VAG TP	129k	205k	75k	549k	180	0	6.7	55	2995
	VAG MP	129k	226k	122k	1143k	261	37	9.3	54	11880
	VAG MP no elim.	129k	226k	75/96k	549/916k	182	1	8.0	64	5417
	HFV	517k	1315k	798k	4.0M	182	1	12.6	104	145704
	VAG CVFE	517k	798k	280k	2.0M	180	0	9.7	59	25403
	VAG TP	517k	798k	280k	2.0M	180	0	10.9	103	50472
	VAG MP no elim	517k	840k	280/322k	2.0/2.7M	199	12	12.8	122	98390

s α K,σ , K ∈ M σ ) refer to the same nodal (resp. facial) saturation values. Similarly, the values s α σ,s , σ ∈ F Γ,s are identical. This is however not true for general heterogeneous matrix and fracture domains.

We define the accumulation terms by

and the VAG fluxes with TP phase potential upwinding of the mobilities by

The VAG TP discretization capturing the saturation jumps at rock type interfaces looks for u nw D and τ D satisfying the conservation equations ( 9) together with the Dirichlet boundary conditions

It will be termed VAG TP m-upwind discretization in the following. The VAG MP mupwind discretization can also be defined as previously using the MP upwind flux

for s ∈ V Γ with the interface mobility

, assuming that rt L = rt K for all L ∈ K s . This assumption can always be verified by setting new interface face(s) between the different rock types in K s . This discretization will be termed VAG MP m-upwind discretization in the following.

A comparison of the f-upwind and m-upwind models with reference equi-dimensional solutions can be found in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF] and [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF]. Basically, it concludes that, thanks to the saturation jump capturing at mf interfaces, the m-upwind model provides a better approximation than the f-upwind model as long as the fractures are not fully filled with the non-wetting phase. When the fractures are filled, the m-upwind model overestimates the fracture capillary pressure and underestimates the capillary barrier effect. In that case the fupwind model provides a better approximation. In the following numerical section, the VAG TP and MP m-upwind discretizations are compared both in terms of solutions and CPU times.

where ∆ is the triangular submesh of Γ defined as the trace on Γ of the tetrahedral submesh T introduced in (8) (see [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces[END_REF]) for details).

Setting z Dm = (z ν ) ν∈Dm and z D f = (z ν ) ν∈D f , the two-phase VAG fluxes combine the VAG single phase Darcy fluxes including gravity

with the usual Two-Point phase potential upwinding of the mobilities, leading to define

for all s ∈ V Γ , K s ∈ M s , and

Control volumes and accumulation terms: as for the continuous pressure model, porous volumes φ K,ν , ν ∈ Ξ K \ D Dir (resp. φ σ,s , s ∈ V σ \ D Dir ) are obtained by distribution of the the cell K ∈ M (resp. fracture face σ ∈ F Γ ) porous volume. A porous volume φ σ,Ks (resp. φ σ,Kσ ) is also distributed from the fracture face σ to the interface d.o.f. K s (resp. K σ ) for σ ∈ F Γ,Ks , K s ∈ M Γ (resp. for K σ ∈ F Γ ). These interface porous volumes are required to avoid the singularity of the linear systems obtained after Newton linearization. Their influence on the solution is small provided that they are chosen small enough (see [START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF]). Then we set

and we define the accumulations terms by

Conservation equations: the VAG discretization of the discontinuous pressure model solves for u nw D and τ D such that

f and m-upwind discontinuous pressure models: the above discontinuous pressure model, termed mf nonlinear model in the following, leads to difficulties to solve the nonlinear system (24) due to the combination of highly contrasted matrix and fracture rock types and to the small pore volumes at mf interface d.o.f. One possibility to solve this issue, still preserving the ability to take into account fractures acting as drains or barriers, is to linearize the matrix fracture transmission conditions w.r.t. the mf interface unknowns and to apply a f or m-upwind approximation of the mobilities. This idea, developed in [START_REF] Brenner | Hybrid-dimensional modeling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF] for the VAG discretization and in [START_REF] Aghili | Two-phase Discrete Fracture Matrix models with linear and nonlinear transmission conditions[END_REF] for the TPFA discretization, replaces the primary unknowns u nw ν , τ ν at matrix fracture d.o.f. ν ∈ D mf by both phase pressures u nw ν , u w ν , ν ∈ D mf , and the conservation equations at matrix fracture d.o.f. by

Note that the pore volumes φ σ,Ks and φ σ,Kσ are set to zero. Since phase saturations are no longer defined at matrix fracture d.o.f., one need to modify the upwind mobilities in the definition of the fluxes q α K,Ks , K s ∈ M Γ now connecting directly the cell K and the fracture d.o.f. s, and in the definition of q α K,Kσ , K σ ∈ F Γ now connecting the cell K and the fracture d.o.f. σ. These new connectivities modify the fracture conservations equations for σ ∈ F Γ and s ∈ V Γ \ D Dir as follows:

The modified fluxes are defined by

for the m-upwind discontinuous pressure model, and by

for the f-upwind discontinuous pressure model, where a fracture rock type rt s has been assigned to the node s. As for the continuous pressure model, a Multi-Point upwinding can also be introduced for these fluxes using the additional mobility unknowns M α Ks , and M α Kσ , α ∈ {nw, w}. Note that, for fracture acting as drains, these f and m-upwind discontinuous pressure models provide basically the same solutions than respectively the f and m-upwind continuous pressure models. As already mentioned, this is not the case of the mf nonlinear discontinuous pressure model [START_REF] Chave | A Hybrid High-Order Method for Darcy Flows in Fractured Porous Media[END_REF] due to the possible degeneracy of the phase mobilities appearing in the matrix fracture transmission conditions.

Numerical experiments

In this subsection, we compare on the following test case, the mf nonlinear, the m-upwind and the f-upwind models using a reference solution obtained by the equi-dimensional model. The code implementation is the same for all models and described in the beginning of Subsection 3.4. The m-upwind and f-upwind models would require the design of specific preconditioners due to the two independent elliptic pressure unknowns at mf interfaces combined with a single independent elliptic unknown at cells and fracture faces. This explains the use for these two models of the direct linear solver SuperLU from the library [START_REF] Li | SuperLU Users' Guide[END_REF]. The GMRes iterative solver combined with the CPR-AMG preconditioner is still used for the mf nonlinear and equi-dimensional models. It results that the overall numbers of Newton iterations N N ewton N ∆t are more relevant for performance comparison than the CPU times which are not reported for this test case.

We consider a fractured domain as defined in Figure 4.1. The matrix permeability is isotropic of 0.1 Darcy and matrix porosity is 0.2. The two lower fractures are drains (f d ) of isotropic permeability 100.0 Darcy and porosity 0.4. In the upper fracture, acting as a barrier (f b ), the permeability is isotropic of 0.001 Darcy and the porosity is 0.