K Brenner

J Hennicker
email: ju-lian.hennicker@unice.fr

R Masson

P Samier

Hybrid Dimensional Darcy Flow in Fractured Porous Media with discontinuous pressures at the matrix fracture interfaces

 to diphasic flow. We thus provide a model for two phase Darcy flow through fracture networks in porous media, in which the d -1 dimensional flow in the fractures is coupled with the d dimensional flow in the matrix, leading to the so called hybrid dimensional Darcy flow model. It accounts for fractures acting either as drains or as barriers, since it allows pressure jumps at the matrix-fracture interfaces. The model also permits to treat gravity dominated flow as well as discontinuous capillary pressure at the material interfaces. We adapt the Vertex Approximate Gradient (VAG) scheme to this problem, in order to account for anisotropy and heterogeneity aspects as well as for applicability on general meshes. Several test cases are presented to compare our hybrid dimensional model to the hybrid dimensional, continuous pressure model (proposed in [7]) and to the generic equidimensional model, in which fractures have the same dimension as the matrix. This does not only provide quantitative evidence about computational gain, but also leads to deep insight about the quality of the proposed reduced model.

Introduction

This work has two aims: Providing a reduced model for two phase flows in porous media with complex Discrete Fracture Networks (DFN) and validating the reduced model by comparing numerically derived solutions of different test cases with the solutions of the full (non reduced) model. More precisely, we are concerned with the modelling and the discretization of two phase Darcy flows in fractured porous media, for which the fractures are represented as interfaces of codimension one. In this framework, the d -1 dimensional flow in the fractures is coupled with the d dimensional flow in the matrix leading to the so called, hybrid dimensional Darcy flow model. These models are derived from the so called equi-dimensional model, where fractures are represented as geological layers of equal dimension as the matrix, by averaging fracture quantities over the fracture width. We consider the case for which the pressure can be discontinuous at the matrix-fracture (mf ) interfaces in order to account for fractures acting either as drains or as barriers as described in [START_REF] Flauraud | Domain Decomposition for an asymptotic geological fault modeling[END_REF], [START_REF] Jaffré | Modeling fractures and barriers as interfaces for flow in porous media[END_REF], [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF], contrary to the continuous pressure model described in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] developed for highly conductive fractures. A hybrid dimensional discontinuous pressure model for two phase flow in global pressure formulation has been derived in [START_REF] Jaffré | A discrete fracture model for two-phase flow with matrixfracture interaction[END_REF]. The model presented in this work, in pressure-pressure formulation, provides features like an upwind coupling condition for mf mass exchange fluxes and the incorporation of gravitational force in these fluxes, which is a novelty. Subsequently, in this work, we use numerically derived solutions of different test cases, to compare our model with the equi-dimensional model and with the hybrid dimensional model for complex DFN, presented in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF], which assumes pressure continuity accross the fractures.

The discretization of such hybrid dimensional Darcy flow models has been the object of several works. For monophasic Darcy flow, a cell-centered Finite Volume scheme using a Two Point Flux Ap-proximation (TPFA) is proposed in [START_REF] Flauraud | Domain Decomposition for an asymptotic geological fault modeling[END_REF], [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF] assuming the orthogonality of the mesh and isotropic permeability fields. Cell-centered Finite Volume schemes using MultiPoint Flux Approximations (MPFA) have been studied in [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] and [START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF]. In [START_REF] Jaffré | Modeling fractures and barriers as interfaces for flow in porous media[END_REF], a Mixed Finite Element (MFE) method is proposed and a MFE discretization adapted to non-matching fracture and matrix meshes is studied in [START_REF] D'angelo | A mixed finite element method for Darcy flow in fractured porous media with non-matching grids[END_REF]. More recently the Hybrid Finite Volume (HFV) scheme, introduced in [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF], has been extended in [START_REF] Faille | Reduced models for flow in porous media containing faults with discretization using hybrid finite volume schemes[END_REF] for the non matching discretization of two reduced fault models. Also a Mimetic Finite Difference (MFD) scheme is used in [START_REF] Antonietti | Mimetic Finite Difference Approximation of flows in Fractured Porous Media[END_REF] in the matrix domain coupled with a TPFA scheme in the fracture network. Discretizations of the related reduced model [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] assuming a continuous pressure at the matrix fracture interfaces have been proposed in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] using a MFE method and in [START_REF] Brenner | Gradient discretization of Hybrid Dimensional Darcy Flows in Fractured Porous Media[END_REF] using the HFV scheme and an extension of the Vertex Approximate Gradient (VAG) scheme introduced in [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF]. Finally, in [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces[END_REF], the VAG and HFV schemes have been extended to the monophasic counterpart of the model presented in this work. For diphasic Darcy flow, a cell-centered Finite Volume scheme using a TPFA is proposed in [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF]. In [START_REF] Hoteit | Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures[END_REF], a mixed finite element method has been adapted to hybrid dim. two phase flow through fractured pourous media. The hybrid dim. continuous pressure model for diphasic flow is discretized 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] using a Control Volume Finite Element method (CVFE) and in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF] using the VAG scheme. To the author's knowledge, there has not yet appeared a comparison of different hybrid dimensional models with the generic equi-dimensional model for two phase flow. This is one achievement of the present paper.

In this work, we present an adaptation of the VAG scheme to the hybrid dim. discontinuous pressure model, with supplementary unknowns at the mf interfaces, to capture the pressure jumps at the fractures. We choose a vertex based scheme, since it is well adapted for symplectic meshes, which is a necessary feature when dealing with complex geometries. Furthermore, the control volume version of the VAG scheme, presented here, allows to take into account saturation jumps (due to capillary pressure) at heterogeneous unknowns, including the unknowns located at the mf interfaces. This is possible, because fluxes are local to each cell and fracture face, respectively. A third advantage of this method becomes obvious in the test case section of this work. To capture gravitational effects in normal direction within the fractures, the supplementary unknowns at the mf interfaces are needed, which excludes purely cell-centered approaches, where face unknowns are eliminated.

