Gradient discretization of Hybrid Dimensional Darcy Flows in Fractured Porous Media with discontinuous pressures at the matrix fracture interfaces

We investigate the discretization of Darcy flow through fractured porous media on general meshes. We consider a hybrid dimensional model, invoking a complex network of planar fractures. The model accounts for matrix-fracture interactions and fractures acting either as drains or as barriers, i.e. we have to deal with pressure discontinuities at matrix-fracture interfaces. The numerical analysis is performed in the general framework of gradient discretizations which is extended to the model under consideration. Two families of schemes namely the Vertex Approximate Gradient scheme (VAG) and the Hybrid Finite Volume scheme (HFV) are detailed and shown to satisfy the gradient scheme framework, which yields, in particular, convergence. Numerical tests confirm the theoretical results. Gradient Discretization; Darcy Flow, Discrete Fracture Networks, Finite Volume


Introduction
This work deals with the discretization of Darcy flows in fractured porous media for which the fractures are modelized 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. We consider the case for which the pressure can be discontinuous at the matrix fracture interfaces in order to account for fractures acting either as drains or as barriers as described in [10], [12] and [3]. In this paper, we will study the family of models described in [12] and [3].
It is also assumed in the following that the pressure is continuous at the fracture intersections. This corresponds to a ratio between the permeability at the fracture intersection and the width of the fracture assumed to be high compared with the ratio between the tangential permeability of each fracture and its length. We refer to [14] for a more general reduced model taking into account discontinuous pressures at fracture intersections in dimension d = 2.
The discretization of such hybrid dimensional Darcy flow model has been the object of several works. In [10], [11], [3] a cell-centered Finite Volume scheme using a Two Point Flux Approximation (TPFA) is proposed assuming the orthogonality of the mesh and isotropic permeability fields. Cell-centered Finite Volume schemes have been extended to general meshes and anisotropic permeability fields using MultiPoint Flux Approximations (MPFA) in [13], [16], and [2]. In [12], a Mixed Finite Element (MFE) method is proposed and a MFE discretization adapted to non-matching fracture and matrix meshes is studied in [6]. More recently the Hybrid Finite Volume (HFV) scheme, introduced in [8], has been extended in [27] for the non matching discretization of two reduced fault models. Also a Mimetic Finite Difference (MFD) scheme is used in [1] in the matrix domain coupled with a TPFA scheme in the fracture network. Discretizations of the related reduced model [28] assuming a continuous pressure at the matrix fracture interfaces have been proposed in [28] using a MFE method, in [20] using a Control Volume Finite Element method (CVFE), in [19] using the HFV scheme, and in [19,5] using an extension of the Vertex Approximate Gradient (VAG) scheme introduced in [7].
In terms of convergence analysis, the case of continuous pressure models at the matrix fracture interfaces [28] is studied in [19] for a general fracture network but the current state of the art for the discontinuous pressure models at the matrix fracture interfaces is still limited to rather simple geometries. Let us recall that the family of models introduced in [12] and [3] depends on a quadrature parameter denoted by ξ ∈ [ 1 2 , 1] for the approximate integration in the width of the fractures. Existing convergence analysis for such models cover the case of one non immersed fracture separating the domain into two subdomains using a MFE discretization in [12] or a non matching MFE discretization in [6] for the range ξ ∈ ( 1 2 , 1]. In [3], the case of one fully immersed fracture in dimension d = 2 using a TPFA discretization is analysed for the full range of parameters ξ ∈ [ 1 2 , 1].
The main goal of this paper is to study the discretizations of such models and their convergence properties by extension of the gradient scheme framework. The gradient scheme framework has been introduced in [7], [22], [21] to analyse the convergence of numerical methods for linear and nonlinear second order diffusion problems. As shown in [22], this framework accounts for various conforming and non conforming discretizations such as Finite Element methods, Mixed and Mixed Hybrid Finite Element methods, and some Finite Volume schemes like symmetric MPFA, the VAG schemes [7], and the HFV schemes [8].
Our extension of the gradient scheme framework to the hybrid dimensional Darcy flow model will account for general fracture networks including fully, partially and non immersed fractures as well as fracture intersections in a 3D surrounding matrix domain. Each individual fracture will be assumed to be planar. The framework will cover the range of parameters ξ ∈ ( 1 2 , 1] excluding the value ξ = 1 2 in order to allow for a primal variational formulation. Two examples of gradient discretizations will be provided, namely the extension of the VAG and HFV schemes defined in [7] and [8] to the family of hybrid dimensional Darcy flow models. In both cases, it is assumed that the fracture network is conforming to the mesh in the sense that it is defined as a collection of faces of the mesh. The mesh is assumed to be polyhedral with possibly non planar faces for the VAG scheme and planar faces for the HFV scheme. Two versions of the VAG scheme will be studied, the first corresponding to the conforming P 1 finite element on a tetrahedral submesh, and the second to a finite volume scheme using lumping for the source terms as well as for the matrix fracture fluxes. The VAG scheme has the advantage to lead to a sparse discretization on tetrahedral or mainly tetrahedral meshes. It will be compared to the HFV discretization using face and fracture edge unknowns in addition to the cell unknowns. Note that the HFV scheme of [8] has been generalized in [23] as the family of Hybrid Mimetic Mixed methods which which encompasses the family of MFD schemes [24]. In this article, we will focus without restriction on the particular case presented in [8] for the sake of simplicity.
In section 2 we introduce the geometry of the matrix and fracture domains and present the strong and weak formulation of the model. Section 3 is devoted to the introduction of the general framework of gradient discretizations and the derivation of the error estimate 3.3. In section 4 we define and investigate the families of VAG and HFV discretizations. Having in mind applications to multi-phase flow, we also present a Finite Volume formulation involving conservative fluxes, which applies for both schemes. In section 5, the VAG and HFV schemes are compared in terms of accuracy and CPU efficiency for both Cartesian and tetrahedral meshes on hererogeneous isotropic and anisotropic media using a family of analytical solutions.