The outline of this work is as follows. The hybrid dimensional two phase flow model is provided in the first section. The second section is devoted to the VAG discretization and provides a finite volume formulation of the model. In the third section, the model is compared to the equidim. model and to the hybrid dim. model with continuous pressure at the mf interfaces, via numerical solutions derived with the VAG scheme for different test cases.

2 Hybrid dimensional Model in Fractured Porous Media

Geometry and Function Spaces

Let Ω denote 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 in the next section. The adaptations to the case d = 2 are straightforward.

Let Γ = i∈I Γ i and its interior Γ = Γ \ ∂Γ denote the network of fractures Γ i ⊂ Ω, i ∈ I, such that each Γ i is a planar polygonal simply connected open domain included in a plane P i of R d . It is assumed that the angles of Γ i are strictly smaller than 2π, and that Γ i ∩ Γ j = ∅ for all i = j . For all i ∈ I, let us set Σ i = ∂Γ i , with n Σ i as unit vector in P i , normal to Σ i and outward to Γ i . Further

Σ i,j = Σ i ∩Σ j , j ∈ I \{i}, Σ i,0 = Σ i ∩∂Ω, Σ i,N = Σ i \( j∈I\{i} Σ i,j ∪Σ i,0 ), Σ = (i,j)∈I×I,i =j (Σ i,j \Σ i,0 ) and Σ 0 = i∈I Σ i,0 . It is assumed that Σ i,0 = Γ i ∩ ∂Ω. Γ 2 Γ 3 Γ 1 Σ 1,0 Σ 2,0 Σ Ω Σ 3,N n a + (2) n a -(2) Ω n a + (3) Γ n a + (1)
Figure 1: Example of a 2D domain Ω and 3 intersecting fractures Γ i , i = 1, 2, 3. We define the fracture plane orientations by a ± (i) ∈ χ for Γ i , i ∈ I.

denote by H 1 Σ 0 (Γ). On Ω\Γ, the gradient operator from H 1 (Ω\Γ) to L 2 (Ω) d is denoted by ∇. On the fracture network Γ, the tangential gradient, acting from H 1 (Γ) to L 2 (Γ) d-1 , is denoted by ∇ τ , and such that

∇ τ v = (∇ τ i v i ) i∈I ,
where, for each i ∈ I, the tangential gradient

∇ τ i is defined from H 1 (Γ i ) to L 2 (Γ i ) d-1
by fixing a reference Cartesian coordinate system of the plane P i containing Γ i . We also denote by div τ i the divergence operator from H div (Γ i ) to L 2 (Γ i ) and by γ τ i the tangential trace operator from

H 1 (Ω) d to L 2 (Γ i ) d-1
, where we adapt the style of above and denote γ τ for (γ τ i ) i∈I .

We define the two unit normal vectors n a ± (i) at each planar fracture Γ i , such that n a + (i) +n a -(i) = 0 (cf. figure 1). We define the set of indices χ = {a + (i), a -(i) | i ∈ I}, such that #χ = 2#I. For ease of notation, we use the convention Γ a + (i) = Γ a + (i) = Γ i . Then, for a = a ± (i) ∈ χ, we can define the trace operator on Γ a :

γ a : H 1 (Ω \ Γ) → L 2 (Γ a ),
and the normal trace operator on Γ a outward to the side a:

γ n,a : H div (Ω \ Γ) → D (Γ a ),
that satisfy γ a (h) = γ a (h ωa ) and γ n,a (p) = γ n,a (p ωa ), where

ω a = {x ∈ Ω | (x -y) • n a < 0, ∀y ∈ Γ i }.
We now define the hybrid dimensional function spaces that will be used as variational spaces for the Darcy flow models in the next subsection:

V = H 1 (Ω \ Γ) × H 1 (Γ),
and its subspace

V 0 = H 1 ∂Ω (Ω \ Γ) × H 1 Σ 0 (Γ)
, where (with γ ∂Ω : H 1 (Ω\Γ) → L 2 (∂Ω) denoting the trace operator on ∂Ω)

H 1 ∂Ω (Ω \ Γ) = {v ∈ H 1 (Ω\Γ) | γ ∂Ω v = 0 on ∂Ω}, as well as W = W m × W f ,
where

W m = q m ∈ H div (Ω \ Γ) | γ n,a q m ∈ L 2 (Γ a
) for all a ∈ χ and

W f = {q f = (q f,i ) i∈I | q f,i ∈ H div (Γ i ) for all i ∈ I and i∈Γ Γ i ∇ τ v • q f,i + v • div τ i q f,i dτ (x) = 0 for all v ∈ H 1 Σ 0 (Γ)}.
In the following, we will use the notation div τ p f = div τ i p f,i on Γ i for all i ∈ I and

p f = (p f,i ) i∈I ∈ W f .

Hybrid Dimensional Two Phase Darcy Flow Model

In the matrix domain Ω \ Γ, let us denote by Λ m ∈ L ∞ (Ω) d×d the permeability tensor such that there exist λ m ≥ λ m > 0 with

λ m |ζ| 2 ≤ (Λ m (x)ζ, ζ) ≤ λ m |ζ| 2 for all ζ ∈ R d , x ∈ Ω,
Analogously, in the fracture network Γ, we denote by Λ f ∈ L ∞ (Γ) (d-1)×(d-1) the tangential permeability tensor, and assume that there exist λ f ≥ λ f > 0, such that holds