Hybrid dimensional Darcy Flow 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 For all i ∈ I, let us set Σ i = ∂Γ i , with n Σ i as unit vector in P i , normal to Σ i and outward Figure 1: Example of a 2D domain Ω and 3 intersecting fractures Γ i , i = 1, 2, 3. We might define the fracture plane orientations by α We will denote by dτ (x) the d − 1 dimensional Lebesgue measure on Γ. On the fracture network Γ, we define the function space 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 ).
We assume that there exists a finite family (Γ α ) α∈χ such that for all α ∈ χ holds: Γ α ⊂ Γ and there exists a lipschitz domain ω α ⊂ Ω \ Γ, such that Γ α = ∂ω α ∩ Γ. For α ∈ χ and an apropriate choice of I α ⊂ I we assume that Γ α = i∈Iα Γ i . Furthermore should hold Γ = α∈χ Γ α . We also assume that each Γ i ⊂ Γ is contained in Γ α for exactly two α ∈ χ and that we can define a unique mapping figure 1). For all i ∈ I, α ± (i) defines the two sides of the fracture Γ i in Ω \ Γ and we can introduce the corresponding unit normal vectors n α ± (i) at Γ i outward to ω α ± (i) , such that n α + (i) + n α − (i) = 0. We therefore obtain for α ∈ χ and a.e. x ∈ Γ α a unique unit normal vector n α (x) outward to ω α . A simple choice of (Γ α ) α∈χ is given by both sides of each fracture i ∈ I but more general choices are also possible such as for example the one exhibited in figure 1.
Then, for α ∈ χ, we can define the trace operator on Γ α : and the normal trace operator on Γ α outward to the side α: We now define the hybrid dimensional function spaces that will be used as variational spaces for the Darcy flow model in the next subsection: , where (with γ ∂Ω : H 1 (Ω\Γ) → L 2 (∂Ω) denoting the trace operator on ∂Ω) On V , we define the positive semidefinite, symmetric bilinear form Note that (·, ·) V is a scalar product and | · | V is a norm on V 0 , denoted by · V 0 in the following.
We define for all (p m , p f ), (q m , q f ) ∈ W the scalar product which induces the norm (q m , q f ) W , and where we have used the notation div Using similar arguments as in the proof of [15], example II.3.4, one can prove the following Poincaré type inequality.
Proof We apply the ideas of the proof of [15], example II.3.4 and assume that the statement of the proposition is not true. Then we can define a sequence (v l ) l∈N in V 0 , such that where, for this proof, · H 1 = · H 1 (Ω\Γ) + · H 1 (Γ) . The imbedding is compact, provided that Ω\Γ has the cone property (see [18], theorem 6.2). Thus, there is a On the other hand it follows from (2) that Remark 2.1 With the precedent proof it is readily seen that inequality (1) holds for all functions v ∈ V whose trace vanishes on a subset of ∂(Ω\Γ) with positive surface measure. The requirement is that v has to be in a closed subspace of (V, · H 1 ) for which · V 0 is a well defined norm.
The convergence analysis presented in section 4 requires some results on the density of smooth subspaces of V and W , which we state below.
is the set of functions ϕ, such that for all x ∈ Ω there exists r > 0, such that for all connected components ω of {x + y ∈ R d | |y| < r} ∩ (Ω \ Γ) one has ϕ ∈ C ∞ (ω).