λ f |ζ| 2 ≤ (Λ f (x)ζ, ζ) ≤ λ f |ζ| 2 for all ζ ∈ R d-1 , x ∈ Γ.
At the fracture network Γ, we introduce the orthonormal system (τ 1 (x), τ 2 (x), n(x)), defined a.e. on Γ. Inside the fractures, the normal direction is assumed to be a permeability principal direction. The normal permeability λ f,n ∈ L ∞ (Γ) is such that λ f,n ≤ λ f,n (x) ≤ λ f,n for a.e. x ∈ Γ with 0 < λ f,n ≤ λ f,n . We also denote by d f ∈ L ∞ (Γ) the width of the fractures, assumed to be such that there exist

d f ≥ d f > 0 with d f ≤ d f (x) ≤ d f for a.e.
x ∈ Γ. The half normal transmissibility in the fracture network is denoted by

T f = 2λ f,n d f .
Furthermore, φ m and φ f are the matrix and fracture porosities, respectively, ρ α denotes the density of phase α (with α = 1 the non-wetting and α = 2 the wetting phase) and g ∈ R d is the gravitational vector field. (k α m , k α f ) and (S α m , S α f ) are the matrix and fracture phase mobilities and saturations, respectively. We suppose that the matrix and the fracture network consist of a finite number of geological layers, that define finite partitions of Ω \ Γ and Γ. To identify the geological layers mathematically, we attribute a proper rocktype rt to each open set ω rt of these partitions. Then, we assume that on each ω rt , (k α m , k α f ) and (S α m , S α f ) are not explicitly space dependent. Moreover, on ω rt , (S 1 m (q m ), S 1 f (q f )) ∈ [0, 1] 2 for all (q m , q f ) ∈ R 2 and S 1 m , S 1 f are non-decreasing lipschitz continuous functions on R, and k α m , k α f are continuous, non-negative valued functions on [0, 1], for α = 1, 2. To simplify, we consider no sources or sinks.

The PDEs model writes:

find (u α m , u α f ) ∈ L 2 (0, T ; V 0 m ) × L 2 (0, T ; V 0 f ), (q α m , q α f ) ∈ L 2 (0, T ; W m ) × L 2 (0, T ; W f ), α = 1, 2, such that:        φ m ∂ t S α m (x, p m ) + div(q α m ) = 0 on Ω \ Γ q α m = -k α m (x, S α m (x, p m )) Λ m (∇u α m -ρ α g) on Ω \ Γ φ f d f ∂ t S α f (x, p f ) + div τ (q α f ) -a∈χ γ n,a q α m = 0 on Γ q α f = -d f k α f (x, S α f (x, p f )) Λ f (∇ τ u f -ρ α γ τ g) on Γ (1)
together with the coupling condition on Γ a , a ∈ χ

γ n,a q α m = k α f (x, S α f (x, γ a p m ))T f (γ a u α m -u α f - ρ α d f 2 γ n,a g) + + k α f (x, S α f (x, p f ))T f (γ a u α m -u α f - ρ α d f 2 γ n,a g) - (2)
where h + = max{0, h} and h -= -(-h) + (for any h), and initial values

(p m , p f )| t=0 = (p 0 m , p 0 f ) ∈ V 0 on (Ω \ Γ) × Γ.
Moreover, the saturations of both phases are coupled by the equation

(S 2 m , S 2 f ) = 1 -(S 1 m , S 1 f )
and the capillary pressure satisfies

(p m , p f ) = (u 1 m -u 2 m , u 1 f -u 2 f ). γ n a q α m q α m q α f,n a d f q α f S α f (p f ) S α f (γ a p m ) S α m (γ a p m )
Figure 2: Illustration of the coupling condition. It can be seen as an upwind two point approximation of q α f,na together with the flux conservation condition q α f,na = γ n,a q α m . The upwinding takes into account the saturation jumps due to discontinuous capillary pressure rules at the matrix-fracture interfaces.

Up to now, the only existing, comparable hybrid dimensional two-phase flow model is [START_REF] Jaffré | A discrete fracture model for two-phase flow with matrixfracture interaction[END_REF], which is presented in global pressure formulation and for only one fracture dividing the matrix domain. We adapted here a pressure-pressure formulation, accounting for complex fracture networks and general invertible capillary pressure functions. Another difference is, that the model presented here uses an upwind coupling condition for the matrix-fracture normal fluxes (see [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF]). This upwinding is necessary to transport the saturations from the matrix to the fractures. The coupling condition (2) also takes into account gravitational force inside the fractures for the matrix-fracture mass exchange. In the test cases below, we see that this is an important feature for the simulation of gravity dominant flow.

Vertex Approximate Gradient Scheme

In this section, the VAG scheme introduced in [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF] for diffusive problems on heterogeneous anisotropic media is extended to the hybrid dimensional model. We consider a finite volume version using lumping both for the source terms and the matrix fracture fluxes. Hence, the underlying discretization is non conforming w.r.t. the function space V 0 .

VAG Discretization

Generalized polyhedral meshes: Following [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF], we consider generalized polyhedral meshes of Ω. Let M be the set of cells that are disjoint open 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 . Let F denote the set of (not necessarily planar) faces of the mesh. We denote by V the set of vertices of the mesh. Let V K , F K , V σ respectively denote the set of the vertices of K ∈ M, faces of K, and vertices of σ ∈ F. For any face σ ∈ F K , we have V σ ⊂ V K . Let M s (resp. F s ) denote the set of the cells (resp. faces) sharing the vertex s ∈ V. The set of edges of the mesh is denoted by E and E σ denotes the set of edges of the face σ ∈ F. Let M σ denote the set of cells sharing the face σ ∈ F. We denote by F ext the subset of faces σ ∈ F such that M σ has only one element, and we set V ext = σ∈Fext V σ . The mesh is assumed to be conforming in the sense that for all σ ∈ F \ F ext , the set M σ contains exactly two cells. It is assumed that for each face σ ∈ F, there exists a so-called "center" of the face x σ such that

x σ = s∈Vσ β σ,s x s , with s∈Vσ β σ,s = 1,
where β σ,s ≥ 0 for all s ∈ V σ . The face σ is assumed to match with the union of the triangles T σ,e defined by the face center x σ and each of its edges e ∈ E σ . The mesh is assumed to be conforming w.r.t. the fracture network Γ in the sense that there exist subsets

F Γ i , i ∈ I of F such that Γ i = σ∈F Γ i σ.
We will denote by F Γ the set of fracture faces i∈I F Γ i . Similarly, we will denote by V Γ the set of fracture vertices σ∈F Γ V σ . We also define a submesh T of tetrahedra, where each tetrahedron D K,σ,e is the convex hull of the cell center x K of K, the face center x σ of σ ∈ F K and the edge e ∈ E σ . Similarly we define a triangulation ∆ of Γ, such that we have:

T = K∈F ,σ∈F K ,e∈Eσ {D K,σ,e } and ∆ = σ∈F Γ ,e∈Eσ {T σ,e }.
The mesh is also assumed to be conforming w.r.t. the partitions {ω rt } rt of Ω \ Γ and Γ, defined by the homogeneous geological layers covering the matrix and fracture domains, respectively (as described in the previous section). We therefore have a well defined rocktype rt for each cell and for each fracture face.

Degrees of freedom:

The set of matrix × fracture degrees of freedom is denoted by dof Dm ×dof D f . The real vector spaces X Dm and X D f of discrete unknowns in the matrix and in the fracture network respectively are then defined by

X Dm = span{e ν | ν ∈ dof Dm } X D f = span{e ν | ν ∈ dof D f }, where e ν = (δ νµ ) µ∈dof Dm for ν ∈ dof Dm (δ νµ ) µ∈dof D f for ν ∈ dof D f .
For u Dm ∈ X Dm and ν ∈ dof Dm we denote by u ν the νth component of u Dm and likewise for u D f ∈ X D f and ν ∈ dof D f . We also introduce the product of these vector spaces

X D = X Dm × X D f ,
for which we have dim X D = #dof Dm + #dof D f . To account for zero boundary conditions on ∂Ω and Σ 0 we introduce the subsets dof Dirm ⊂ dof Dm , and dof Dir f ⊂ dof D f , and we set dof Dir = dof Dirm × dof Dir f , and

X 0 D = {u ∈ X D | u ν = 0 for all ν ∈ dof Dir }.
Concretely, we consider the set of d.o.f. as illustrated in figure 3. Formally, for the matrix nodal unknowns, we first establish an equivalence relation on each M s , 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 Γ . Then we define (cf. figure 3)

dof Dm = M ∪ K σ σ ∈ F Γ , K ∈ M σ ∪ K s s ∈ V, K s ∈ M s , dof D f = F Γ ∪ V Γ , dof Dirm := K s s ∈ V ext , K s ∈ M s , dof Dir f = V Γ ∩ V ext .
Discrete gradients: The matrix discrete gradient ∇ Dm is defined on X Dm as the V 0 m conforming P 1 Finite Element gradient reconstruction on the tetrahedral submesh T , using barycentric interpolation to eliminate the d.o.f. at the non-fracture faces σ ∈ F\F Γ . The fracture discrete gradient ∇ D f is defined on X D f as the V 0 f conforming P 1 Finite Element gradient reconstruction on the triangulation ∆ of the DFN.

u M s = u N s u K σ u σ u L s u s u K u N u L u K s u M u L σ u K u s u K σ u K s u σ
∈ M let dof K = {K s , s ∈ V K } ∪ {K σ , σ ∈ F K ∩ F Γ }. Analogously, in the fracture domain, for σ ∈ F Γ let dof σ = V σ . We introduce, for any given K ∈ M, a partition {ω ν K } ν∈{K}∪dof K \dof Dir of K. Similarly, we define for any given σ ∈ F Γ a partition {ω ν σ } ν∈{σ}∪dofσ\dof Dir of σ. For each ν ∈ dof Dm , we define the open set ω ν = int K∈M ω ν K , with the convention ω ν K = ∅, if ν ∈ dof K . For each ν ∈ dof D f , we define the open set ω ν = int σ∈F ω ν σ , where ω ν σ = ∅, if ν ∈ dof σ .
We thus obtain the partitions {ω ν } ν∈dof Dm \dof Dirm of Ω and {ω ν } ν∈dof D f \dof Dir f of Γ. We also introduce for each T = T σ,s,s ∈ ∆ a partition {T i } i=1,...,3 of T , which we need for the definition of the VAG-CV matrix-fracture interaction operators. We assume that holds

|T 1 | = |T 2 | = |T 3 | = 1 3
|T | in order to preserve the first order accuracy of the scheme. Note that porosity is constant per cell and per fracture face, since we have a well defined value for the porosity for each rocktype and since the mesh is assumed to be conform with the partition in rocktypes (as described above). Therefore, in the numerical scheme, we do not need to reconstruct the just introduced partitions explicitly, but only have to define their corresponding volumes. Finally, we need a mapping between the degrees of freedom of the matrix domain, which are situated on one side of the fracture network, and the set of indices χ. For K σ ∈ dof Dm we have the one-element set

χ(K σ ) = {a ∈ χ | (x K -x σ ) • n a < 0} and therefore the notation a(K σ ) = a ∈ χ(K σ ).
The VAG-CV scheme's reconstruction operators are

• A function reconstruction operator on the matrix domain:

Π Dm u Dm = ν∈dof Dm \dof Dirm u ν 1 ων
• A function reconstruction operator on the fracture network:

Π D f u D f = ν∈dof D f \dof Dir f u ν 1 ων
• Reconstruction operators of the jump at the matrix fracture intersections on Γ a for a ∈ χ:

u D a,D = T σ,s,s ∈∆ K∈Mσ ((u Kσ -u σ )1 T 1 + (u Ks -u s )1 T 2 + (u K s -u s )1 T 3 )δ a(Kσ)a 1 Γa
• Reconstruction operators of the trace on Γ a for a ∈ χ:

T a Dm u Dm = T σ,s,s ∈∆ K∈Mσ (u Kσ 1 T 1 + u Ks 1 T 2 + u K s 1 T 3 )δ a(Kσ)a 1 Γa .

Finite Volume Formulation

Recall the definitions

dof K = {K s , s ∈ V K } ∪ {K σ , σ ∈ F K ∩ F Γ } for K ∈ M and dof σ = V σ for σ ∈ F Γ .
We introduce the family of rocktypes (rt ν ) ν∈M∪dof D f . Then, for any ν ∈ dof K the discrete matrix-matrix -fluxes are defined as

F α Kν (u 1 Dm , u 2 Dm ) = k α m (rt K , S α m (rt K , p K )) • f α Kν (u α Dm ) + + k α m (rt K , S α m (rt K , p ν )) • f α Kν (u α Dm ) -, where f α Kν (u α Dm ) = ν ∈dof K T νν K (u α K -u α ν -ρ α (x K -x ν ) • g),
with transmissivities

T νν K = K Λ m ∇ Dm e ν ∇ Dm e ν dx. It holds Ω Λ m (∇ Dm u α Dm -ρ α g)∇ Dm v Dm dx = K∈M ν∈dof K f α Kν (u α Dm )(v K -v ν )
, and where For all ν ∈ dof σ the discrete fracture-fracture-fluxes are defined as

F α σν (u 1 D f , u 2 D f ) = k α f (rt σ , S α f (rt σ , p σ )) • f α σ,ν (u α D f ) + + k α f (rt σ , S α f (rt σ , p ν )) • f α σ,ν (u α D f ) -, where f α σν (u α D f ) = ν ∈dofσ T νν σ (u α σ -u α ν -ρ α (x σ -x ν ) • g),
with transmissivities

T νν σ = σ Λ f ∇ D f e ν ∇ D f e ν dτ f (x). It holds Γ Λ f (∇ D f u α D f -ρ α γ τ g)∇ D f v D f dτ f (x) = σ∈F Γ ν∈dofσ f α σν (u α D f )(v σ -v ν ).
Let us further introduce the set of matrix-fracture (mf ) connectivities

C = {(ν m , ν f ) | ν m ∈ dof Γ Dm , ν f ∈ dof D f s.t. x νm = x ν f } with dof Γ Dm = {ν ∈ dof Dm | x ν ∈ Γ}.
The mf -fluxes are defined as

F α νmν f (u 1 D , u 2 D ) = k α f (rt ν f , S α f (rt ν f , p νm )) • f α νmν f (u α Dm , u α D f ) + + k α f (rt ν f , S α f (rt ν f , p ν f )) • f α νmν f (u α Dm , u α D f ) -, where f α νmν f (u α Dm , u α D f ) = T νmν f (u α νm -u α ν f - ρ α d f 2 γ n,a g),
with transmissivities

T νmν f = a∈χ Γa T f (T a Dm e νm ) 2 dτ (x). It holds (νm,ν f )∈C f α νmν f (u α Dm , u α D f )(v νm -v ν f ) = a∈χ Γa T f ( u α D a,D - ρ α d f 2 γ n,a g) v D a,D dτ (x), for all (v Dm , v D f ) ∈ X D .
We observe that for the VAG-CV scheme, the fluxes F νmν f are two point flux approximations.

Let 0 = t 0 < t 1 < • • • < t N = T , with ∆t n = t nt n-1 be a time discretization. Given p 0 D ∈ X D , the Finite Volume formulation of (1) reads as follows: Find

{u α D } α=1,2 ∈ (X 0 D ) 2N such that for all n ∈ {1, . . . , N } F K σ ,σ F K,K σ F σ,s u s u K s F K s ,s u σ u K σ u K s ′ u K F K,K s F K,K s ′
                                                                                                           for all K ∈ M : |ω K |φ K S α m (rt K ,p n K )-S α m (rt K ,p n-1 K ) ∆t n + ν∈dof K F α Kν (u 1,n Dm , u 2,n Dm ) = 0 for all ν m ∈ dof Dm \ (M ∪ dof Γ Dm ∪ dof Dirm ) : K∈Mν m |ω K,νm |φ K S α m (rt K ,p n νm )-S α m (rt K ,p n-1 νm ) ∆t n -F α Kνm (u 1,n Dm , u 2,n Dm ) = 0 for all ν m ∈ dof Γ Dm \ dof Dirm : ν f ∈dof D f s.t. (νm,ν f )∈C F α νmν f (u 1,n D , u 2,n D ) + K∈Mν m |ω K,νm |φ K S α m (rt K ,p n νm )-S α m (rt K ,p n-1 νm ) ∆t n -F α Kνm (u 1,n Dm , u 2,n Dm ) = 0
for all σ ∈ F Γ :

|ω σ |φ σ S α f (rtσ,p n σ )-S α f (rtσ,p n-1 σ ) ∆t n + ν∈dofσ F α σν (u 1,n D f , u 2,n D f ) - νm∈dof Dm s.t. (νm,σ)∈C F α νmσ (u 1,n D , u 2,n D ) = H α σ for all ν f ∈ dof D f \ (F Γ ∪ dof Dir f ) : σ∈F Γ,ν f |ω σ,ν f |φ σ S α f (rtσ,p n ν f )-S α f (rtσ,p n-1 ν f ) ∆t n -F α σν f (u 1,n D f , u 2,n D f ) - νm∈dof Dm s.t. (νm,ν f )∈C F α νmν f (u 1,n D , u 2,n D ) = 0. (3)
Here, M νm stands for the set of indices

{K ∈ M | ν m ∈ dof K } and F Γ,ν f stands for the set {σ ∈ F Γ | ν f ∈ dof σ }.

Two Phase Flow Test Cases