C
Let us first state the following Lemma that will be used to prove the density of Proof Firstly, for all q m ∈ C ∞ 0 (Ω\Γ) d , we have and therefore v m ∈ H 1 (Ω\Γ) and ∇v m = G.
Proof Since W f is a closed subspace of the Hilbert space i∈I H div (Γ i ), any linear form l ∈ W f is the restriction to W f of a linear form still denoted by l in i∈I H div (Γ i ) . Then, for some f ∈ L 2 (Γ) and g ∈ L 2 (Γ) d−1 holds for all q f ∈ W f . Let us assume now that < l, ϕ >= 0 for all ϕ ∈ C ∞ W f . Corresponding to for all q m ∈ W m . Furthermore, let us assume that < l, ϕ >= 0 for all ϕ ∈ C ∞ Wm . From Lemma 2.1 we deduce that f ∈ H 1 ∂Ω (Ω \ Γ), that g = ∇f and that h α = γ α f (α ∈ χ). Using this, we conclude, again by the rule of partial integration, that < l, q m >= 0 for all q m ∈ W m .

Strong formulation
In the matrix domain Ω \ Γ, let us denote by Λ m ∈ L ∞ (Ω) d×d the permeability tensor such that there exist λ m ≥ λ m > 0 with 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 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 . We consider the source terms h m ∈ L 2 (Ω) (resp. h f ∈ L 2 (Γ)) in the matrix domain Ω \ Γ (resp. in the fracture network Γ). The half normal transmissibility in the fracture network is denoted by

Weak formulation
The hybrid dimensional weak formulation amounts to find (u m , u f ) ∈ V 0 satisfying the following variational equality for all The following proposition states the well posedness of the variational formulation (5).
Proof Using that for all ξ ∈ ( 1 2 , 1] and for all (a, b) ∈ R 2 one has the Lax-Milgram Theorem applies, which ensures the statement of the proposition.

Gradient Scheme Framework
A gradient discretization D of hybrid dimensional Darcy flow models is defined by a vector space of degrees of freedom X D = X Dm × X D f , its subspace satisfying ad hoc homogeneous boundary conditions X 0 D = X 0 Dm × X 0 D f , and the following gradient and reconstruction operators: • Gradient operator on the matrix domain: • Gradient operator on the fracture network: • A function reconstruction operator on the matrix domain: • Two function reconstruction operators on the fracture network: • Reconstruction operators of the trace on Γ α for α ∈ χ: which is assumed to define a norm on X 0 D .
The following properties of gradient discretizations are crucial for the convergence analysis of the corresponding numerical schemes: Coercivity: Let D be a gradient discretization and and be two sequences of gradient discretisations of (5) and let us assume that (D l ) l∈N is coercive, consistent and limit conforming. Let us furthermore assume that the sequence (ζ D l ,D l ) l∈N , defined by and that there is a constant C ∈ R independent of l such that for all v D l ∈ X 0 D l , l ∈ N. Then (D l ) l∈N is coercive, consistent and limit conforming.
In the last inequality we have used that v D D ≤ max(1, C) v D D , which follows from (7).
Consistency: Let l ∈ N be fixed and D = D l , D = D l . We first choose, for a given u = which holds for all v D ∈ X D , we obtain which implies that v D l D l is uniformly bouded and therefore S D l (u) → 0 as l → ∞.
Limit Conformity: Let again l ∈ N be fixed and Taking (7) into account, we derive Therefore W D l (q m , q f ) tends to zero as l goes to infinity.
Proposition 3.1 (Regularity at the Limit) Let (D l ) l∈N be a coercive and limit conforming sequence of gradient discretizations and let Proof By definition of the norm of X 0 D l and by coercivity, Using limit conformity we obtain (by letting l → ∞) The statement of the proposition follows now from Lemma 2.1.
Corollary 3.1 Let (D l ) l∈N be a sequence of gradient discretizations, assumed to be limit con- Then holds the conclusion of Proposition 3.1.