We present in this section a series of test cases for diphasic flow through a fractured 2 dimensional reservoir of geometry as shown in figure 5. The domain Ω is of extension (0, 400)m × (0, 800)m and the fracture width is assumed to be constantly d f = 4m. This choice for the width is motivated by the robustness of the equidimensional scheme (see below) and rather corresponds to the width of a fault (although we will keep the terminology fracture in the following). We consider isotropic permeability in the matrix and in the fractures. All tests have in common that initially, the reservoir is saturated with water (density 1000 kg m 3 , viscosity 0.001P a.s) and oil (density 700 kg m 3 , viscosity 0.005P a.s) is injected in the bottom fracture, which is managed by imposing non-homogeneous Neumann conditions at the injection location. The oil then rises by gravity, thanks to it's lower density compared to water and by the overpressure induced by the imposed injection rate. Also, Dirichlet boundary conditions are imposed at the upper boundary of the domain. Elsewhere, we have homogeneous Neumann conditions. The following test cases present a variety of geological and physical configurations in regard to matrix and fracture permeabilities and capillary pressure curves. We use the VAG discretization to obtain solutions for three different models for this two phase flow process. In the first model, fractures are represented as geological layers of equal dimension as the matrix and therefore, we refer to this model as the equi-dimensional model. The second model is the model we presented in the first part of this paper, refered to as discontinuous hybrid dimensional model, since pressure jumps at the matrix-fracture interfaces are allowed. The third model is the continuous hybrid dimensional model, presented in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF], which assumes pressure continuity accross the fractures.

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. Hence, we double the number of nodal unknowns, while keeping the number of cell and face unknowns constant (cf. table 1). This has to be kept in mind, when interpreting the computational cost. In order to account for the stratification of saturation in normal direction inside the fractures, which can play a major role in the flow process (see below), we need at least two layers of cells in the fractures for the equidim. model, to obtain valid reference solutions. Obviously, the larger number of cells for the equi-dimensional mesh is due to the need of tiny cells inside the DFN. In this regard, it is worth to mention that, with the hybrid dimensional model, the scale of fracture faces does not have a maximum constraint, other than with the equidimensional model, where the fracture width imposes an upper bound for the scale of faces between the matrix and the fracture, due to mesh regularity. However, all meshes are at fracture scale, here. The mesh for the two hybrid dimensional models is the same, but the number of degrees of freedom differs. The supplementary degrees of freedom for the discontinuous model are located at the matrix-fracture intersections and capture the pressure discontinuities, as described in the previous section.

The discrete problem is solved implicitly, where the non-linear system of equations occurring at each time step is solved via the Newton algorithm with relaxation. The stopping criterion is crit rel N ewton on the (L 1 ) relative residual. To ensure well defined values for the capillary pressure, after each Newton iteration, we project the (oil) saturation on the intervall [0, 1 -], with > 0 as small as desired. The resolution of the linear systems is performed by the GMRes solver (with stopping criterion crit rel GM Res on the relative residual), preconditioned by CPR-AMG. The time loop uses adaptive time stepping, i.e. the objective for the (max per d.o.f.) change in saturation per time step, ∆τ obj , is given and from this the time step is deduced under the condition that it does neither exceed a given maximal time step ∆t max nor 1.2 times the time step of the previous iteration. Also, if at a given time iteration the Newton algorithm does not converge after 35 iterations, then the actual time step is divided by 2 and the time iteration is repeated. The number of time step failures at the end of a simulation is indicated by N Chop .

Model

Comparisons between the equi and hybrid dimensional solutions for gravity dominated flow with zero capillary pressure

In this test case, we neglect capillary effects by setting the capillary pressure to zero. The following geological configuration is considered. In the matrix domain, permeability is isotropic of 0.1 Darcy and porosity is 0.2. In the DFN, permeability is isotropic of 100.0 Darcy and porosity is 0.4. Table 3: N ∆t is the number of successful time steps; N N ewton is the total number of Newton iterations (for successful time steps); N GM Res is the total number of GMRes iterations (for successful time steps); N Chop is the number of time step chops; CPU denotes the total cpu time in seconds.

This test case shows impressively, how the incorporation of normal fluxes at the mf intersections of the disc. hybrid dim. model allows to get much closer to the equidimensional reference solution than the cont. hybrid dimensional model does. The supplementary unknowns at the mf interfaces enables us to capture the segregation of saturation inside the DFN (due to gravity, here). In this view, the supplementary d.o.f. appear as a mesh refinement at the mf interfaces, that allows to reproduce the transport in normal direction to the DFN. In the gravity dominated test case shown in figure 6, this becomes particularly important, when gravitational acceleration is in a steep angle to the fracture network, which can be observed at the upper fracture. The drawback of this feature is that we have to deal with small volumes at the mf intesections, which is reflected in terms of computational cost, but the hybrid dim. model is still much cheaper than the full equidim. model. The absence of capillarity, of course, emphasisis this gap between the two hybrid dim. models, since at the mf interfaces, the matrix does not behave as a capillary barrier (saturation does not jump) and nothing holds back the oil from leaving the DFN. Also no capillary diffusion inside the fracture prevents the gravity segregation effect in the normal direction of the fracture. Therefore, in the next series of tests with capillary pressure, we might expect better match of the cont. hybrid dim. solution.

Comparisons between the equi and hybrid dimensional solutions for gravity dominated flow with discontinuous capillary pressure

The tests presented here account for capillarity. Inside the matrix domain the capillary pressure function is given by Corey's law p m = -a m log S 1 m . Inside the fracture network, we suppose p f = -a f log S 1 f . The hybrid dimensional model presented in the previous part of this paper is build to account for saturation jumps at the matrix-fracture interfaces (cf. figure 2). To treat the degenerated case of a f = 0, we adapt a novel variable switch technique presented 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]. This consists of introducing generalized variables as primary unknowns, that are used to parametrize the saturation and capillary pressure curves in order to avoid singularities at the matrix-fracture interfaces. As a counterpart, we had to regularize the accumulation terms at the mf interfaces, by adding small accumulation terms, driven by the fracture saturation

ν f ∈dof D f s.t. (νm,ν f )∈C K∈Mν m θ|ω K,νm |φ K S α f (rt ν f , p n νm ) -S α f (rt ν f , p n-1 νm ) ∆t n , θ ∈ R + , to the equations for ν m ∈ dof Γ Dm in (3)
, while keeping the conservation of volume.

drain-matrix permeability ratio of 1000

The geological setting is as follows. In the matrix domain, permeability is isotropic of 0.1 Darcy and porosity is 0.2. In the DFN, permeability is isotropic of 100.0 Darcy and porosity is 0.4. The Corey parameters are a m = 10 5 and a f = 0. Table 4: N ∆t is the number of successful time steps; N N ewton is the total number of Newton iterations (for successful time steps); N GM Res is the total number of GMRes iterations (for successful time steps); N Chop is the number of time step chops; CPU denotes the total cpu time in seconds.

drain-matrix permeability ratio of 100

The geological setting is as follows. In the matrix domain, permeability is isotropic of 0.1 Darcy and porosity is 0.2. In the DFN, permeability is isotropic of 10.0 Darcy and porosity is 0.4. We observe a degradation of the hybrid dimensional solutions w.r.t. the equidim. solution, when the mf permeability ratio decreases. This can be explained by the importance of fracture conductivity w.r.t. the diffusion in normal direction inside the fracture network, added artificially by the averaging procedure in the derivation of the reduced models. The more conductive the fractures, the more dominant the convective f f and mf fluxes and the less important the artificial diffusion. Figures 11 and12 reveil that the matching of equi-and hybrid dimensional solutions can be enhanced by adding capillarity in the DFN. More precisely, we note that the hybrid dim. solutions change insignificantly, but the equidim. solution changes towards the hybrid dim. solutions. Capillarity has a diffusive effect and smoothens out the stratification in the DFN, as shown in figure 13, which agrees better with the hybrid dimensional approach of averaging physical quantities over the fracture width. We also observe that the cont. hybrid dim. model simulates the global behaviour of the flow process quite well and is very competitive in view of cputime and robustness. However, compared to the disc. hybrid dim. model, the resolution at the DFN neighbourhood is much lower and local patterns can not be reproduced. It has been checked, that this is not an issue of mesh refinement. Rather, the adapted manner to approximate normal fluxes through the DFN of the disc. hybrid dim. model (as discussed in the previous test case) might play a role, here.

Comparisons between the equi and hybrid dimensional solutions for grav-

ity dominated flow with discontinuous capillary pressure at the matrix-drain interfaces and an upper barrier of matrix rocktype

In the matrix domain, permeability is isotropic of 0.1 Darcy and porosity is 0.2. The two lower fractures are drains of isotropic permeability 100.0 Darcy and porosity 0.4. In the upper fracture, permeability is isotropic of 0.001 Darcy and porosity is 0.2. Note that the continuous hybrid dimensional model does not incorporate a normal permeability in the DFN. We conducted the test case also for this model and observed, as expected, the unability to reproduce the barrier behaviour of the upper fracture. The Corey parameters are a m = a barrier = 10 5 and a drain = 0. Table 7: N ∆t is the number of successful time steps; N N ewton is the total number of Newton iterations (for successful time steps); N GM Res is the total number of GMRes iterations (for successful time steps); N Chop is the number of time step chops; CPU denotes the total cpu time in seconds.

Conclusion