Application to (5)
The non conforming discrete variational formulation of the model problem is defined by: Proposition 3.2 Let ξ ∈ ( 1 2 , 1] and D be a gradient discretization, then (9) has a unique solution (u Dm , u D f ) ∈ X 0 D satisfying the a priori estimate Proof The Lax-Milgram Theorem applies, which ensures this result.
The main theoretical result for gradient schemes is stated by the following proposition: , and λ f,n such that one has the following error estimate: Proof From the definition of W D , and using the definitions (4) of the solution u, q and (9) of the discrete solution u D , it holds for all Let (10). Then holds with a constant C > 0 depending only on ξ, λ m , λ f ,λ m , λ f , d f , d f , λ f,n , and λ f,n . Taking coercivity into account leads to the statement of the proposition.

Two Examples of Gradient Schemes
Following [7], we consider generalised 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 faces of the mesh. The faces are not assumed to be planar for the VAG discretization, hence the term "generalised polyhedral cells", but they need to be planar for the HFV discretization. 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 F e denote the set of faces sharing the edge e ∈ E, and 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 E ext = σ∈Fext E σ , and 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 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 edge e ∈ E σ .
The mesh is assumed to be conforming w.r.t. the fracture network Γ in the sense that there exist subsets We will denote by F Γ the set of fracture faces i∈I F Γ i . Similarly, we will denote by E Γ the set of fracture edges σ∈F Γ E σ and 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: We introduce for D ∈ T the diameter h D of D and set h T = max D∈T h D . The regularity of our polyhedral mesh will be measured by the shape regularity of the tetrahedral submesh defined by 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 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 To account for our homogeneous 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

Vertex Approximate Gradient Discretization
In this subsection, the VAG discretization introduced in [7] for diffusive problems on heterogeneous anisotropic media is extended to the hybrid dimensional model. We consider the P 1 finite element construction as well as a finite volume version using lumping both for the source terms and the matrix fracture fluxes.
We first establish an equivalence relation on each M s , s ∈ V, by K ≡ Ms L ⇐⇒ there exists n ∈ N and a sequence ( We thus have Now we can introduce the piecewise affine interpolators (or reconstruction operators) which act linearly on X Dm and X D f , such that Π T u Dm is affine on each D K,σ,e ∈ T and satisfies on each cell K ∈ M where x ν ∈ Ω is the grid point associated with the degree of freedom ν ∈ dof Dm ∪ dof D f . The discrete gradients on X Dm and X D f are subsequently defined by We define the VAG-FE scheme's reconstruction operators by For the family of VAG-CV schemes, reconstruction operators are piecewise constant. We introduce, for any given K ∈ M, a partition Similarly, we define for any given σ ∈ F Γ a partition Similarly, for all s ∈ V Γ \ V ext we define ω s by We obtain the partitions We also introduce for each T = T σ,s,s ∈ ∆ a partition T = 3 i=1 T i , 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 convergence of the scheme.
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 σ ) = {α ∈ χ | n K,σ = n α on σ} and therefore the notation α(K σ ) = α ∈ χ(K σ ).
The VAG-CV scheme's reconstruction operators are Proposition 4.1 Let us consider a sequence of meshes (M l ) l∈N and let us assume that the sequence (T l ) l∈N of tetrahedral submeshes is shape regular, i.e. θ T l is uniformly bounded. We also assume that lim l→∞ h T l = 0. Then, the corresponding sequence of gradient discretizations (D l ) l∈N , defined by (12), (13), (14), is coercive, consistent and limit conforming.
Proof The VAG-FE scheme's reconstruction operators are conforming, i.e. V D ⊂ V 0 . Therefore we deduce coercivity from Proposition 2.1. Furthermore we have by partial integration W D (q m , q f ) = 0 for all (q m , q f ) ∈ W . Hence (D l ) l∈N is limit conforming.
To prove consistency, we need the following prerequisites. We define the linear mapping P Dm : C ∞ Ω → X 0 Dm such that for all ψ m ∈ C ∞ Ω and any cell K ∈ M one has Likewise, we define the linear mapping P D f : C ∞ Γ → X 0 D f such that for all ψ f ∈ C ∞ Γ holds (P D f ψ f ) ν = ψ f (x ν ) for all ν ∈ dof D f . It follows from the classical Finite Element approximation theory and from the fact that the interpolation s∈Vσ β σ,s (P Dm ψ m ) Ks at the point x σ , σ ∈ F K \F Γ is exact on cellwise affine functions, that for all (ψ m , The trace inequality implies that for all v ∈ H 1 ∂Ω (Ω\Γ) holds We can then calculate for (u m , u f ) ∈ C ∞ Ω × C ∞ Γ : Proposition 4.2 Let us consider a sequence of meshes (M l ) l∈N and let us assume that the sequence (T l ) l∈N of tetrahedral submeshes is shape regular, i.e. θ T l is uniformly bounded. We also assume that lim l→∞ h T l = 0. Then, any corresponding sequence of gradient discretizations (D l ) l∈N , defined by (12), (13), (15), is coercive, consistent and limit conforming. For the following, we define F α = i∈Iα F Γ i and V α = σ∈F α V σ . To ease the notation in the proof, we will use, for α ∈ χ, the uniquely identified mapping µ α : V α ∪ F α ⊂ dof D f → dof Dm , defined by µ α (σ) = K σ (such that χ(K σ ) = {α}) and µ α (s) = K s (for a cell K such that K ∈ M σ with σ ∈ F α ∩ F s and χ(K σ ) = {α}). Let now α ∈ χ be fixed. Since the mesh is conforming with respect to the fracture network, there is for every σ ∈ F α , e = ss ∈ E σ a ν(σ, e) ∈ {σ, s, s }, such that Then we have We have to check (6) now. It can be verified that [4], Lemma 3.4 applies to our case, both, in the matrix domain, where face unknowns might occur, as well as in the fracture network, a domain of codimension 1. This means that we can state that there exist constants C m (θ T ), C f (θ T ) > 0, such that (17) For the following calculation we take into account [4], Lemmata 3.2 and 3.4. We also use that the mesh is conforming with respect to the fracture network and that for σ ∈ F and K ∈ M σ (or equivalently for K ∈ M, σ ∈ F K ) holds: h K is asymptotically equivalent to h σ and |K| is asymptotically equivalent to h σ |σ|, where h K := max T D⊂K h D and h σ := max ∆ T ⊂σ h T . Let α ∈ χ, σ ∈ F α and K ∈ M σ , such that χ(K σ ) = {α}. Then we have Altogether we obtain with a constant C depending only on #χ and θ T . This proves that (6) is satisfied.
However, we can prove a higher order of convergence, i.e.

while (16) grants that holds
Limit Conformity: For all T ∈ ∆ and for all u Dm ∈ X Dm we have that Introducing the linear operator P : L 2 (Γ α ) → L 2 (Γ α ) such that P (ϕ) = 1 |T | T ϕdτ (x) on T for all T ∈ ∆, we first calculate for any q m ∈ C ∞ Wm γ n,α q m − P (γ n,α q m ) 2 We proceed: for all q m ∈ C ∞ Wm , where we have used (19) in the last inequality. We can now conclude by calculating for all for q = (q m , where we have taken into account the conformity of D in the first equation and (17), (18) in the last inequality.

Remark 4.2
The proofs of Propositions 4.1 and 4.2 show that for solutions (u m , u f ) ∈ V 0 and (q m , q f ) ∈ W of (4) such that u m ∈ C 2 (K), u f ∈ C 2 (σ), q m ∈ (C 1 (K)) d , q f ∈ (C 1 (σ)) d−1 for all K ∈ M and all σ ∈ Γ f , the VAG schemes are consistent and limit conforming of order 1, and therefore convergent of order 1.

Hybrid Finite Volume Discretization
In this subsection, the HFV scheme introduced in [8] is extended to the hybrid dimensional Darcy flow model. We assume here that the faces are planar and that x σ is the barycenter of σ for all σ ∈ F.
The set of indices dof Dm × dof D f for the unknowns is defined by (cf. figure 3) The discrete gradients in the matrix (respectively in the fracture domain) are defined in each cell (respectively in each face) by the 3D (respectively 2D) discrete gradients ∇ Dm (resp. ∇ D f ) as proposed in [8], pp. 8-9.
The function reconstruction operators are piecewise constant on a partition of the cells and of the fracture faces. These partitions are respectively denoted, for all K ∈ M, by and, for all σ ∈ F Γ , by With each σ ∈ F \ F ext and K σ ∈ M σ we associate an open set ω Kσ , s.t.
Similarly, for all e ∈ E Γ \ E ext we define ω e by We obtain the partitions Ω = We also 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 σ ∈ F Γ and K σ ∈ M σ holds by definition K σ = {K} for a K ∈ M σ and hence n Kσ = n K,σ is well defined. We obtain the one-element set χ(K σ ) = {α ∈ χ | n Kσ = n α on σ} and therefore the notation α(K σ ) = α ∈ χ(K σ ).
We define the HFV scheme's reconstruction operators by Proposition 4.3 Let us consider a sequence of meshes (M l ) l∈N and let us assume that the sequence (T l ) l∈N of tetrahedral submeshes is shape regular, i.e. θ T l is uniformly bounded. We also assume that lim l→∞ h T l = 0. Then, any corresponding sequence of gradient discretizations (D l ) l∈N , defined by (20), (21) and definition (22), is coercive, consistent and limit conforming.
Proof Let us denote in the following by Π M and Π F = Π F the HFV matrix and fracture reconstruction operators for the special case that ω Kσ = ∅ = ω e for all K σ ∈ σ∈F M σ and e ∈ E Γ . We start our numerical analysis for HFV by proving the proposition for these special choices and then use Lemma 3.1 for generalizing the results.
Coercivity: We first prove that limit conformity against regular test functions, as proved below, implies coercivity. Assume that the sequence of discretizations (D l ) l∈N is not coercive. Then we can find a sequence ( Then follows from a compactness result of [21] that there exists a u = (u m , u f ) ∈ L 2 (Ω)×L 2 (Γ), s.t. up to a subsequence and therefore u m L 2 (Ω) + u f L 2 (Γ) = 1. On the other hand follows from the discretizations' limit conformity against regular test functions (see below) by Proposition 3.1 and Corollary 3.1 that (u m , u f ) ∈ V 0 and that up to a subsequence Since by construction holds (u D l m , u D l f ) D l → 0, we obtain (u m , u f ) V 0 = 0. But · V 0 is a norm on V 0 , which contradicts the fact that u m L 2 (Ω) + u f L 2 (Γ) = 1.
Consistency: For (ϕ m , ϕ f ) ∈ C ∞ Ω × C ∞ Γ let us define the projection P Dm ϕ m ∈ X 0 Dm such that for all cell K ∈ M one has where C ϕm := max Ω ∇ϕ m . Summing over K ∈ M yields Analogously we can derive where c ϕ f := max Γ ∇ τ ϕ f . Furthermore, it follows from Lemma 4.3 of [8] that there exists C > 0 depending only on θ T and ϕ such that we see that the treated discretisation is consistent.
Limit Conformity: Let ϕ m ∈ C ∞ Wm and for all K ∈ M, σ ∈ F K let ϕ K := 1 |K| K ϕ m dx and ϕ K,σ := 1 |σ| σ γ n K,σ ϕ m dτ (x). In exactly the same manner as [19], (29)-(31) are proved, we can show that holds where with the definition of the gradient stabilization term R K,σ (u Dm ) as in [8], pp. 8-9. Therefore, applying Cauchy-Schwarz inequality to (25), using the regularity of ϕ m , and the estimate (24), we deduce that there exists C depending only on ϕ m , θ T , such that Taking into account the result [19] (33), i.e. for all ϕ ∈ C ∞ W f exists a constant C > 0 depending only on θ T , such that we obtain all together This result is shown above to imply coercivity, which is needed to conclude now. Finally, using that C ∞ Wm × C ∞ W f is dense in W and the coercivity of the scheme, we derive limit conformity on the whole space of test functions.
Generalization to arbitrary HFV discretizations: We want to apply Lemma 3.1. From [8] Lemma 4.1 and [21], it follows that there are positive constants C m and C f only depending on θ T and d, such that for all u D ∈ X D holds The remaining conditions of Lemma 3.1 are trivially satisfied, from what follows the statement of the proposition.

Remark 4.3
The precedent proof shows that for solutions (u m , u f ) ∈ V 0 and (q m , q f ) ∈ W of (4) such that u m ∈ C 2 (K), u f ∈ C 2 (σ), q m ∈ (C 1 (K)) d , q f ∈ (C 1 (σ)) d−1 for all K ∈ M and all σ ∈ Γ f , the HFV schemes are consistent and limit conforming of order 1, and therefore convergent of order 1.

Finite Volume Formulation for VAG and HFV Schemes
For K ∈ M let Analogously, in the fracture domain, for σ ∈ F Γ let Then, for any ν ∈ dof K the discrete matrix-matrix -fluxes are defined as . For all ν ∈ dof σ the discrete fracture-fracture-fluxes are defined as To take interactions of the matrix and the fracture domain into account we introduce the set of matrix-fracture (mf ) connectivities The mf -fluxes are built such that For all σ ∈ F Γ and K ∈ M σ , let us denote by α(K, σ) the unique α ∈ χ such that σ ∈ F α and n α = n K,σ . Let us also set for all σ ∈ F Γ , (χ × χ) σ = {(α(K, σ), α(L, σ)), (α(L, σ), α(K, σ))} with M σ = {K, L}. Then, holds For all σ ∈ F Γ , K ∈ M σ and x ∈ σ, let us notice that, for the VAG scheme, one has Π for all σ ∈ F Γ , M σ = {K, L} , and by for all s ∈ V Γ , Q s ∈ M s . Similarly the HFV matrix fracture fluxes are defined by We observe that for the VAG-CV scheme (since σ T f ( Π D f e s )( Π D f e s )dτ (x) = 0 for s = s and σ T f ( Π D f e σ )( Π D f e s )dτ (x) = 0) as well as for the HFV scheme, the fluxes F νmν f only involves the d.o.f. located at the point The discrete source terms are defined by The following Finite Volume formulation of (5) is equivalent to the discrete variational formulation (9) Here, M νm stands for the set of indices {K ∈ M | ν m ∈ dof K } and F Γ,ν f stands for the set It is important to note that, using the equation in each cell, the cell unknowns u K , K ∈ M, can be eliminated without fill-in.

Numerical Results
The objective of this numerical section is to compare the VAG-FE, VAG-CV, and the HFV schemes in terms of accuracy and CPU efficiency for both Cartesian and tetrahedral meshes on heterogeneous isotropic and anisotropic media. For that purpose a family of analytical solutions is built for the fixed value of the parameter ξ = 1. We refer to [12], [3], [2] for a comparison of the solutions obtained with different values of the parameter ξ ∈ [ 1 2 , 1] with the solution obtained with a 3D representation of the fractures. Table 1 exhibits for the Cartesian and tetrahedral meshes, as well as for both the VAG and HFV schemes, the number of degrees of freedom (Nb dof), the number of d.o.f. after elimination of the cell and Dirichlet unknowns (nb dof el.), and the number of nonzero element in the linear system after elimination without any fill-in of the cell unknowns (Nb Jac).
In all test cases, the linear system obtained after elimination of the cell unknowns is solved using the GMRes iterative solver with the stopping criteria 10 −10 . The GMRes solver is preconditioned by ILUT [25], [26] using the thresholding parameter 10 −4 chosen small enough in such a way that all the linear systems can be solved for both schemes and for all meshes. In tables 2 and 3, we report the number of GMRes iterations Iter and the CPU time taking into account the elimination of the cell unknowns, the ILUT factorization, the GMRes iterations, and the computation of the cell values.
We ran the program on a 2,6 GHz Intel Core i5 processor with 8 GB 1600 MHz DDR3 memory.
Derivation: For (u m , u f ) ∈ V , we denote u m (x, y, z) = u i (x, y, z) on Ω i , i = 1, . . . , 4 and u f (x, y, z) = u ij (y, z) on Γ ij , ij ∈ J, where we have introduced J = {12, 23, 34, 14}. We assume that a solution of the discontinuous pressure model writes in the fracture network u ij (y, z) = α f (z) + β ij (z)γ ij (y), ij ∈ J and in the matrix domain On γ ij , ij ∈ J we assume γ ij (0) = 0, such that the continuity of u f is well established at the fracture-fracture intersection, as well as γ ij (0) = 1, to ease the following calculations. For i = 1, . . . , 4 let K i = Λ m Ω i and for ij ∈ J let T ij = T f Γ ij . From the conditions γ n,α q m = T f (γ α u m − u f ) on Γ α , α ∈ χ, we then get, after some effort in computation,    Iter is the number of solver iterations; CPU refers to the solver CPU time in seconds; err sol , err grad , err jump are the respective L 2 -errors as defined above; α sol , α grad , α jump are the orders of convergence of the solution, of the gradient and of the jump, respectively.  Table 3: Anisotropic test case. Key refers to the mesh defined in table 1; Iter is the number of solver iterations; CPU refers to the solver CPU time in seconds; err sol , err grad , err jump are the respective L 2 -errors as defined above; α sol , α grad , α jump are the orders of convergence w.r.t. #M − 1 3 of the solution, of the gradient and of the jump, respectively.
The test case shows that, on cartesian grids, we obtain, as classically expected, convergence of order 2 for both, the solution and it's gradient. For tetrahedral grids, we obtain convergence of order 2 for the solution and convergence of order 1 for it's gradient. We observe that the VAG scheme is more efficient then the HFV scheme and this observation gets more obvious with increasing anisotropy. Comparing the precision of the discrete solution (and it's gradient) for VAG and HFV on a given mesh, we see that on hexahedral meshes, the advantage is on the side of VAG, whereas on tetrahedral meshes HFV is more precise (but much more expensive). On a given mesh, HFV is usually (see [19]) more accurate than VAG both for tetrahedral and hexahedral meshes. This is not the case for our test cases on Cartesian meshes maybe due to the higher number for VAG than for HFV of d.o.f. at the interfaces Γ α on the matrix side. It is also important to notice that there is literally no difference between VAG with finite element respectively lumped mf -fluxes concerning accuracy and convergence rate.

Conclusion
In this work, we extended the framework of gradient schemes (see [7]) to the model problem (4) of stationary Darcy flow through fractured porous media and gave numerical analysis results for this general framework.
The model problem (an extension to a network of fractures of a PDE model presented in [10], [12] and [3]) takes heterogeneities and anisotropy of the porous medium into account and involves a complex network of planar fractures, which might act either as barriers or as drains.
We also extended the VAG and HFV schemes to our model, where fractures acting as barriers force us to allow for pressure jumps across the fracture network. We developed two versions of VAG schemes, the conforming finite element version and the non-conforming control volume version, the latter particularly adapted for the treatment of material interfaces (cf. [9]). We showed, furthermore, that both versions of VAG schemes, as well as the proposed nonconforming HFV schemes, are incorporated by the gradient scheme's framework. Then, we applied the results for gradient schemes on VAG and HFV to obtain convergence, and, in particular, convergence of order 1 for "piecewise regular" solutions.
For implementation purposes and in view of the application to multi-phase flow, we also proposed a uniform Finite Volume formulation for VAG and HFV schemes. The numerical experiments on a family of analytical solutions show that the VAG scheme offers a better compromise between accuracy and CPU time than the HFV scheme especially for anisotropic problems.