The hybrid dimensional model for two phase flow through fractured porous media with pressure discontinuities at the mf intersecions presented here completes the literature by a pressure-pressure formulation of the problem with upwind mf fluxes -to account at the continuous model level for the transport from the matrix to the fracture -that integrate gravitational force inside the fractures -to account for gravity dominated flow. The only comparable predecessor [START_REF] Jaffré | A discrete fracture model for two-phase flow with matrixfracture interaction[END_REF] (in global pressure formulation with only one fracture dividing the matrix domain) has thus been extended to complex DFN and by the aforementioned features of mf fluxes. The Vertex Aproximate Gradient (VAG) scheme, as introduced in [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces[END_REF] for the monophasic stationary hybrid dimensional model, has been presented in a finite volume formulation for the two phase flow model. The VAG scheme is used to compare the numerically derived solutions of three different models for a 2D flow process through a fractured reservoir. More precisely, the discontinuous hybrid dimensional solution (model presented in this paper) has been compared to the continuous hybrid dimensional solution (cf. [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF]) w.r.t. a reference solution given by the equidimensional model (full model with fractures represented as heterogeneous layers), for a variety of geological and physical configurations in regard to matrix and fracture permeabilities and capillary pressure curves. Since the stratification in normal direction inside the fractures can play a major role, it is worth to mention that more than one layer of fracture cells is necessary in order to get valid reference solutions. The test cases show, that in terms of cpu time and robustness, the cont. hybrid dim. model has an advantage. Yet, the disc. hybrid dim. model still is much cheaper and more robust than the equidim. model. Moreover, for fracture apertures less than d f = 4m and fracture tangential permeabilities higher than λ f = 100.0 Darcy, the equidim. model is unpracticable. We observed that for high ratios of fracture and matrix permeabilities, the equi-and hybrid dimensional solutions match quite well and that for lower ratios, they differ more. This can be explained by the importance of fracture conductivity w.r.t. the diffusion in normal direction inside the fracture network, added artificially by the averaging procedure in the derivation of the reduced models. The more conductive the fractures, the more dominant the convective f f and mf fluxes and the less important the artificial diffusion. On the other hand, by adding capillarity in the DFN, the hybrid dimensional solutions fit much more to the equidim. solution. In fact, the equidim. solution moves towards the hybrid dim. solutions. In the first test case, gravitational segregation has a major influence on the global flow behaviour. This effect cannot be reproduced by the cont. hybrid dim. model, with single unknowns at the DFN. This remark applies to any cell centered scheme. The supplementary unknowns at the mf interfaces of the disc. hybrid dim. model enables us to capture gravitational segregation inside the DFN, which allows to be much more precise on the transport accross the DFN. This advantage becomes most striking, when acceleration acts in normal direction to fractures. Due to the assumption of pressure continuity at the mf interfaces (and the induced absence of λ f,n as a model parameter), the cont. hybrid dim. model is unusable, when it comes to the simulation of barriers. In the barrier test case presented here, we see that the disc. hybrid dim. model performs well, both, in terms of accuracy and computational cost. In any case, we observed a significant gain in precision for the disc. hybrid dim. solution w.r.t. the equidim. reference solution, compared to the cont. hybrid dim solution.

Figure 3 :

 3 Figure 3: Left: Illustration of d.o.f. in 2D for four cells intersected by three fractures (thick lines). Right: 3D cell K touching a fracture face σ. Illustration of the simplices on which ∇ Dm is constant (red) and ∇ D f is constant (grey). The facial d.o.f. in light grey is eliminated by barycentric interpolation.

Figure 4 :

 4 Figure 4: VAG mm-fluxes (red), mf -fluxes (dark red) and ff -fluxes (black) on a 3D cell touching a fracture

Figure 5 :

 5 Figure 5: Geometry of the reservoir under consideration. DFN in red and matrix domain in blue. Ω = (0, 400)m × (0, 800)m and d f = 4m.

Figure 6 :Figure 7 :

 67 Figure 6: Comparison of the equi dimensional (first line), discontinuous hybrid dimensional (mid line) and continuous hybrid dimensional (last line) numerical solutions for oil saturation at times t = 360, 1800, 3600, 5400 days (from left to right).

ModelN

  ∆t N N ewton N GM Res N Chop CPU equi dim

Figure 8 :

 8 Figure 8: Comparison of the equi dimensional (first line), discontinuous hybrid dimensional (mid line) and continuous hybrid dimensional (last line) numerical solutions for oil saturation at times t = 360, 1800, 3600, 5400 days (from left to right).

  IN THE FRACTURES equidim. model hybrid dim. disc. model hybrid dim. cont. model

Figure 9 :

 9 Figure 9: Comparison of the equi dimensional and hybrid dimensional matrix and fracture volumes occupied by oil as a function of time.

Figure 10 :Figure 11 :

 1011 Figure 10: Comparison of the equi dimensional (first line), discontinuous hybrid dimensional (mid line) and continuous hybrid dimensional (last line) numerical solutions for oil saturation at times t = 360, 1800, 4320, 5400 days (from left to right).

Figure 13 :

 13 Figure 13: Zoom on bottom DFN. Comparison of the equi dimensional oil saturation stratification in the fractures for Corey parameters a f = 0 (left) and a f = 1.E4 (right) at time t = 360 days.

Figure 14 :

 14 Figure 14: Comparison of the equi dimensional (first line) and discontinuous hybrid dimensional (second line) numerical solutions for oil saturation at times t = 360, 1800, 3600, 5400 days (from left to right).

Figure 15 :

 15 Figure 15: Comparison of the equi dimensional and discontinuous hybrid dimensional numerical liquid pressure at time t = 5400 days (from left to right).

Figure 16 :

 16 Figure 16: Comparison of the equi dimensional and hybrid dimensional matrix and fracture volumes occupied by oil as a function of time.

Table 1 :

 1 Nb Cells Nb dof Nb dof el. Nb Cells is the number of cells of the mesh; Nb dof is the number of discrete unknowns; Nb dof el. is the number of discrete unknowns after elimination of cell unknowns without fill-in.The numerical parameters are chosen as follows:

	equi dim.		22477	45315	22838
	disc. hybrid		16889	35355	18466
	cont. hybrid	16889	34291	17402
	Model	crit rel N ewton	crit rel GM Res	∆τ obj ∆t max
	equi dim.		1.E-5	1.E-6	0.5	10d
	disc. hybrid		1.E-6	1.E-6	0.5	10d
	cont. hybrid		1.E-6	1.E-6	0.5	60d

Table 2 :

 2 Numerical parameters.

Table 5 :

 5 N ∆t is the number of successful time steps; N N ewton is the total number of Newton iterations (for successful time steps); N GM Res is the total number of GMRes iterations (for successful time steps); N Chop is the number of time step chops; CPU denotes the total cpu time in seconds.non-zero capillary pressure in the DFN: The Corey parameters are a m = 10 5 and a f = 10 4 . Comparison of the equi dimensional and hybrid dimensional matrix and fracture volumes occupied by oil as a function of time.

	Model	N ∆t N N ewton N GM Res N Chop CPU
	equi dim.	933	6552	82477	30	5048
	disc. hybrid 1182	5619	47697	19	3244
	cont. hybrid 149	1082	12559	0	553

Table 6 :

 6 N ∆t is the number of successful time steps; N N ewton is the total number of Newton iterations (for successful time steps); N GM Res is the total number of GMRes iterations (for successful time steps); N Chop is the number of time step chops; CPU denotes the total cpu time in seconds.

We will denote by dτ (x) the d -1 dimensional Lebesgue measure on Γ. On the fracture network Γ, we define the function space L 2 (Γ) = {v = (v i ) i∈I , v i ∈ L 2 (Γ i ), i ∈ I}, and its subspace H 1 (Γ) consisting of functions v = (v i ) i∈I such that v i ∈ H 1 (Γ i ), i ∈ I with continuous traces at the fracture intersections Σ i,j , j ∈ I \ {i}. We also define it's subspace with vanishing traces on Σ 0 , which we