# Gradient discretization of hybrid dimensional Darcy flows in fractured porous media

Konstantin Brenner, Mayya Groza, Cindy Guichard, Gilles Lebeau, Roland Masson

► To cite this version:

Konstantin Brenner, Mayya Groza, Cindy Guichard, Gilles Lebeau, Roland Masson. Gradient discretization of hybrid dimensional Darcy flows in fractured porous media. Numerische Mathematik, 2015, pp. 527-535. 10.1007/s00211-015-0782-x. hal-01097704v2

HAL Id: hal-01097704

https://hal.science/hal-01097704v2

Submitted on 1 Sep 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gradient discretization of Hybrid Dimensional Darcy Flows in Fractured Porous Media

K. Brenner,\* M. Groza,† C. Guichard,‡ G. Lebeau,§ R. Masson¶

September 1, 2015

## Abstract

This article deals with the discretization of hybrid dimensional Darcy flows in fractured porous media. These models couple the flow in the fractures represented as surfaces of codimension one with the flow in the surrounding matrix. The convergence analysis is carried out in the framework of gradient schemes which accounts for a large family of conforming and nonconforming discretizations. The Vertex Approximate Gradient (VAG) scheme and the Hybrid Finite Volume (HFV) scheme are extended to such models and are shown to verify the gradient scheme framework. Our theoretical results are confirmed by numerical experiments performed on tetrahedral, Cartesian and hexahedral meshes in heterogeneous isotropic and anisotropic porous media.

## 1 Introduction

This article 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 focus on the particular case where the pressure is continuous at the interfaces between the fractures and the matrix domain. This type of Darcy flow model introduced in [3] corresponds physically to 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. Note that it does not cover the case of fractures acting as barriers for which the pressure is discontinuous at the matrix fracture interfaces (see [17], [21], [4] for discontinuous pressure models). It is

<footnote>\*Laboratoire de Mathématiques J.A. Dieudonné, UMR 7351 CNRS, University Nice Sophia Antipolis, and team COFFEE, INRIA Sophia Antipolis Méditerranée, Parc Valrose 06108 Nice Cedex 02, France, konstantin.brenner@unice.fr</footnote>

<footnote>†Laboratoire de Mathématiques J.A. Dieudonné, UMR 7351 CNRS, University Nice Sophia Antipolis, and team COFFEE, INRIA Sophia Antipolis Méditerranée, Parc Valrose 06108 Nice Cedex 02, France, mayya.groza@unice.fr</footnote>

<footnote>‡Sorbonne Universités, UPMC Univ Paris 06, UMR 7598, CNRS, Laboratoire Jacques-Louis Li-ons, F-75005, Paris, and INRIA, ANGE project-team, Rocquencourt - B.P. 105, F78153 Le Chesnay Cedex, and CEREMA, ANGE project-team, 134 rue de Beauvais, F-60280 Margny-Lès-Compiègne, France, cindy.guichard@ann.jussieu.fr</footnote>

<footnote>§Laboratoire de Mathématiques J.A. Dieudonné, UMR 7351 CNRS, University Nice Sophia Antipolis, Parc Valrose 06108 Nice Cedex 02, France, gilles.lebeau@unice.fr</footnote>

<footnote>¶Laboratoire de Mathématiques J.A. Dieudonné, UMR 7351 CNRS, University Nice Sophia Antipolis, and team COFFEE, INRIA Sophia Antipolis Méditerranée, Parc Valrose 06108 Nice Cedex 02, France, roland.masson@unice.fr</footnote>
also assumed in our model that the pressure is continuous at the fracture intersections. It corresponds to a high ratio assumption between the permeability at the fracture intersections and the width of the fracture compared with the ratio between the tangential permeability of each fracture and its length. We refer to [20] for a more general reduced model taking into account discontinuous pressures at fracture intersections in dimension $d = 2$.

The discretization of the hybrid dimensional Darcy flow model with continuous pressures has been the object of several works. In [19] 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 can be extended to general meshes and anisotropic permeability fields using MultiPoint Flux Approximations (MPFA) following the ideas introduced in [30], [28], and [2] for discontinuous pressure models. In [3], a Mixed Finite Element (MFE) method is proposed, and Control Volume Finite Element Methods (CVFE) using nodal unknowns have been introduced for such models in [25] and [24]. A MFE discretization adapted to non-matching fracture and matrix grids is also studied in [10].

The main goal of this paper is to extend the gradient scheme framework to the case of hybrid dimensional Darcy flow models. This framework has been introduced in [14], [16], [11] to analyse the convergence of numerical methods for linear and nonlinear second order diffusion problems. As shown in [16], 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, Vertex Approximate Gradient (VAG) schemes [14], and Hybrid Finite Volume (HFV) schemes [13].

The extension of the gradient scheme framework to the hybrid dimensional Darcy flow model is defined by a vector space of degrees of freedom, two discrete gradient reconstruction operators and two discrete function reconstruction operators in the matrix and the fracture domains. The gradient discretization of the hybrid dimensional Darcy flow model is then based on a primal non conforming variational formulation using the previous operators. In the spirit of [14],[16] the well posedness and convergence of the gradient scheme is obtained assuming that the gradient discretization satisfies the so-called coercivity, consistency, and limit conformity assumptions.

Two examples of gradient discretization are given, namely we extend the VAG and HFV schemes defined in [14] and [13] to the hybrid dimensional Darcy flow model. 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 VAG scheme uses nodal and fracture face unknowns in addition to the cell unknowns which can be eliminated without any fill-in. It leads to a sparse discretization on tetrahedral or mainly tetrahedral meshes. It has the advantage, compared with CVFE approaches to avoid the mixing of the control volumes at the fracture matrix interfaces, which is a key feature for its application to multiphase Darcy flows (see [6]). It will be compared to the HFV discretization using face and fracture edge unknowns in addition to the cell unknowns which can be as well eliminated without any fill-in.

The proof that both the VAG and HFV schemes satisfy the coercivity, consistency, and limit conformity assumptions of the gradient scheme framework is based on a key result providing the density of smooth functions subspaces in both the variational space and in the flux space of the model. These density results are shown to hold for a general 3D network of possibly intersecting, immersed or non immersed planar fractures.

The outline of the paper is the following, in Section 2 we introduce the general 3D network of planar fractures, the function spaces, as well as the primal variational formulation of the hybrid dimensional Darcy flow model. Section 3 defines the gradient discretization framework stating the coercivity, consistency, limit conformity, and compactness assumptions. Then, the gradient scheme is introduced for the hybrid dimensional model and its well posedness and
convergence is shown to hold under the coercivity, consistency, and limit conformity assumptions. Section 4 extends the VAG and HFV schemes to our model and prove that each of them satisfies the gradient scheme assumptions. Section 5 proves the density of smooth functions subspaces in both the variational space and in the flux space which is the key ingredient to show that the VAG and HFV schemes satisfy the gradient scheme assumptions. Section 6 provides a numerical comparison of the VAG and HFV schemes on 3D analytical solutions using Cartesian, hexahedral and tetrahedral meshes. Both heterogeneous and anisotropic test cases are considered.

# 2 Hybrid dimensional Darcy Flow Model in Fractured Porous Media

## 2.1 Discrete Fracture Network and functional setting

Let $\Omega$ denote a bounded domain of $\mathbb{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.

We consider the asymptotic model introduced in [3] where fractures are represented as interfaces of codimension 1. Let $I$ be a finite set and let $\bar{\Gamma} = \bigcup_{i \in I} \bar{\Gamma}_i$ and its interior $\Gamma = \bar{\Gamma} \setminus \partial \bar{\Gamma}$ denote the network of fractures $\Gamma_i \subset \Omega$, $i \in I$, such that each $\Gamma_i$ is a planar polygonal simply connected open domain included in an oriented plane $\mathcal{P}_i$ of $\mathbb{R}^d$. It is assumed that the angles of $\Gamma_i$ are strictly smaller than $2\pi$ and that $\Gamma_i \cap \bar{\Gamma}_j = \emptyset$ for all $i \neq j$. For all $i \in I$, let us set $\Sigma_i = \partial \Gamma_i$, $\Sigma_{i,j} = \Sigma_i \cap \Sigma_j$, $j \in I \setminus \{i\}$, $\Sigma_{i,0} = \Sigma_i \cap \partial \Omega$, $\Sigma_{i,N} = \Sigma_i \setminus (\bigcup_{j \in I \setminus \{i\}} \Sigma_{i,j} \cup \Sigma_{i,0})$, and $\Sigma = \bigcup_{(i,j) \in I \times I, i \neq j} \Sigma_{i,j}$. It is assumed that $\Sigma_{i,0} = \bar{\Gamma}_i \cap \partial \Omega$. Let us define the following trace operator

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of a 2D domain with 3 intersecting fractures $\Gamma_1, \Gamma_2, \Gamma_3$ and 2 connected components $\Omega_1, \Omega_2$.</figcaption></figure>

$\gamma_{\Sigma_i} : H^1(\Gamma_i) \to L^2(\Sigma_i)$. We will denote by $d\tau(\mathbf{x})$ the $d-1$ dimensional Lebesgue measure on $\Gamma$. On the fracture network $\Gamma$, we define the function space $L^2(\Gamma) = \{v = (v_i)_{i \in I}, v_i \in L^2(\Gamma_i), i \in I\}$, endowed with the norm $\|v\|_{L^2(\Gamma)} = (\sum_{i \in I} \|v_i\|_{L^2(\Gamma_i)}^2)^{\frac{1}{2}}$. Its subspace $H^1(\Gamma)$ is defined as the space of functions $v = (v_i)_{i \in I}$ such that $v_i \in H^1(\Gamma_i)$, $i \in I$ with continuous traces at the fracture intersections i.e. $\gamma_{\Sigma_i} v_i = \gamma_{\Sigma_j} v_j$ on $\Sigma_{i,j}$ for all $i \neq j$ such that $\Sigma_{i,j}$ has a non zero $d-2$ Lebesgue measure. The space $H^1(\Gamma)$ is endowed with the norm $\|v\|_{H^1(\Gamma)} = (\sum_{i \in I} \|v_i\|_{H^1(\Gamma_i)}^2)^{\frac{1}{2}}$ and its subspace with vanishing traces on $\Sigma_0 = \bigcup_{i \in I} \Sigma_{i,0}$ is denoted by $H_{\Sigma_0}^1(\Gamma)$.

Let us also consider the trace operator $\gamma_i$ from $H^1(\Omega)$ to $L^2(\Gamma_i)$ as well as the trace operator $\gamma$ from $H^1(\Omega)$ to $L^2(\Gamma)$ such that $(\gamma v)_i = \gamma_i(v)$ for all $i \in I$.
On $\Omega$, the gradient operator from $H^1(\Omega)$ to $L^2(\Omega)^d$ is denoted by $\nabla$. On the fracture network $\Gamma$, the tangential gradient $\nabla_{\tau}$ acting from $H^1(\Gamma)$ to $L^2(\Gamma)^{d-1}$ is defined by

$$ \nabla_{\tau} v = (\nabla_{\tau_i} v_i)_{i \in I}, $$

where, for each $i \in I$, the tangential gradient $\nabla_{\tau_i}$ is defined from $H^1(\Gamma_i)$ to $L^2(\Gamma_i)^{d-1}$ by fixing a reference Cartesian coordinate system of the plane $\mathcal{P}_i$ containing $\Gamma_i$. We also denote by $\text{div}_{\tau_i}$ the divergence operator from $H_{\text{div}}(\Gamma_i)$ to $L^2(\Gamma_i)$.

The function spaces arising in the variational formulation of the hybrid dimensional Darcy flow model are

$$ V = \{v \in H^1(\Omega) \text{ such that } \gamma v \in H^1(\Gamma)\}, $$

and its subspace

$$ V^0 = \{v \in H_0^1(\Omega) \text{ such that } \gamma v \in H_{\Sigma_0}^1(\Gamma)\}. $$

From Poincaré inequality on $H_0^1(\Omega)$ and the continuity of the trace operator $\gamma$, we deduce the following inequality.

**Proposition 2.1** There exists $C_P > 0$ such that for all $v \in V^0$ one has

$$ \|v\|_{L^2(\Omega)} + \|\gamma v\|_{L^2(\Gamma)} \leq C_P \|\nabla v\|_{L^2(\Omega)^d}. $$

Thus, the space $V^0$ is endowed with the Hilbertian norm

$$ \|v\|_{V^0} = \left( \|\nabla v\|_{L^2(\Omega)^d}^2 + \|\nabla_{\tau} \gamma v\|_{L^2(\Gamma)^{d-1}}^2 \right)^{1/2}, $$

$$ \text{and the space } V \text{ with the Hilbertian norm } \|v\|_V = \left( \|v\|_{V^0}^2 + \|v\|_{L^2(\Omega)}^2 + \|\gamma v\|_{L^2(\Gamma)}^2 \right)^{1/2}. $$

The following density result is proved in Subsection 5.1.

**Proposition 2.2** The spaces $C^{\infty}(\bar{\Omega})$ and $C_c^{\infty}(\Omega)$ are dense subspaces of respectively $V$ and $V^0$.

Let $\Omega_\alpha, \alpha \in \mathcal{A}$ denote the connected components of $\Omega \setminus \bar{\Gamma}$, and let us define the space $H_{\text{div}}(\Omega \setminus \bar{\Gamma}) = \{\mathbf{q}_m = (\mathbf{q}_{m,\alpha})_{\alpha \in \mathcal{A}} | \mathbf{q}_{m,\alpha} \in H_{\text{div}}(\Omega_\alpha)\}$. For all $i \in I$, we can define the two sides $\pm$ of the fracture $\Gamma_i$ and the corresponding unit normal vector $\mathbf{n}_i^\pm$ at $\Gamma_i$ outward to the sides $\pm$. Each side $\pm$ corresponds to the subdomain $\alpha_i^\pm \in \mathcal{A}$ with possibly $\alpha_i^+ = \alpha_i^-$. For all $\mathbf{q}_m \in H_{\text{div}}(\Omega \setminus \bar{\Gamma})$, let $\mathbf{q}_{m,\alpha_i^\pm} \cdot \mathbf{n}_i^\pm|_{\Gamma_i}$ denote the two normal traces at the fracture $\Gamma_i$ and let us define the jump operator $H_{\text{div}}(\Omega \setminus \bar{\Gamma}) \to \mathcal{D}'(\Gamma_i)$ in the sense of distributions by $[\mathbf{q}_m \cdot \mathbf{n}_i] = \mathbf{q}_{m,\alpha_i^+} \cdot \mathbf{n}_i^+|_{\Gamma_i} + \mathbf{q}_{m,\alpha_i^-} \cdot \mathbf{n}_i^-|_{\Gamma_i}$.

For all $i \in I$, we denote by $\mathbf{n}_{\Sigma_i}$ the unit vector normal to $\Sigma_i$ outward (and tangent) to $\Gamma_i$.

Let us define the function space $H(\Omega, \Gamma)$ by

$$ H(\Omega, \Gamma) = \left\{ \begin{array}{ll} \mathbf{q}_m = (\mathbf{q}_{m,\alpha})_{\alpha \in \mathcal{A}}, & \mathbf{q}_f = (\mathbf{q}_{f,i})_{i \in I} | \mathbf{q}_m \in H_{\text{div}}(\Omega \setminus \bar{\Gamma}), \\ \mathbf{q}_{f,i} \in L^2(\Gamma_i)^{d-1}, & \text{div}_{\tau_i}(\mathbf{q}_{f,i}) - [\mathbf{q}_m \cdot \mathbf{n}_i] \in L^2(\Gamma_i), i \in I \end{array} \right\}. $$

It is an Hilbert space endowed with the scalar product

$$ \langle (\mathbf{p}_m, \mathbf{p}_f), (\mathbf{q}_m, \mathbf{q}_f) \rangle_H = \sum_{\alpha \in \mathcal{A}} \int_{\Omega_\alpha} (\mathbf{p}_{m,\alpha} \cdot \mathbf{q}_{m,\alpha} + \text{div}(\mathbf{p}_{m,\alpha})\text{div}(\mathbf{q}_{m,\alpha})) d\mathbf{x} \\ + \sum_{i \in I} \int_{\Gamma_i} (\mathbf{p}_{f,i} \cdot \mathbf{q}_{f,i} + (\text{div}_{\tau_i}(\mathbf{p}_{f,i}) - [\mathbf{p}_m \cdot \mathbf{n}_i]) (\text{div}_{\tau_i}(\mathbf{q}_{f,i}) - [\mathbf{q}_m \cdot \mathbf{n}_i])) d\tau(\mathbf{x}). $$
and the norm

$$ \|(\mathbf{q}_m, \mathbf{q}_f)\|_H = \langle (\mathbf{q}_m, \mathbf{q}_f), (\mathbf{q}_m, \mathbf{q}_f) \rangle_H^{1/2}. $$

On $H(\Omega, \Gamma) \times (L^2(\Omega) \times L^2(\Gamma) \times L^2(\Omega)^d \times L^2(\Gamma)^{d-1})$ we define the continuous bilinear form

$$ a_{\Sigma}((\mathbf{q}_m, \mathbf{q}_f), (v_m, v_f, \mathbf{g}_m, \mathbf{g}_f)) = \sum_{\alpha \in \mathcal{A}} \int_{\Omega_{\alpha}} (\mathbf{q}_{m,\alpha} \cdot \mathbf{g}_m + \operatorname{div}(\mathbf{q}_{m,\alpha})v_m) d\mathbf{x} \\ + \sum_{i \in I} \int_{\Gamma_i} (\mathbf{q}_{f,i} \cdot \mathbf{g}_f + (\operatorname{div}_{\tau_i}(\mathbf{q}_{f,i}) - [\mathbf{q}_m \cdot \mathbf{n}_i])v_f) d\tau(\mathbf{x}). \quad (1) $$

For all $\alpha \in \mathcal{A}$ let us denote by $C_b^\infty(\Omega_\alpha)$ the set of functions $\varphi$ such that for all $\mathbf{x} \in \bar{\Omega}_\alpha$, there exists $r > 0$ such that for all connected component $\omega$ of the domain $\{\mathbf{x} \in \mathbb{R}^d \mid |\mathbf{x}| < r\} \cap \Omega_\alpha$, one has $\varphi|_{\omega} \in C^\infty(\bar{\omega})^d$.

For all $(\mathbf{q}_m, \mathbf{q}_f) \in H(\Omega, \Gamma)$ with $\mathbf{q}_{m,\alpha} \in C_b^\infty(\Omega_\alpha)$, $\alpha \in \mathcal{A}$, and $\mathbf{q}_{f,i} \in C^\infty(\bar{\Gamma}_i)^{d-1}$, $i \in I$, and for all $v \in V^0$, it is readily checked that

$$ a_{\Sigma}((\mathbf{q}_m, \mathbf{q}_f), (v, \gamma v, \nabla v, \nabla_{\tau} \gamma v)) = \int_{\Sigma \setminus \Sigma_0} \gamma v (\sum_{i \in I} \mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i}) dl(\mathbf{x}) \\ + \sum_{i \in I} \int_{\Sigma_{i,N}} \gamma v (\mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i}) dl(\mathbf{x}). $$

This lead us to the definition of the following closed Hilbert subspace of $H(\Omega, \Gamma)$

$$ H_{\Sigma}(\Omega, \Gamma) = \left\{ \begin{array}{ll} (\mathbf{q}_m, \mathbf{q}_f) \in H(\Omega, \Gamma) & | \\ a_{\Sigma}((\mathbf{q}_m, \mathbf{q}_f), (v, \gamma v, \nabla v, \nabla_{\tau} \gamma v)) = 0 & \text{for all } v \in V^0 \end{array} \right\} \quad (2) $$

corresponding to impose in a weak sense the conditions $\sum_{i \in I} \mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i} = 0$ on $\Sigma \setminus \Sigma_0$ and $\mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i} = 0$ on $\Sigma_{i,N}, i \in I$.

Finally let us define a subspace of smooth functions in $H_{\Sigma}(\Omega, \Gamma)$ by

$$ C_{H_{\Sigma}}^{\infty}(\Omega, \Gamma) = \left\{ \sum_{i \in I} \mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i} = 0 \text{ on } \Sigma \setminus \Sigma_0, \mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i} = 0 \text{ on } \Sigma_{i,N}, i \in I \right\}. \quad (3) $$

The proof of the following density result is given in Subsection 5.2.

**Proposition 2.3** The space $C_{H_{\Sigma}}^{\infty}(\Omega, \Gamma)$ is a dense subspace of $H_{\Sigma}(\Omega, \Gamma)$.

## 2.2 Hybrid dimensional Darcy Flow Model

In the matrix domain $\Omega \setminus \bar{\Gamma}$ (resp. in the fracture network $\Gamma$), let us denote by $\Lambda_m \in L^\infty(\Omega)^{d \times d}$ (resp. $\Lambda_f \in L^\infty(\Gamma)^{(d-1) \times (d-1)}$) the permeability tensor such that there exist $\bar{\lambda}_m \ge \bar{\lambda}_f > 0$ (resp. $\bar{\lambda}_f \ge \bar{\lambda}_m > 0$) with

$$ \bar{\lambda}_m |\boldsymbol{\xi}|^2 \le (\Lambda_m(\mathbf{x})\boldsymbol{\xi}, \boldsymbol{\xi}) \le \bar{\lambda}_m |\boldsymbol{\xi}|^2 \text{ for all } \boldsymbol{\xi} \in \mathbb{R}^d, \mathbf{x} \in \Omega, $$

$$ (\text{resp. } \bar{\lambda}_f |\boldsymbol{\xi}|^2 \le (\Lambda_f(\mathbf{x})\boldsymbol{\xi}, \boldsymbol{\xi}) \le \bar{\lambda}_f |\boldsymbol{\xi}|^2 \text{ for all } \boldsymbol{\xi} \in \mathbb{R}^{d-1}, \mathbf{x} \in \Gamma). $$

We also denote by $d_f \in L^\infty(\Gamma)$ the width of the fractures assumed to be such that there exist $\bar{d}_f \ge \underline{d}_f > 0$ with $\underline{d}_f \le d_f(\mathbf{x}) \le \bar{d}_f$ for all $\mathbf{x} \in \Gamma$. Let us define the weighted Lebesgue $d-1$ dimensional measure on $\Gamma$ by $d\tau_f(\mathbf{x}) = d_f(\mathbf{x})d\tau(\mathbf{x})$. We consider the source terms $h_m \in L^2(\Omega)$ (resp. $h_f \in L^2(\Gamma)$) in the matrix domain $\Omega \setminus \bar{\Gamma}$ (resp. in the fracture network $\Gamma$).
The strong formulation of the hybrid dimensional Darcy flow model amounts to find $u \in V^0$ and $(\mathbf{q}_m, \mathbf{q}_f) \in H_\Sigma(\Omega, \Gamma)$ such that:

$$ egin{cases} 
\operatorname{div}(\mathbf{q}_{m,\alpha}) = h_m & \text{on } \Omega_\alpha, \alpha \in \mathcal{A}, \\
\mathbf{q}_{m,\alpha} = -\Lambda_m \nabla u & \text{on } \Omega_\alpha, \alpha \in \mathcal{A}, \\
\operatorname{div}_{\tau_i}(\mathbf{q}_{f,i}) - [\mathbf{q}_m \cdot \mathbf{n}_i] = d_f h_f & \text{on } \Gamma_i, i \in I, \\
\mathbf{q}_{f,i} = -d_f \Lambda_f \nabla_{\tau_i} \gamma_i u & \text{on } \Gamma_i, i \in I.
\end{cases} 
\quad (4) $$

The weak formulation of (4) amounts to find $u \in V^0$ satisfying the following variational equality for all $v \in V^0$:

$$ egin{cases} 
\int_{\Omega} \Lambda_m(\mathbf{x}) \nabla u(\mathbf{x}) \cdot \nabla v(\mathbf{x}) d\mathbf{x} + \int_{\Gamma} \Lambda_f(\mathbf{x}) \nabla_{\tau} \gamma u(\mathbf{x}) \cdot \nabla_{\tau} \gamma v(\mathbf{x}) d\tau_f(\mathbf{x}) \\
- \int_{\Omega} h_m(\mathbf{x}) v(\mathbf{x}) d\mathbf{x} - \int_{\Gamma} h_f(\mathbf{x}) \gamma v(\mathbf{x}) d\tau_f(\mathbf{x}) = 0.
\end{cases} 
\quad (5) $$

The following proposition is a direct application of the Lax-Milgram theorem and Proposition 2.1.

**Proposition 2.4** The variational problem (5) has a unique solution $u \in V^0$ which satisfies the a priori estimate

$$ \|u\|_{V^0} \le \frac{C_P}{\min(\lambda_m, \lambda_f \underline{d}_f)} (\|h_m\|_{L^2(\Omega)} + \|d_f h_f\|_{L^2(\Gamma)}). $$

In addition $(\mathbf{q}_m, \mathbf{q}_f) = (-\Lambda_m \nabla u, -d_f \Lambda_f \nabla_{\tau} \gamma u)$ belongs to $H_\Sigma(\Omega, \Gamma)$.

# 3 Gradient discretization of the hybrid dimensional model

In this section we extend the gradient scheme framework introduced in [14], [16] for elliptic and parabolic problems to our hybrid dimensional Darcy flow model.

## 3.1 Gradient discretization

A gradient discretization $\mathcal{D}$ of (5) is defined by a vector space of degrees of freedom $X_\mathcal{D}$, its subspace associated with homogeneous Dirichlet boundary conditions $X_\mathcal{D}^0$, and the following set of linear operators:

* Two discrete gradient operators on the matrix and fracture domains:

$$ \nabla_{\mathcal{D}_m} : X_{\mathcal{D}} \to L^2(\Omega)^d \text{ and } \nabla_{\mathcal{D}_f} : X_{\mathcal{D}} \to L^2(\Gamma)^{d-1}. $$

* Two function reconstruction operators on the matrix and fracture domains:

$$ \Pi_{\mathcal{D}_m} : X_{\mathcal{D}} \to L^2(\Omega) \text{ and } \Pi_{\mathcal{D}_f} : X_{\mathcal{D}} \to L^2(\Gamma). $$

The vector space $X_\mathcal{D}$ is endowed with the semi-norm

$$ \|v_\mathcal{D}\|_\mathcal{D} = \left( \|\nabla_{\mathcal{D}_m} v_\mathcal{D}\|_{L^2(\Omega)^d}^2 + \|\nabla_{\mathcal{D}_f} v_\mathcal{D}\|_{L^2(\Gamma)^{d-1}}^2 \right)^{\frac{1}{2}}, $$

which is assumed to define a norm on $X_\mathcal{D}^0$.

In the following, the gradient discretization of the hybrid dimensional model with homogeneous Dirichlet boundary conditions will be denoted by the quintuplet

$$ \mathcal{D} = (X_{\mathcal{D}}^0, \Pi_{\mathcal{D}_m}, \Pi_{\mathcal{D}_f}, \nabla_{\mathcal{D}_m}, \nabla_{\mathcal{D}_f}). $$
Next, we define the coercivity, consistency, limit conformity and compactness properties for sequences of gradient discretizations. Note that the compactness property is useful for the convergence analysis of nonlinear models and hence will not be used for the convergence analysis of our model. Nevertheless, for the sake of completeness, it is stated in this section and will be proved to be verified for the VAG and HFV schemes in section 4.

**Coercivity:** Let $C_D > 0$ be defined by

$$C_D = \max_{0 \neq v_D \in X_D^0} \frac{\|\Pi_{D_m} v_D\|_{L^2(\Omega)} + \|\Pi_{D_f} v_D\|_{L^2(\Gamma)}}{\|v_D\|_D} \quad (6)$$

Then, a sequence of gradient discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ is said to be coercive if there exist $\bar{C}_P > 0$ such that $C_{\mathcal{D}^l} \le \bar{C}_P$ for all $l \in \mathbb{N}$.

**Consistency:** For all $u \in V^0$ and $v_D \in X_D^0$ let us define

$$S_D(u, v_D) = \|\nabla_{D_m} v_D - \nabla u\|_{L^2(\Omega)^d} + \|\nabla_{D_f} v_D - \nabla_{\tau} \gamma u\|_{L^2(\Gamma)^{d-1}} \\ + \|\Pi_{D_m} v_D - u\|_{L^2(\Omega)} + \|\Pi_{D_f} v_D - \gamma u\|_{L^2(\Gamma)}, \quad (7)$$

and

$$S_D(u) = \min_{v_D \in X_D^0} S_D(u, v_D). \quad (8)$$

Then, a sequence of gradient discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ is said to be consistent if for all $u \in V^0$ one has $\lim_{l \to +\infty} S_{\mathcal{D}^l}(u) = 0$.

**Limit Conformity:** For all $(\mathbf{q}_m, \mathbf{q}_f) \in H_\Sigma(\Omega, \Gamma)$ and $v_D \in X_D^0$, let us define

$$W_D(\mathbf{q}_m, \mathbf{q}_f, v_D) = a_\Sigma \left( (\mathbf{q}_m, \mathbf{q}_f), (\Pi_{D_m} v_D, \Pi_{D_f} v_D, \nabla_{D_m} v_D, \nabla_{D_f} v_D) \right) \quad (9)$$

and

$$W_D(\mathbf{q}_m, \mathbf{q}_f) = \max_{0 \neq v_D \in X_D^0} \frac{|W_D(\mathbf{q}_m, \mathbf{q}_f, v_D)|}{\|v_D\|_D} \quad (10)$$

Then, a sequence of gradient discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ is said to be limit conforming if for all $(\mathbf{q}_m, \mathbf{q}_f) \in H_\Sigma(\Omega, \Gamma)$ one has $\lim_{l \to +\infty} W_{\mathcal{D}^l}(\mathbf{q}_m, \mathbf{q}_f) = 0$.

**Compactness:** A sequence of gradient discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ is said to be compact if for all sequences $(v_{\mathcal{D}^l})_{l \in \mathbb{N}}$ with $v_{\mathcal{D}^l} \in X_{\mathcal{D}^l}^0$ for all $l \in \mathbb{N}$ such that there exists $C > 0$ with $\|v_{\mathcal{D}^l}\|_{\mathcal{D}^l} \le C$ for all $l \in \mathbb{N}$, then there exist $u_m \in L^2(\Omega)$ and $u_f \in L^2(\Gamma)$ such that one has up to a subsequence

$$\lim_{l \to +\infty} \|\Pi_{D_m^l} v_{\mathcal{D}^l} - u_m\|_{L^2(\Omega)} = 0 \quad \text{and} \quad \lim_{l \to +\infty} \|\Pi_{D_f^l} v_{\mathcal{D}^l} - u_f\|_{L^2(\Gamma)} = 0.$$

The following proposition states a property of limit conforming and coercive sequences of gradient discretizations.

**Proposition 3.1 Regularity at the limit.** Let $(\mathcal{D}^l)_{l \in \mathbb{N}}$ be a family of discretizations assumed to be limit conforming and coercive. Let $(v_{\mathcal{D}^l})_{l \in \mathbb{N}}$ with $v_{\mathcal{D}^l} \in X_{\mathcal{D}^l}^0$ for all $l \in \mathbb{N}$ be a bounded sequence in the sense that there exists $C$ such that $\|v_{\mathcal{D}^l}\|_{\mathcal{D}^l} \le C$ for all $l \in \mathbb{N}$. Then, there exist $v \in V^0$ and a subsequence still denoted by $(v_{\mathcal{D}^l})_{l \in \mathbb{N}}$ such that

$$\left\{ \begin{array}{l} \Pi_{D_m^l} v_{\mathcal{D}^l} \to v \text{ weakly in } L^2(\Omega), \\ \nabla_{\mathcal{D}^l} v_{\mathcal{D}^l} \to \nabla v \text{ weakly in } L^2(\Omega)^d, \\ \Pi_{D_f^l} v_{\mathcal{D}^l} \to \gamma v \text{ weakly in } L^2(\Gamma), \\ \nabla_{D_f} v_{\mathcal{D}^l} \to \nabla_{\tau} \gamma v \text{ weakly in } L^2(\Gamma)^{d-1}. \end{array} \right.$$
**Proof:** From the boundedness of the sequence $||v_{\mathcal{D}^l}||_{\mathcal{D}^l}$, $l \in \mathbb{N}$ and the coercivity assumption, there exist $v_m \in L^2(\Omega)$, $v_f \in L^2(\Gamma)$, $g_m \in L^2(\Omega)^d$, and $g_f \in L^2(\Gamma)^{d-1}$ such that $\Pi_{\mathcal{D}_m^l} v_{\mathcal{D}^l}$ weakly converges to $v_m$ in $L^2(\Omega)$, $\Pi_{\mathcal{D}_f^l} v_{\mathcal{D}^l}$ weakly converges to $v_f$ in $L^2(\Gamma)$, $\nabla_{\mathcal{D}_m^l} v_{\mathcal{D}^l}$ weakly converges to $g_m$ in $L^2(\Omega)^d$, and $\nabla_{\mathcal{D}_f^l} v_{\mathcal{D}^l}$ weakly converges to $g_f$ in $L^2(\Gamma)^{d-1}$. From the limit conformity it follows that $a_\Sigma((q_m, q_f), (v_m, v_f, g_m, g_f)) = 0$ for all $(q_m, q_f) \in H_\Sigma(\Omega, \Gamma)$. From Lemma 5.5, we deduce that $v = v_m \in V^0$, $v_f = \gamma v$, $g_m = \nabla v$ and $g_f = \nabla_\tau \gamma v$. $\square$

The following Lemma will be used in the next sections to prove the coercivity, consistency, limit conformity and compactness of sequences of families of gradient discretizations.

**Lemma 3.1** Let $(\mathcal{D}^l)_{l \in \mathbb{N}}$ be a sequence of gradient discretizations with

$$ \mathcal{D}^l = (X_{\mathcal{D}^l}^0, \Pi_{\mathcal{D}_m^l}, \Pi_{\mathcal{D}_f^l}, \nabla_{\mathcal{D}_m^l}, \nabla_{\mathcal{D}_f^l}), $$

and let for all $l \in \mathbb{N}$, $\tilde{\Pi}_{\mathcal{D}_m^l}$, $\tilde{\Pi}_{\mathcal{D}_f^l}$ be a couple of linear mappings from $X_{\mathcal{D}^l}^0$ to $L^2(\Omega)$ and $L^2(\Gamma)$ respectively such that there exists a real sequence $(\xi_{\mathcal{D}^l})_{l \in \mathbb{N}}$ satisfying $\lim_{l \to \infty} \xi_{\mathcal{D}^l} = 0$ and

$$ \|\Pi_{\mathcal{D}_m^l} v_{\mathcal{D}^l} - \tilde{\Pi}_{\mathcal{D}_m^l} v_{\mathcal{D}^l}\|_{L^2(\Omega)} + \|\Pi_{\mathcal{D}_f^l} v_{\mathcal{D}^l} - \tilde{\Pi}_{\mathcal{D}_f^l} v_{\mathcal{D}^l}\|_{L^2(\Gamma)} \le \xi_{\mathcal{D}^l} \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l}. $$

for all $v_{\mathcal{D}^l} \in X_{\mathcal{D}^l}^0$ and all $l \in \mathbb{N}$. Let us define the following new sequence of gradient discretizations $(\tilde{\mathcal{D}}^l)_{l \in \mathbb{N}}$ with $\tilde{\mathcal{D}}^l = (X_{\tilde{\mathcal{D}}^l}^0, \tilde{\Pi}_{\tilde{\mathcal{D}}_m^l}, \tilde{\Pi}_{\tilde{\mathcal{D}}_f^l}, \nabla_{\tilde{\mathcal{D}}_m^l}, \nabla_{\tilde{\mathcal{D}}_f^l})$. Then, each property (coercivity or consistency or limit conformity or compactness) for the sequence $(\mathcal{D}^l)_{l \in \mathbb{N}}$ is equivalent to the same property for the sequence $(\tilde{\mathcal{D}}^l)_{l \in \mathbb{N}}$.

**Proof:** By symmetry it suffices to show that each property for the sequence $(\mathcal{D}^l)_{l \in \mathbb{N}}$ implies the same property for the sequence $(\tilde{\mathcal{D}}^l)_{l \in \mathbb{N}}$. Assuming the coercivity of $(\mathcal{D}^l)_{l \in \mathbb{N}}$, the coercivity property of the sequence $(\tilde{\mathcal{D}}^l)_{l \in \mathbb{N}}$ derives from $\|v_{\mathcal{D}^l}\|_{\mathcal{D}^l} = \|v_{\tilde{\mathcal{D}}^l}\|_{\tilde{\mathcal{D}}^l}$ for all $v_{\mathcal{D}^l} \in X_{\mathcal{D}^l}^0$ and from the estimate

$$ \|\tilde{\Pi}_{\mathcal{D}_m^l} v_{\mathcal{D}^l}\|_{L^2(\Omega)} + \|\tilde{\Pi}_{\mathcal{D}_f^l} v_{\mathcal{D}^l}\|_{L^2(\Gamma)} \le (C_{\mathcal{D}^l} + \xi_{\mathcal{D}^l}) \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l}. $$

Let $u \in V^0$, for all $v_{\mathcal{D}^l} \in X_{\mathcal{D}^l}^0$ one has the estimates

$$ \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l} \le \|\nabla u\|_{L^2(\Omega)^d} + \|\nabla_\tau \gamma u\|_{L^2(\Gamma)^{d-1}} + S_{\mathcal{D}^l}(u, v_{\mathcal{D}^l}), $$

and

$$ S_{\tilde{\mathcal{D}}^l}(u, v_{\mathcal{D}^l}) \le \xi_{\mathcal{D}^l} \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l} + S_{\mathcal{D}^l}(u, v_{\mathcal{D}^l}). $$

We deduce that

$$ S_{\tilde{\mathcal{D}}^l}(u, v_{\mathcal{D}^l}) \le \xi_{\mathcal{D}^l} (\|\nabla u\|_{L^2(\Omega)^d} + \|\nabla_\tau \gamma u\|_{L^2(\Gamma)^{d-1}}) + (1 + \xi_{\mathcal{D}^l}) S_{\mathcal{D}^l}(u, v_{\mathcal{D}^l}) $$

and hence the consistency of the sequence $\tilde{\mathcal{D}}^l$, $l \in \mathbb{N}$ derives from the consistency of the sequence $(\mathcal{D}^l)_{l \in \mathbb{N}}$.

For all $(q_m, q_f) \in H_\Sigma(\Omega, \Gamma)$ and all $v_{\mathcal{D}^l} \in X_{\mathcal{D}^l}^0$, it follows from (9) that

$$ W_{\tilde{\mathcal{D}}^l}(q_m, q_f, v_{\mathcal{D}^l}) \le W_{\mathcal{D}^l}(q_m, q_f, v_{\mathcal{D}^l}) + \xi_{\mathcal{D}^l} \|(q_m, q_f)\|_H \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l}, $$

from which we deduce that the limit conformity of the sequence $(\tilde{\mathcal{D}}^l)_{l \in \mathbb{N}}$ derives from the limit conformity of the sequence $(\mathcal{D}^l)_{l \in \mathbb{N}}$.

Finally, using the following estimates

$$ \|\tilde{\Pi}_{\mathcal{D}_m^l} v_{\mathcal{D}^l} - u_m\|_{L^2(\Omega)} \le \|\Pi_{\mathcal{D}_m^l} v_{\mathcal{D}^l} - u_m\|_{L^2(\Omega)} + \xi_{\mathcal{D}^l} \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l}, $$

$$ \|\tilde{\Pi}_{\mathcal{D}_f^l} v_{\mathcal{D}^l} - u_f\|_{L^2(\Gamma)} \le \|\Pi_{\mathcal{D}_f^l} v_{\mathcal{D}^l} - u_f\|_{L^2(\Gamma)} + \xi_{\mathcal{D}^l} \|v_{\mathcal{D}^l}\|_{\mathcal{D}^l}, $$

it is clear that the compactness of the sequence of gradient discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ implies the compactness of the sequence of gradient discretizations $(\tilde{\mathcal{D}}^l)_{l \in \mathbb{N}}$. $\square$
## 3.2 Application to the hybrid dimensional model

The gradient discretization of the hybrid dimensional model (4) is based on the primal variational formulation (5). It is defined by: find $u_D \in X_D^0$ such that for all $v_D \in X_D^0$:

$$ \int_{\Omega} \Lambda_m(\mathbf{x}) \nabla_{\mathcal{D}_m} u_D(\mathbf{x}) \cdot \nabla_{\mathcal{D}_m} v_D(\mathbf{x}) d\mathbf{x} + \int_{\Gamma} \Lambda_f(\mathbf{x}) \nabla_{\mathcal{D}_f} u_D(\mathbf{x}) \cdot \nabla_{\mathcal{D}_f} v_D(\mathbf{x}) d\tau_f(\mathbf{x}) \\ - \int_{\Omega} h_m(\mathbf{x}) \Pi_{\mathcal{D}_m} v_D(\mathbf{x}) d\mathbf{x} - \int_{\Gamma} h_f(\mathbf{x}) \Pi_{\mathcal{D}_f} v_D(\mathbf{x}) d\tau_f(\mathbf{x}) = 0. \quad (11) $$

**Proposition 3.2** Let $\mathcal{D}$ be a gradient discretization of (5). Then (11) has a unique solution $u_D \in X_D^0$ satisfying the a priori estimate

$$ \|u_D\|_{\mathcal{D}} \le \frac{C_{\mathcal{D}}}{\min(\underline{\lambda}_m, \overline{\lambda}_f \underline{d}_f)} \left( \|h_m\|_{L^2(\Omega)} + \|d_f h_f\|_{L^2(\Gamma)} \right). $$

**Proof:** For any solution $u_D \in X_D^0$ of (11), setting $v_D = u_D$ in (11), and using the Cauchy Schwarz inequality, the definition (6) of $C_{\mathcal{D}}$, and the assumption that $\|\cdot\|_{\mathcal{D}}$ defines a norm on $X_D^0$, we obtain the estimates

$$ \begin{aligned} \min(\underline{\lambda}_m, \overline{\lambda}_f \underline{d}_f) \|u_D\|_{\mathcal{D}}^2 &\le \int_{\Omega} h_m(\mathbf{x}) \Pi_{\mathcal{D}_m} u_D(\mathbf{x}) d\mathbf{x} + \int_{\Gamma} h_f(\mathbf{x}) \Pi_{\mathcal{D}_f} u_D(\mathbf{x}) d\tau_f(\mathbf{x}) \\ &\le C_{\mathcal{D}} \left( \|h_m\|_{L^2(\Omega)} + \|d_f h_f\|_{L^2(\Gamma)} \right) \|u_D\|_{\mathcal{D}}, \end{aligned} $$

from which we deduce the a priori estimate and hence the uniqueness and existence of a solution.

<figure><img src="image_3.png" /></figure>

**Proposition 3.3** Error estimates. Let $u \in V^0$ be the solution of (5) and let us set $(\mathbf{q}_m, \mathbf{q}_f) = (-\Lambda_m \nabla u, -d_f \Lambda_f \nabla_{\tau} \gamma u) \in H_{\Sigma}(\Omega, \Gamma)$. Let $\mathcal{D}$ be a gradient discretization of (5), and let $u_D \in X_D^0$ be the solution of (11). Then, there exist $C_1, C_2$ depending only on $\bar{\lambda}_m, \underline{\lambda}_m, \bar{\lambda}_f, \underline{\lambda}_f, \bar{d}_f, \underline{d}_f$, and $C_3, C_4$ depending only on $C_{\mathcal{D}}, \bar{\lambda}_m, \underline{\lambda}_m, \bar{\lambda}_f, \underline{\lambda}_f, \bar{d}_f, \underline{d}_f$, such that one has the following error estimates:

$$ \left\{ \begin{array}{l} \|\nabla u - \nabla_{\mathcal{D}_m} u_D\|_{L^2(\Omega)^d} + \|\nabla_{\tau} \gamma u - \nabla_{\mathcal{D}_f} u_D\|_{L^2(\Gamma)^{d-1}} \\ \le C_1 \mathcal{S}_{\mathcal{D}}(u) + C_2 \mathcal{W}_{\mathcal{D}}(\mathbf{q}_m, \mathbf{q}_f), \\ \\ \|\Pi_{\mathcal{D}_m} u_D - u\|_{L^2(\Omega)} + \|\Pi_{\mathcal{D}_f} u_D - \gamma u\|_{L^2(\Gamma)} \le C_3 \mathcal{S}_{\mathcal{D}}(u) + C_4 \mathcal{W}_{\mathcal{D}}(\mathbf{q}_m, \mathbf{q}_f). \end{array} \right. $$

**Proof:** Using the definition of $\mathcal{W}_D$ and the definitions of the solution $u_D$ of (11) and of the solution $u, (\mathbf{q}_m, \mathbf{q}_f)$ of (4), we obtain that for all $v_D \in X_D^0$

$$ \left| \int_{\Omega} \left( \Lambda_m \nabla_{\mathcal{D}_m} v_D \cdot (\nabla u - \nabla_{\mathcal{D}_m} u_D) \right) d\mathbf{x} + \int_{\Gamma} \left( \Lambda_f \nabla_{\mathcal{D}_f} v_D \cdot (\nabla_{\tau} \gamma u - \nabla_{\mathcal{D}_f} u_D) \right) d\tau_f(\mathbf{x}) \right| \le \|v_D\|_{\mathcal{D}} \mathcal{W}_D(\mathbf{q}_m, \mathbf{q}_f). $$

Let us introduce $w_D \in X_D^0$ defined as

$$ w_D = \underset{v_D \in X_D^0}{\operatorname{argmin}} S_D(u, v_D), $$

and let us set in the previous estimate $v_D = w_D - u_D$. Applying the Cauchy Schwarz inequality, we obtain the first estimate. In addition, from the definition of $C_{\mathcal{D}}$, we have that

$$ \|\Pi_{\mathcal{D}_m} w_D - \Pi_{\mathcal{D}_m} u_D\|_{L^2(\Omega)} + \|\Pi_{\mathcal{D}_f} w_D - \Pi_{\mathcal{D}_f} u_D\|_{L^2(\Gamma)} \le C_{\mathcal{D}} \|w_D - u_D\|_{\mathcal{D}}, $$

which proves the second estimate using the definition of $w_D$. $\square$
# 4 Two examples of Gradient Discretizations

Following [14], we consider generalised polyhedral meshes of $\Omega$ which allow for non planar faces. Let us stress that this general definition is used for the VAG scheme introduced in subsection 4.1 while the HFV scheme described in subsection 4.2 requires in addition that the faces are planar and that the face center $\mathbf{x}_\sigma$ is the center of gravity of the face $\sigma$.

**Definition 4.1 (Polyhedral mesh)** Let $\mathcal{M}$ be the set of cells that are disjoint open polyhedral subsets of $\Omega$ such that $\bigcup_{K \in \mathcal{M}} \bar{K} = \bar{\Omega}$. For all $K \in \mathcal{M}$, $\mathbf{x}_K$ denotes the so-called “center” of the cell $K$ under the assumption that $K$ is star-shaped with respect to $\mathbf{x}_K$. We then denote by $\mathcal{F}_K$ the set of interfaces of non zero $d-1$ dimensional measure among the interior faces $\bar{K} \cap \bar{L}$, $L \in \mathcal{M} \setminus \{K\}$, and the boundary interface $\bar{K} \cap \partial\Omega$, which possibly splits in several boundary faces. Let us denote by

$$ \mathcal{F} = \bigcup_{K \in \mathcal{M}} \mathcal{F}_K $$

the set of all faces of the mesh. Remark that the faces are not assumed to be planar, hence the term “generalised polyhedral mesh”. For $\sigma \in \mathcal{F}$, let $\mathcal{E}_\sigma$ be the set of interfaces of non zero $d-2$ dimensional measure among the interfaces $\bar{\sigma} \cap \bar{\sigma}'$, $\sigma' \in \mathcal{F} \setminus \{\sigma\}$. Then, we denote by

$$ \mathcal{E} = \bigcup_{\sigma \in \mathcal{F}} \mathcal{E}_\sigma $$

the set of all edges of the mesh. Let $\mathcal{V}_\sigma = \bigcup_{e,e' \in \mathcal{E}_\sigma, e \neq e'} (e \cap e')$ be the set of vertices of $\sigma$. For each $K \in \mathcal{M}$ we define $\mathcal{V}_K = \bigcup_{\sigma \in \mathcal{F}_K} \mathcal{V}_\sigma$, and we also denote by

$$ \mathcal{V} = \bigcup_{K \in \mathcal{M}} \mathcal{V}_K $$

the set of all vertices of the mesh. It is then assumed that for each face $\sigma \in \mathcal{F}$, there exists a so-called “center” of the face $\mathbf{x}_\sigma \in \sigma \setminus \bigcup_{e \in \mathcal{E}_\sigma} e$ such that $\mathbf{x}_\sigma = \sum_{s \in \mathcal{V}_\sigma} \beta_{\sigma,s} \mathbf{x}_s$, with $\sum_{s \in \mathcal{V}_\sigma} \beta_{\sigma,s} = 1$, and $\beta_{\sigma,s} \ge 0$ for all $\mathbf{s} \in \mathcal{V}_\sigma$; moreover the face $\sigma$ is assumed to be defined by the union of the triangles $T_{\sigma,e}$ defined by the face center $\mathbf{x}_\sigma$ and each edge $e \in \mathcal{E}_\sigma$.

The mesh is also supposed to be conforming w.r.t. the fracture network $\Gamma$ in the sense that for each $i \in I$ there exists a subset $\mathcal{F}_{\Gamma_i}$ of $\mathcal{F}$ such that $\bar{\Gamma}_i = \bigcup_{\sigma \in \mathcal{F}_{\Gamma_i}} \bar{\sigma}$. We will denote by $\mathcal{F}_\Gamma$ the set of fracture faces $\bigcup_{i \in I} \mathcal{F}_{\Gamma_i}$, and by $\mathcal{E}_\Gamma$ the set of fracture edges $\bigcup_{\sigma \in \mathcal{F}_\Gamma} \mathcal{E}_\sigma$.

A tetrahedral sub-mesh of $\mathcal{M}$ is defined by

$$ \mathcal{T} = \{T_{K,\sigma,e}, e \in \mathcal{E}_\sigma, \sigma \in \mathcal{F}_K, K \in \mathcal{M}\}, $$

where $T_{K,\sigma,e}$ is the tetrahedron joining the cell center $\mathbf{x}_K$ to the triangle $T_{\sigma,e}$ (see Figure 2 for examples of such tetrahedra).

Let $\rho_T$ denote the insphere diameter of a given tetrahedron $T$, $h_T$ its diameter, and $h_T = \max_{T \in \mathcal{T}} h_T$. We will assume in the convergence analysis that the family of tetrahedral sub-meshes $\mathcal{T}$ is shape regular. Hence let us set

$$ \theta_T = \max_{T \in \mathcal{T}} \frac{h_T}{\rho_T}. $$
## 4.1 Vertex Approximate Gradient Discretization

The VAG discretization has been introduced in [14] for diffusive problems on heterogeneous anisotropic media. Its extension to the hybrid dimensional Darcy model is based on the following vector space of degrees of freedom:

$$X_D = \{v_K, v_s, v_\sigma \in \mathbb{R}, K \in \mathcal{M}, s \in \mathcal{V}, \sigma \in \mathcal{F}_\Gamma\}, \quad (12)$$

and its subspace with homogeneous Dirichlet boundary conditions on $\partial\Omega$:

$$X_D^0 = \{v_D \in X_D | v_s = 0 \text{ for } s \in \mathcal{V}_{\text{ext}}\}. \quad (13)$$

where $\mathcal{V}_{\text{ext}} = \mathcal{V} \cap \partial\Omega$ denotes the set of boundary vertices, and $\mathcal{V}_{\text{int}} = \mathcal{V} \setminus \partial\Omega$ denotes the set of interior vertices.

<figure><img src="image_4.png" /><figcaption>Figure 2: Degrees of freedom of the VAG scheme: cell unknowns $v_K, v_L$, fracture face unknown $v_\sigma$, and node unknowns $v_s, v_{s1}, v_{s2}, v_{s3}, v_{s4}$. The fracture faces of $\mathcal{F}_\Gamma$ are in bold. The value of $v_{\sigma'}$ is obtained by interpolation of the node unknowns $v_{s1}, v_{s2}, v_{s3}, v_{s4}$ of the face $\sigma' \in \mathcal{F} \setminus \mathcal{F}_\Gamma$ while $v_\sigma$ is kept as an unknown for $\sigma \in \mathcal{F}_\Gamma$.</figcaption></figure>

A $\mathbb{P}_1$ finite element discretization of $V$ is built using the tetrahedral sub-mesh $\mathcal{T}$ of $\mathcal{M}$ and a second order interpolation at the face centers $\mathbf{x}_\sigma$, $\sigma \in \mathcal{F} \setminus \mathcal{F}_\Gamma$ defined by the operator $I_\sigma: X_D \to \mathbb{R}$ such that

$$I_\sigma(v_D) = \sum_{s \in \mathcal{V}_\sigma} \beta_{\sigma,s} v_s.$$

For a given $v_D \in X_D$, we define the function $\Pi_\tau v_D \in V$ as the continuous piecewise affine function on each tetrahedron of $\mathcal{T}$ such that $\Pi_\tau v_D(\mathbf{x}_K) = v_K$, $\Pi_\tau v_D(\mathbf{x}_s) = v_s$, $\Pi_\tau v_D(\mathbf{x}_\sigma) = v_\sigma$, and $\Pi_\tau v_D(\mathbf{x}_{\sigma'}) = I_{\sigma'}(v_D)$ for all $K \in \mathcal{M}, s \in \mathcal{V}, \sigma \in \mathcal{F}_\Gamma$, and $\sigma' \in \mathcal{F} \setminus \mathcal{F}_\Gamma$. The discrete gradient operators are obtained from this finite element discretization of $V$, setting

$$\nabla_{\mathcal{D}_m} = \nabla \Pi_\tau \text{ and } \nabla_{\mathcal{D}_f} = \nabla_\tau \gamma \Pi_\tau. \quad (14)$$

In addition to this conforming finite element discretization of $V$, the VAG discretization uses two non conforming piecewise constant reconstructions of functions from $X_D$ into respectively $L^2(\Omega)$ and $L^2(\Gamma)$ based on a partition of the cells and of the fracture faces. These partitions are respectively denoted, for all $K \in \mathcal{M}$, by

$$\bar{K} = \bar{\omega}_K \cup \left( \bigcup_{s \in \mathcal{V}_K \cap \mathcal{V}_{\text{int}}} \bar{\omega}_{K,s} \right) \cup \left( \bigcup_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_\Gamma} \bar{\omega}_{K,\sigma} \right),$$
and, for all $\sigma \in \mathcal{F}_\Gamma$, by

$$ \bar{\sigma} = \bar{\omega}_\sigma \cup \left( \bigcup_{s \in \mathcal{V}_\sigma \cap \mathcal{V}_{\text{int}}} \bar{\omega}_{\sigma,s} \right). $$

Then, the function reconstruction operators are defined by

$$ \Pi_{\mathcal{D}_m} v_{\mathcal{D}}(\mathbf{x}) = \begin{cases} v_K & \text{for all } \mathbf{x} \in \omega_K, K \in \mathcal{M}, \\ v_s & \text{for all } \mathbf{x} \in \omega_{K,s}, \mathbf{s} \in \mathcal{V}_K \cap \mathcal{V}_{\text{int}}, K \in \mathcal{M}, \\ v_\sigma & \text{for all } \mathbf{x} \in \omega_{K,\sigma}, \sigma \in \mathcal{F}_K \cap \mathcal{F}_\Gamma, K \in \mathcal{M}, \end{cases} \quad (15) $$

and

$$ \Pi_{\mathcal{D}_f} v_{\mathcal{D}}(\mathbf{x}) = \begin{cases} v_\sigma & \text{for all } \mathbf{x} \in \omega_\sigma, \sigma \in \mathcal{F}_\Gamma, \\ v_s & \text{for all } \mathbf{x} \in \omega_{\sigma,s}, \mathbf{s} \in \mathcal{V}_\sigma \cap \mathcal{V}_{\text{int}}, \sigma \in \mathcal{F}_\Gamma. \end{cases} \quad (16) $$

It is shown below that the above VAG discretization defines a coercive, consistent, limit conforming and compact gradient discretization whatever the choice of these partitions.

**Properties of VAG discretization:** we state without proof two results that can be readily adapted from [5] Lemmas 3.4 and 3.7 noticing that the shape regularity of $\mathcal{T}$ implies the shape regularity of the triangular submesh of $\Gamma$ defined by $\mathcal{T} \cap \Gamma$.

**Lemma 4.1** There exists $C > 0$ depending only on $\theta_{\mathcal{T}}$ such that, for all $v_{\mathcal{D}} \in X_{\mathcal{D}}$, one has the estimates

$$ \begin{aligned} \|\Pi_{\mathcal{D}_m} v_{\mathcal{D}} - \Pi_{\mathcal{T}} v_{\mathcal{D}}\|_{L^2(\Omega)} &\leq C h_{\mathcal{T}} \|\nabla_{\mathcal{D}_m} v_{\mathcal{D}}\|_{L^2(\Omega)^d}, \\ \|\Pi_{\mathcal{D}_f} v_{\mathcal{D}} - \gamma \Pi_{\mathcal{T}} v_{\mathcal{D}}\|_{L^2(\Gamma)} &\leq C h_{\mathcal{T}} \|\nabla_{\mathcal{D}_f} v_{\mathcal{D}}\|_{L^2(\Gamma)^{d-1}}. \end{aligned} $$

For any continuous function $\varphi \in C^0(\bar{\Omega})$, let us introduce the operator $P_{\mathcal{D}} : C^0(\bar{\Omega}) \to X_{\mathcal{D}}$ such that

$$ (P_{\mathcal{D}}\varphi)_K = \varphi(\mathbf{x}_K), (P_{\mathcal{D}}\varphi)_s = \varphi(\mathbf{x}_s), (P_{\mathcal{D}}\varphi)_\sigma = \varphi(\mathbf{x}_\sigma) $$

for all $K \in \mathcal{M}$, $s \in \mathcal{V}$ and $\sigma \in \mathcal{F}_\Gamma$.

We have the following classical finite element approximation result for the finite element interpolation operator $\Pi_{\mathcal{T}} P_{\mathcal{D}}$.

**Proposition 4.1** For all $\varphi \in C^\infty(\bar{\Omega})$, there exists $C_\varphi > 0$ depending only on $\varphi$, $\theta_{\mathcal{T}}$ such that one has the error estimates

$$ \|\varphi - \Pi_{\mathcal{T}} P_{\mathcal{D}} \varphi\|_{L^2(\Omega)} + \|\gamma \varphi - \gamma \Pi_{\mathcal{T}} P_{\mathcal{D}} \varphi\|_{L^2(\Gamma)} \leq h_{\mathcal{T}}^2 C_\varphi, $$

and

$$ \|\nabla \varphi - \nabla \Pi_{\mathcal{T}} P_{\mathcal{D}} \varphi\|_{L^2(\Omega)^d} + \|\nabla \gamma \varphi - \nabla_{\mathcal{T}} \gamma \Pi_{\mathcal{T}} P_{\mathcal{D}} \varphi\|_{L^2(\Gamma)^{d-1}} \leq h_{\mathcal{T}} C_\varphi. $$

Let us now state our main result concerning the VAG discretization.

**Proposition 4.2 (Main result on VAG)** Let us consider a family of meshes $(\mathcal{M}^l)_{l \in \mathbb{N}}$ as defined above. It is assumed that the family of tetrahedral submeshes $\mathcal{T}^l$ of $\mathcal{M}^l$ is shape regular in the sense that there exists $\theta > 0$ such that $\theta_{\mathcal{T}^l} \leq \theta$ for all $l \in \mathbb{N}$. It is also assumed that $h_{\mathcal{T}^l}$ tends to zero when $l \to +\infty$. Then, the sequence of VAG discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ with $\mathcal{D}^l = (X_{\mathcal{D}^l}^0, \Pi_{\mathcal{D}^l}, \Pi_{\mathcal{D}^l}, \nabla_{\mathcal{D}^l}, \nabla_{\mathcal{D}^l})$ defined by (13), (15), (16), (14) is coercive, consistent, limit conforming and compact.
**Proof:** Let us denote by

$$ \mathcal{D}_{FE} = (X_{\mathcal{D}}^{0}, \Pi_{\tau}, \gamma\Pi_{\tau}, \nabla_{\mathcal{D}_{m}}, \nabla_{\mathcal{D}_{f}}), $$

the conforming $\mathbb{P}_1$ finite element VAG discretization. It results from Lemma 2.1 that

$$ \|\Pi_{\tau} v_{\mathcal{D}}\|_{L^2(\Omega)} + \|\gamma\Pi_{\tau} v_{\mathcal{D}}\|_{L^2(\Gamma)} \leq C_P \|\nabla_{\mathcal{D}_m} v_{\mathcal{D}}\|_{L^2(\Omega)^d} \quad (17) $$

for all $v_{\mathcal{D}} \in X_{\mathcal{D}}^0$. On the other hand for all $(\mathbf{q}_m, \mathbf{q}_f) \in H_{\Sigma}(\Omega, \Gamma)$ and all $v_{\mathcal{D}} \in X_{\mathcal{D}}^0$ one has

$$ W_{\mathcal{D}_{FE}}(\mathbf{q}_m, \mathbf{q}_f, v_{\mathcal{D}}) = 0. \quad (18) $$

We deduce from (17) and (18) that the sequence of conforming gradient discretizations $(\mathcal{D}_{FE}^l)_{l \in \mathbb{N}}$ is coercive and limit conforming. The consistency of $(\mathcal{D}_{FE}^l)_{l \in \mathbb{N}}$ results from Proposition 4.1 and from the density of $C_c^{\infty}(\Omega)$ in $V^0$ given by Proposition 2.2. The following estimates

$$ \|\Pi_{\tau} v_{\mathcal{D}}\|_{H^1(\Omega)} \leq C_1 \|\nabla \Pi_{\tau} v_{\mathcal{D}}\|_{L^2(\Omega)^d}, $$

and

$$ \|\gamma_i \Pi_{\tau} v_{\mathcal{D}}\|_{H^1(\Gamma_i)} \leq C_2 \left( \|\nabla \Pi_{\tau} v_{\mathcal{D}}\|_{L^2(\Omega)^d} + \|\nabla_{\tau} \gamma_i \Pi_{\tau} v_{\mathcal{D}}\|_{L^2(\Gamma)^{d-1}} \right) $$

for constants $C_1$ and $C_2$ independent on the mesh and on $v_{\mathcal{D}} \in X_{\mathcal{D}}^0$ are deduced from the Poincaré inequality and the Trace theorem. Then, thanks to the Rellich Compactness Theorem, one obtains the compactness of $(\mathcal{D}_{FE}^l)_{l \in \mathbb{N}}$. From Lemma 3.1 and Lemma 4.1 we deduce that the sequence $(\mathcal{D}_{FE}^l)_{l \in \mathbb{N}}$ is also coercive, consistent, limit conforming and compact. $\square$

## 4.2 Hybrid Finite Volume Discretization

In this subsection, the HFV scheme introduced in [13] is extended to the hybrid dimensional Darcy flow model. Let us recall that the HFV scheme of [13] has been generalised in [12] as the family of Hybrid Mimetic Mixed methods which encompasses the family of Mimetic Finite Difference schemes [8]. In the following, we focus on the particular case presented in [13] for the sake of simplicity.

Let us recall that, for the HFV scheme, the faces $\sigma \in \mathcal{F}$ are assumed to be planar and $\mathbf{x}_\sigma$ is assumed to be the center of gravity of the face $\sigma$. We also denote by $\mathbf{x}_e$ the center of the edge $e \in \mathcal{E}$. Let $\mathcal{F}_{int} \subset \mathcal{F}$ (resp. $\mathcal{E}_{int} \subset \mathcal{E}$) denote the subset of interior faces (resp. interior edges). The vector space of degrees of freedom $X_{\mathcal{D}}$ is defined by

$$ X_{\mathcal{D}} = \{u_K \in \mathbb{R}, u_{\sigma} \in \mathbb{R}, u_e \in \mathbb{R} \text{ for all } K \in \mathcal{M}, \sigma \in \mathcal{F}, e \in \mathcal{E}_{\Gamma}\}, \quad (19) $$

and its subspace $X_{\mathcal{D}}^0$ is defined by

$$ X_{\mathcal{D}}^0 = \{u_{\mathcal{D}} \in X_{\mathcal{D}} | u_{\sigma} = 0, u_e = 0 \text{ for all } \sigma \in \mathcal{F} \setminus \mathcal{F}_{\text{int}} \text{ and } e \in \mathcal{E}_{\Gamma} \setminus \mathcal{E}_{\text{int}}\}. \quad (20) $$

For any continuous function $\varphi \in C^0(\bar{\Omega})$, let us define its projection $P_{\mathcal{D}}\varphi$ onto $X_{\mathcal{D}}$ such that $(P_{\mathcal{D}}\varphi)_K = \varphi(\mathbf{x}_K)$, $(P_{\mathcal{D}}\varphi)_\sigma = \varphi(\mathbf{x}_\sigma)$, $(P_{\mathcal{D}}\varphi)_e = \varphi(\mathbf{x}_e)$ for $K \in \mathcal{M}, \sigma \in \mathcal{F}, e \in \mathcal{E}_{\Gamma}$.
<figure><img src="image_5.png" /><figcaption>Figure 3: Degrees of freedom of the HFV scheme: cell unknowns $u_K$, $u_L$, fracture face unknown $u_\sigma$, matrix face unknown $u_{\sigma'}$ and fracture edge unknown $u_e$.</figcaption></figure>

For each cell $K$ and $u_D \in X_D$, let us define

$$ \nabla_K u_D = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K} |\sigma|(u_\sigma - u_K)\mathbf{n}_{K,\sigma}, $$

where $|K|$ is the volume of the cell $K$, $|\sigma|$ is the surface of the face $\sigma$, and $\mathbf{n}_{K,\sigma}$ is the unit normal vector of the face $\sigma \in \mathcal{F}_K$ outward to the cell $K$. We recall from [13] that $\nabla_K u_D$ is exact on affine functions $\varphi$ in the sense that $\nabla_K P_D \varphi = \nabla \varphi$. Also note that $\nabla_K u_D$ does not depend on $u_K$ since $\sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K,\sigma} = 0$. Hence a stabilised discrete gradient is defined as follows:

$$ \nabla_{K,\sigma} u_D = \nabla_K u_D + R_{K,\sigma}(u_D) \mathbf{n}_{K,\sigma}, \sigma \in \mathcal{F}_K, $$

with

$$ R_{K,\sigma}(u_D) = \frac{\sqrt{d}}{d_{K,\sigma}} \left( u_{\sigma} - u_K - \nabla_K u_D \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_K) \right), $$

setting $d_{K,\sigma} = \mathbf{n}_{K,\sigma} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_K)$ which leads to the definition of the matrix discrete gradient

$$ \nabla_{D_m} u_D(\mathbf{x}) = \nabla_{K,\sigma} u_D \text{ on } K_{\sigma} \text{ for all } K \in \mathcal{M}, \sigma \in \mathcal{F}_K, \quad (21) $$

where $K_\sigma$ is the cone joining the face $\sigma$ to the cell center $\mathbf{x}_K$. The fracture discrete gradient is defined similarly by

$$ \nabla_{D_f} u_D(\mathbf{x}) = \nabla_{\sigma,e} u_D \text{ on } \sigma_e \text{ for all } \sigma \in \mathcal{F}_\Gamma, e \in \mathcal{E}_\sigma, \quad (22) $$

with

$$ \nabla_{\sigma,e} u_D = \nabla_{\sigma} u_D + R_{\sigma,e}(u_D) \mathbf{n}_{\sigma,e}, $$

and

$$ \nabla_{\sigma} u_D = \frac{1}{|\sigma|} \sum_{e \in \mathcal{E}_{\sigma}} |e|(u_e - u_{\sigma})\mathbf{n}_{\sigma,e}, $$

$$ R_{\sigma,e}(u_D) = \frac{\sqrt{d-1}}{d_{\sigma,e}} \left( u_e - u_{\sigma} - \nabla_{\sigma} u_D \cdot (\mathbf{x}_e - \mathbf{x}_{\sigma}) \right), $$

where $\mathbf{n}_{\sigma,e}$ is the unit normal vector to the edge $e$ in the tangent plane of the face $\sigma$ and outward to the face $\sigma$, $d_{\sigma,e} = \mathbf{n}_{\sigma,e} \cdot (\mathbf{x}_e - \mathbf{x}_{\sigma})$, and $\sigma_e$ is the triangle of base $e$ and vertex $\mathbf{x}_{\sigma}$.

The matrix and fracture discrete gradients are both consistent in the sense that for any affine function $\varphi \in C^0(\bar{\Omega})$ one has $\nabla_{D_m} P_D \varphi = \nabla \varphi$, and for any function $\varphi \in C^0(\bar{\Omega})$ affine on
the fracture $\Gamma_i$ one has $\nabla_{\mathcal{D}_f} P_{\mathcal{D}}\varphi = \nabla_{\tau_i}\varphi$ on $\Gamma_i$. We recall also from [13] that for all $u_{\mathcal{D}} \in X_{\mathcal{D}}$ one has

$$ \sum_{\sigma \in \mathcal{F}_K} |\sigma| d_{K,\sigma} R_{K,\sigma}(u_{\mathcal{D}}) \mathbf{n}_{K,\sigma} = 0 \quad \text{and} \quad \sum_{e \in \mathcal{E}_{\sigma}} |e| d_{\sigma,e} R_{\sigma,e}(u_{\mathcal{D}}) \mathbf{n}_{\sigma,e} = 0. \qquad (23) $$

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 \in \mathcal{M}$, by

$$ \bar{K} = \bar{\omega}_K \cup \left( \bigcup_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{int}} \bar{\omega}_{K,\sigma} \right), $$

and, for all $\sigma \in \mathcal{F}_{\Gamma}$, by

$$ \bar{\sigma} = \bar{\omega}_{\sigma} \cup \left( \bigcup_{e \in \mathcal{E}_{\sigma} \cap \mathcal{E}_{int}} \bar{\omega}_{\sigma,e} \right). $$

Then, the function reconstruction operators are defined for all $v_{\mathcal{D}} \in X_{\mathcal{D}}$ by

$$ \Pi_{\mathcal{D}_m} v_{\mathcal{D}}(\mathbf{x}) = \begin{cases} v_K & \text{for all } \mathbf{x} \in \omega_K, K \in \mathcal{M}, \\ v_{\sigma} & \text{for all } \mathbf{x} \in \omega_{K,\sigma}, \sigma \in \mathcal{F}_K \cap \mathcal{F}_{int}, K \in \mathcal{M}, \end{cases} \qquad (24) $$

and

$$ \Pi_{\mathcal{D}_f} v_{\mathcal{D}}(\mathbf{x}) = \begin{cases} v_{\sigma} & \text{for all } \mathbf{x} \in \omega_{\sigma}, \sigma \in \mathcal{F}_{\Gamma}, \\ v_e & \text{for all } \mathbf{x} \in \omega_{\sigma,e}, e \in \mathcal{E}_{\sigma} \cap \mathcal{E}_{int}, \sigma \in \mathcal{F}_{\Gamma}. \end{cases} \qquad (25) $$

As for the VAG scheme, it is shown below that the above HFV discretization defines a coercive, consistent, limit conforming and compact gradient discretization whatever the choice of these partitions.

Let us define the two piecewise constant mappings $\Pi_{\mathcal{M}}$ (resp. $\Pi_{\mathcal{F}}$) from $X_{\mathcal{D}}$ to $L^2(\Omega)$ (resp. $L^2(\Gamma)$) such that for all $v_{\mathcal{D}} \in X_{\mathcal{D}}$

$$ \Pi_{\mathcal{M}} v_{\mathcal{D}}|_{K} = v_{K} \quad \text{for all } K \in \mathcal{M} \text{ (resp } \Pi_{\mathcal{F}} v_{\mathcal{D}}|_{\sigma} = v_{\sigma} \text{ for all } \sigma \in \mathcal{F}_{\Gamma}). \qquad (26) $$

Following the proof of Lemma 4.1 from [13], one can obtain the following Lemma.

**Lemma 4.2** There exists $C > 0$ depending only on $\theta_{\mathcal{T}}$ such that for all $u_{\mathcal{D}} \in X_{\mathcal{D}}$ one has

$$ \|\Pi_{\mathcal{D}_m} u_{\mathcal{D}} - \Pi_{\mathcal{M}} u_{\mathcal{D}}\|_{L^2(\Omega)} + \|\Pi_{\mathcal{D}_f} u_{\mathcal{D}} - \Pi_{\mathcal{F}} u_{\mathcal{D}}\|_{L^2(\Gamma)} \leq C h_{\mathcal{T}} \|u_{\mathcal{D}}\|_{\mathcal{D}}. $$

**Properties of HFV discretizations:** Let us first consider the HFV discretization

$$ \mathcal{D} = (X_{\mathcal{D}}^0, \Pi_{\mathcal{M}}, \Pi_{\mathcal{F}}, \nabla_{\mathcal{D}_m}, \nabla_{\mathcal{D}_f}), $$

defined by the vector space $X_{\mathcal{D}}^0$ from (20), the discrete gradient operators $\nabla_{\mathcal{D}_m}$ from (21) and $\nabla_{\mathcal{D}_f}$ from (22), and the function reconstruction operators $\Pi_{\mathcal{M}}$, $\Pi_{\mathcal{F}}$ from (26). From Lemma 5.3 and Lemma 4.1 of [13] and Lemma 1.51 of [11], one has the following discrete Poincaré estimates

$$ \begin{aligned} \|\Pi_{\mathcal{M}} u_{\mathcal{D}}\|_{L^2(\Omega)} &\leq C_{\mathcal{D},m} \|\nabla_{\mathcal{D}_m} u_{\mathcal{D}}\|_{L^2(\Omega)^d}, \\ \|\Pi_{\mathcal{F}} u_{\mathcal{D}}\|_{L^2(\Gamma)} &\leq C_{\mathcal{D},f} \left( \|\nabla_{\mathcal{D}_m} u_{\mathcal{D}}\|_{L^2(\Omega)^d} + \|\nabla_{\mathcal{D}_f} u_{\mathcal{D}}\|_{L^2(\Gamma)^{d-1}} \right), \end{aligned} \qquad (27) $$

for all $u_{\mathcal{D}} \in X_{\mathcal{D}}^0$ with $C_{\mathcal{D},m}$ and $C_{\mathcal{D},f}$ depending only on $\theta_{\mathcal{T}}$.

It follows from Lemma 4.3 of [13] that for all $\varphi \in C^{\infty}(\bar{\Omega})$ there exists $C > 0$ depending only on $\theta_{\mathcal{T}}$ and $\varphi$ such that

$$ \|\nabla_{\mathcal{D}_m} P_{\mathcal{D}}\varphi - \nabla_{\varphi}\|_{L^2(\Omega)} + \|\nabla_{\mathcal{D}_f} P_{\mathcal{D}}\varphi - \nabla_{\varphi}\|_{L^2(\Gamma)} \leq C h_{\mathcal{T}} \qquad (28) $$
It is easy to show that for all $\varphi \in C^\infty(\bar{\Omega})$, there exists $C > 0$, depending only on $\varphi$ and $\theta_T$ such that

$$ \|\Pi_M P_D \varphi - \varphi\|_{L^2(\Omega)} + \|\Pi_F P_D \varphi - \varphi\|_{L^2(\Gamma)} \leq Ch_T. \quad (29) $$

**Proposition 4.3** Let $(\varphi_m, \varphi_f) \in C_{H_\Sigma}^\infty(\Omega, \Gamma)$, there exist $C$ depending only on $(\varphi_m, \varphi_f)$ and $\theta_T$ such that for all $u_D \in X_D^0$

$$ W_D(\varphi_m, \varphi_f, u_D) \leq Ch_T \|u_D\|_D. $$

**Proof:** Let us define $\varphi_K = \frac{1}{|K|} \int_K \varphi_m dx$ for all $K \in \mathcal{M}$, and $\varphi_{K,\sigma} = \lim_{\epsilon \to 0^+} \frac{1}{|\sigma|} \int_{\sigma} \varphi_m(x - n_{K,\sigma}\epsilon) d\tau(x)$ for all $\sigma \in \mathcal{F}_K$.

Let us define

$$ A_1^D = A_{11}^D + A_{12}^D = \int_{\Omega} \nabla_{D_m} u_D \cdot \varphi_m dx, $$

with

$$ A_{11}^D = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma|(u_\sigma - u_K)\varphi_K \cdot n_{K,\sigma}, $$

and

$$ A_{12}^D = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} R_{K,\sigma}(u_D)n_{K,\sigma} \cdot \int_{K_\sigma} \varphi_m dx. $$

Using (23), one has

$$ A_{12}^D = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} R_{K,\sigma}(u_D)n_{K,\sigma} \cdot \int_{K_\sigma} (\varphi_m - \varphi_K)dx. $$

We can deduce as in Lemma 4.2 from [13] that there exists $C$ depending only on $\varphi_m$, $\theta_T$ such that

$$ |A_{12}^D| \leq Ch_T \|\nabla_{D_m} u_D\|_{L^2(\Omega)^d}. \quad (30) $$

Let us consider the term $A_2^D$ defined by

$$ A_2^D = \sum_{\alpha \in \mathcal{A}} \int_{\Omega_\alpha} (\Pi_M u_D) \mathrm{div}(\varphi_{m,\alpha}) dx = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| u_K \varphi_{K,\sigma} \cdot n_{K,\sigma}. $$

Adding and subtracting $\sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| u_\sigma \varphi_{K,\sigma} \cdot n_{K,\sigma}$ to $A_2^D$ and using that $\sum_{K \in \mathcal{M}_\sigma} |\sigma| \varphi_{K,\sigma} \cdot n_{K,\sigma} = 0$ for all $\sigma \in \mathcal{F} \setminus \mathcal{F}_\Gamma$, leads to

$$ A_2^D = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma|(u_K - u_\sigma)\varphi_{K,\sigma} \cdot n_{K,\sigma} + \sum_{\sigma \in \mathcal{F}_\Gamma} \sum_{K \in \mathcal{M}_\sigma} |\sigma| u_\sigma \varphi_{K,\sigma} \cdot n_{K,\sigma}. $$

It results that

$$ \begin{aligned} & A_{11}^D + A_2^D - \sum_{i \in I} \int_{\Gamma_i} (\Pi_F u_D) \llbracket \varphi_m \cdot n_i \rrbracket d\tau(x) \\ & = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma|(u_K - u_\sigma)(\varphi_{K,\sigma} - \varphi_K) \cdot n_{K,\sigma} \end{aligned} \quad (31) $$
Therefore, applying Cauchy-Schwartz inequality to (31), using the regularity of $\varphi_m$, and the estimate (30), we deduce that there exists $C$ depending only on $\varphi_m$, $\theta_T$ such that

$$ \left| \sum_{\alpha \in \mathcal{A}} \int_{\Omega_{\alpha}} (\nabla_{\mathcal{D}_m} u_{\mathcal{D}} \cdot \varphi_{m,\alpha} + (\Pi_{\mathcal{M}} u_{\mathcal{D}}) \mathrm{div}(\varphi_{m,\alpha})) d\mathbf{x} - \sum_{i \in I} \int_{\Gamma_i} (\Pi_{\mathcal{F}} u_{\mathcal{D}}) [\varphi_m \cdot \mathbf{n}_i] d\tau(\mathbf{x}) \right| \leq C h_{\tau} \|\nabla_{\mathcal{D}_m} u_{\mathcal{D}}\|_{L^2(\Omega)^d}. \quad (32) $$

Next, we define for all $\sigma \in \mathcal{F}_{\Gamma}$ and $e \in \mathcal{E}_{\sigma}$

$$ \varphi_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \varphi_{f} d\tau(\mathbf{x}) \quad \text{and} \quad \varphi_{\sigma,e} = \lim_{\epsilon \to 0^+} \frac{1}{|e|} \int_{e} \varphi_{f}(\mathbf{x} - \mathbf{n}_{\sigma,e}\epsilon) dl(\mathbf{x}). $$

Let us set

$$ B_{1}^{\mathcal{D}} = B_{11}^{\mathcal{D}} + B_{12}^{\mathcal{D}} = \sum_{i \in I} \int_{\Gamma_i} \nabla_{\mathcal{D}_f} u_{\mathcal{D}} \cdot \varphi_{f,i} d\tau(\mathbf{x}), $$

with

$$ B_{11}^{\mathcal{D}} = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e|(u_e - u_{\sigma})\varphi_{\sigma} \cdot \mathbf{n}_{\sigma,e}, $$

and

$$ B_{12}^{\mathcal{D}} = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} R_{\sigma,e}(u_{\mathcal{D}}) \mathbf{n}_{\sigma,e} \cdot \int_{\sigma_e} \varphi_f d\tau(\mathbf{x}). $$

Using (23), one has

$$ B_{12}^{\mathcal{D}} = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} R_{\sigma,e}(u_{\mathcal{D}}) \mathbf{n}_{\sigma,e} \cdot \int_{\sigma_e} (\varphi_f - \varphi_{\sigma}) d\tau(\mathbf{x}). $$

We can deduce as in [13] that there exists $C$ depending only on $\varphi_f$, $\theta_T$ such that

$$ |B_{12}^{\mathcal{D}}| \leq C h_{\tau} \|\nabla_{\mathcal{D}_f} u_{\mathcal{D}}\|_{L^2(\Gamma)^{d-1}}. \quad (33) $$

Let us consider the term $B_2^D$ defined by

$$ B_2^D = \sum_{i \in I} \int_{\Gamma_i} (\Pi_{\mathcal{F}} u_{\mathcal{D}}) \mathrm{div}(\varphi_{f,i}) d\tau(\mathbf{x}) = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e| u_{\sigma} \varphi_{\sigma,e} \cdot \mathbf{n}_{\sigma,e}. $$

Adding and subtracting $\sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e| u_{\sigma} \varphi_{\sigma,e} \cdot \mathbf{n}_{\sigma,e}$ to $B_2^D$ we obtain that

$$ B_2^D = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e|(u_{\sigma} - u_e)\varphi_{\sigma,e} \cdot \mathbf{n}_{\sigma,e} + \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e| u_e \varphi_{\sigma,e} \cdot \mathbf{n}_{\sigma,e} $$

Taking into account the definition of $\varphi_{f,i}$ and the fact that $u_e = 0$ for all $e \in \mathcal{E}_{\Gamma} \setminus \mathcal{E}_{int}$ we conclude that

$$ \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e| u_e \varphi_{\sigma,e} \cdot \mathbf{n}_{\sigma,e} = 0. $$

It results that

$$ B_{11}^{\mathcal{D}} + B_{2}^{\mathcal{D}} = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \sum_{e \in \mathcal{E}_{\sigma}} |e|(u_{\sigma} - u_e)(\varphi_{\sigma,e} - \varphi_{\sigma}) \cdot \mathbf{n}_{\sigma,e}, $$
from which we can deduce as in [13] and using (33) that there exists C depending only on $\varphi_f$ and $\theta_T$ such that

$$ \left| \sum_{i \in I} \int_{\Gamma_i} \left( \nabla_{\mathcal{D}_f} u_{\mathcal{D}} \cdot \varphi_{f,i} + (\Pi_{\mathcal{F}} u_{\mathcal{D}}) \mathrm{div}(\varphi_{f,i}) \right) d\tau(\mathbf{x}) \right| \le C h_T \| \nabla_{\mathcal{D}_f} u_{\mathcal{D}} \|_{L^2(\Gamma)^{d-1}}. \quad (34) $$

Combining the estimates (32) and (34) concludes the proof of Proposition 4.3. $\square$

**Proposition 4.4 (Main result on HFV)** Let us consider a family of meshes $(\mathcal{M}^l)_{l \in \mathbb{N}}$ as defined above. It is assumed that the family of tetrahedral submeshes $\mathcal{T}^l$ of $\mathcal{M}^l$ is shape regular in the sense that there exist a positive constant $\theta$ such that $\theta_T^l < \theta$ for all $l \in \mathbb{N}$. It is also assumed that $h_T^l$ tends to zero when $l \to +\infty$. Then, the sequence of HFV discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ with $\mathcal{D}^l = (X_{\mathcal{D}^l}^0, \Pi_{\mathcal{D}^l}, \Pi_{\mathcal{D}^l}, \nabla_{\mathcal{D}^l}, \nabla_{\mathcal{D}^l})$ defined by (20), (24), (25), (21), (22) is coercive, consistent, limit conforming and compact.

**Proof:** The coercivity of the sequence of HFV discretizations

$$ (X_{\mathcal{D}^l}^0, \Pi_{\mathcal{M}^l}, \Pi_{\mathcal{F}^l}, \nabla_{\mathcal{D}^l}, \nabla_{\mathcal{D}^l})_{l \in \mathbb{N}}, $$

results from (27). Its consistency is obtained using (28), (29) and the density of $C_c^\infty(\Omega)$ in $V^0$ given by Proposition 2.2. Its limit conformity is obtained by Proposition 4.3 and the density of $C_{H_\Sigma}^\infty(\Omega, \Gamma)$ in $H_\Sigma(\Omega, \Gamma)$ given by Proposition 2.3. Its compactness results from Lemma 5.6 of [13] and Lemma 1.57 of [11]. Then, the coercivity, consistency, limit conformity and compactness of the sequence of HFV discretizations $(\mathcal{D}^l)_{l \in \mathbb{N}}$ results from Lemma 3.1 and Lemma 4.2. $\square$

**Remark 4.1** The proofs of Propositions 4.2, 4.4 and of Lemma 3.1 show that for solutions $u \in V^0$ and $(\mathbf{q}_m, \mathbf{q}_f) \in H_\Sigma(\Omega, \Gamma)$ of (4) such that $u \in C^2(\bar{K})$, $\mathbf{q}_m \in (C^1(\bar{K}))^d$, $\mathbf{q}_f \in (C^1(\bar{\sigma}))^{d-1}$ for all $K \in \mathcal{M}$ and all $\sigma \in \mathcal{F}_\Gamma$, the VAG and HFV schemes are consistent and limit conforming of order 1, and therefore convergent of order 1.

## 4.3 Finite Volume Formulation of the VAG and HFV schemes

Both the VAG and HFV schemes can be formulated as finite volume schemes. Moreover, the definition of the fluxes and of the conservation equations for both schemes can be unified using the following data structure which has been used in the practical implementation of the code. Let us define the set of degrees of freedom (d.o.f.)

$$ \mathrm{dof}_{\mathcal{D}} = \begin{cases} \mathcal{M} \cup \mathcal{V} \cup \mathcal{F}_{\Gamma} & \text{for VAG} \\ \mathcal{M} \cup \mathcal{F} \cup \mathcal{E}_{\Gamma} & \text{for HFV} \end{cases} $$

The subset of d.o.f. located at the boundary of $\Omega$ where Dirichlet boundary conditions are imposed is denoted by $\mathrm{dof}_{\mathrm{Dir}}$ such that

$$ \mathrm{dof}_{\mathrm{Dir}} = \begin{cases} \mathcal{V}_{\mathrm{ext}} & \text{for VAG} \\ (\mathcal{F} \setminus \mathcal{F}_{\mathrm{int}}) \cup (\mathcal{E}_{\Gamma} \setminus \mathcal{E}_{\mathrm{int}}) & \text{for HFV} \end{cases} $$

For each cell $K \in \mathcal{M}$ let us also define the subset $\mathrm{dof}_{\partial K}$ of d.o.f. located at the boundary of K with

$$ \mathrm{dof}_{\partial K} = \begin{cases} \mathcal{V}_K \cup (\mathcal{F}_K \cap \mathcal{F}_\Gamma) & \text{for VAG} \\ \mathcal{F}_K & \text{for HFV} \end{cases} $$
Similarly, we define for each fracture face $\sigma \in \mathcal{F}_\Gamma$ the subset $dof_{\partial\sigma}$ of d.o.f. located at the boundary of $\sigma$

$$dof_{\partial\sigma} = \begin{cases} \mathcal{V}_{\sigma} & \text{for VAG} \\ \mathcal{E}_{\sigma} & \text{for HFV} \end{cases}$$

Finally, the vector space $X_D$ is identified to $\mathbb{R}^{dof_D}$ and its subspace $X_D^0$ to $\mathbb{R}^{dof_D \setminus dof_{Dir}}$ and we denote by $(e_\nu, \nu \in dof_D)$ the canonical basis of $X_D$.

Using these unified notations and following [6] Section 4 for the VAG scheme and [13] Remark 2.7 for the HFV scheme, we can define for both schemes the matrix and fracture fluxes which are exhibited in Figure 4 inside a cell $K$ with a fracture face $\sigma$. The matrix fluxes connect a cell $K \in \mathcal{M}$ to its boundary d.o.f. $\nu \in dof_{\partial K}$ and are defined for all $u_D \in X_D$ by

$$F_{K,\nu}(u_D) = - \int_K \Lambda_m \nabla_{D_m} u_D \cdot \nabla_{D_m} e_\nu d\mathbf{x} = \sum_{\nu' \in dof_{\partial K}} T_{K,\nu'}^{\nu'}(u_K - u_{\nu'}),$$

with

$$T_{K,\nu}^{\nu'} = \int_K \Lambda_m \nabla_{D_m} e_{\nu'} \cdot \nabla_{D_m} e_\nu d\mathbf{x}.$$

Similarly, the fracture fluxes connect a fracture face $\sigma \in \mathcal{F}_\Gamma$ to its boundary d.o.f. $\nu \in dof_{\partial\sigma}$ and are defined for all $u_D \in X_D$ by

$$F_{\sigma,\nu}(u_D) = - \int_{\sigma} d_f \Lambda_f \nabla_{D_f} u_D \cdot \nabla_{D_f} e_\nu d\tau(\mathbf{x}) = \sum_{\nu' \in dof_{\partial\sigma}} T_{\sigma,\nu'}^{\nu'}(u_K - u_{\nu'}),$$

with

$$T_{\sigma,\nu}^{\nu'} = \int_{\sigma} d_f \Lambda_f \nabla_{D_f} e_{\nu'} \cdot \nabla_{D_f} e_\nu d\tau(\mathbf{x}).$$

These fluxes are such that for all $(u_D, v_D) \in X_D \times X_D$ one has

$$\int_{\Omega} \Lambda_m(\mathbf{x}) \nabla_{D_m} u_D(\mathbf{x}) \cdot \nabla_{D_m} v_D(\mathbf{x}) d\mathbf{x} + \int_{\Gamma} \Lambda_f(\mathbf{x}) \nabla_{D_f} u_D(\mathbf{x}) \cdot \nabla_{D_f} v_D(\mathbf{x}) d\tau_f(\mathbf{x}) \\ = \sum_{K \in \mathcal{M}} \sum_{\nu \in dof_{\partial K}} F_{K,\nu}(u_D)(v_K - v_\nu) + \sum_{\sigma \in \mathcal{F}_\Gamma} \sum_{\nu \in dof_{\partial\sigma}} F_{\sigma,\nu}(u_D)(v_\sigma - v_\nu).$$

It follows that the variational formulation (11) is equivalent to the following finite volume formulation: Find $u_D \in X_D$ such that

$$\left\{
\begin{array}{l}
\displaystyle \sum_{\nu \in dof_{\partial K}} F_{K,\nu}(u_D) = \int_{\omega_K} h_m(\mathbf{x}) d\mathbf{x}, \quad \text{for all } K \in \mathcal{M}, \\[2ex]
\displaystyle \sum_{\nu \in dof_{\partial \sigma}} F_{\sigma,\nu}(u_D) + \sum_{K \in \mathcal{M}_{\sigma}} -F_{K,\sigma}(u_D) \\
\displaystyle = \sum_{K \in \mathcal{M}_{\sigma}} \int_{\omega_{K,\sigma}} h_m(\mathbf{x}) d\mathbf{x} + \int_{\omega_{\sigma}} h_f(\mathbf{x}) d\tau_f(\mathbf{x}), \quad \text{for all } \sigma \in \mathcal{F}_{\Gamma}, \\[2ex]
\displaystyle \sum_{K \in \mathcal{M}_{\nu}} -F_{K,\nu}(u_D) + \sum_{\sigma \in \mathcal{F}_{\Gamma,\nu}} -F_{\sigma,\nu}(u_D) = \sum_{K \in \mathcal{M}_{\nu}} \int_{\omega_{K,\nu}} h_m(\mathbf{x}) d\mathbf{x} \\
\displaystyle \quad + \sum_{\sigma \in \mathcal{F}_{\Gamma,\nu}} \int_{\omega_{\sigma,\nu}} h_f(\mathbf{x}) d\tau_f(\mathbf{x}), \quad \text{for all } \nu \in dof_D \setminus (\mathcal{M} \cup \mathcal{F}_{\Gamma} \cup dof_{Dir}), \\[2ex]
u_{\nu} = 0, \quad \text{for all } \nu \in dof_{Dir},
\end{array}
\right.
\tag{35}$$
<figure><img src="image_6.png" /><figcaption>Figure 4: Matrix (in blue) and fracture (in red) fluxes inside a cell K with a fracture face σ (in bold) for the VAG (left) and HFV (right) schemes. The matrix fluxes $F_{K,ν}$ connect the cell K to the d.o.f. $ν ∈ dof_{∂K}$ located at the boundary of K. The fracture fluxes $F_{σ,ν}$ connect the face σ to the d.o.f. $ν ∈ dof_{∂σ}$ located at the boundary of σ.</figcaption></figure>

with $M_ν = \{K \in M | ν \in dof_{∂K}\}$ and $F_{Γ,ν} = \{σ \in F_Γ | ν \in dof_{∂σ}\}$. Following [6], when applying the VAG or HFV discretization to two phase Darcy flow models or to the coupling of the Darcy flow equation with a tracer equation, the choice of the cells and fracture faces partitioning defining the control volumes is done in order to avoid the mixture of heterogeneous properties inside each control volume. In particular, at the matrix fracture interfaces, one simply need to set $ω_{K,ν} = \emptyset$ for all $ν$ such that $x_ν \in Γ$. Note also that, in practice for such models, one does not need to build the partitions but only to choose the volume distribution ratios $α_{K,ν} = \frac{\int_{ω_{K,ν}} dx}{|K|}$, $ν \in dof_{∂K} \setminus dof_{Dir}$, and $α_{σ,ν} = \frac{\int_{ω_{σ,ν}} dτ_f(x)}{|σ|}$, $ν \in dof_{∂σ} \setminus dof_{Dir}$.

# 5 Density results for pressure and flux spaces

This section proves the density results stated in Propositions 2.2 and 2.3.

## 5.1 Proof of Proposition 2.2

We prove in this subsection that the space $C^∞(\bar{\Omega})$ is a dense subspace of $V$. The density of $C_c^∞(\bar{\Omega})$ in $V^0$ can be proved in the same way. We begin with technical lemmas concerning the fractional Sobolev spaces.

Let $s$ be a non-negative real number, for an non-empty open set $\Omega \subseteq \mathbb{R}^d$ we denote by $H^s(\Omega)$ the standard Sobolev space $W^{s,p}(\Omega)$ with exponent $p = 2$. Note that for $0 \le s_1 \le s_2$ one has $H^{s_2}(\Omega) \subset H^{s_1}(\Omega)$, where $H^0(\Omega)$ stands for $L^2(\Omega)$. We also define the spaces $H_0^s(\Omega) = \overline{C_c^∞(\Omega)}^{H^s(\Omega)}$ and $\tilde{H}^s(\Omega) = \overline{C_c^∞(\Omega)}^{H^s(\mathbb{R}^d)}$.

If $\Omega$ is a Lipschitz domain, then the trace operator $γ_{∂\Omega} : C^∞(\bar{\Omega}) \to C^∞(\partial\Omega)$ can be continuously extended to $γ_{∂\Omega} : H^s(\Omega) \to H^{s-1/2}(\partial\Omega)$ for all $\frac{1}{2} < s < \frac{3}{2}$; moreover the trace operator is surjective and has a continuous right inverse [9, Lemma 3.6]. For $s > \frac{1}{2}$ we define

$$H_{\partial\Omega}^s(\Omega) = \{u \in H^s(\Omega) \text{ such that } \gamma_{\partial\Omega}u = 0\}.$$

The following lemma summarizes some relations between the spaces $H^s(\Omega)$, $H_0^s(\Omega)$, $\tilde{H}^s(\Omega)$ and $H_{\partial\Omega}^s(\Omega)$ for Lipschitz domains.
**Lemma 5.1** If $\Omega$ is a Lipschitz domain, then

(i) If $s \ge 0$, then $\tilde{H}^s(\Omega) \subset H_0^s(\Omega)$; moreover if $s - \frac{1}{2} \notin \mathbb{Z}$, then $\tilde{H}^s(\Omega) = H_0^s(\Omega)$.

(ii) If $0 \le s \le \frac{1}{2}$, then $H_0^s(\Omega) = H^s(\Omega)$.

(iii) If $\frac{1}{2} < s < \frac{3}{2}$, then $H_0^s(\Omega) = H_{\partial\Omega}^s(\Omega)$.

**Proof:** See [18, Theorems 1.4.4.5 and 1.4.2.4] for the proof of (i) and (ii). The proof of (iii) for $\frac{1}{2} < s \le 1$ is also given in [22, Theorem 3.40], for $1 < s < \frac{3}{2}$ see [23]. $\square$

Lemma 5.1 implies in particular that for $\frac{1}{2} < s < \frac{3}{2}$ one has $H_0^s(\Omega) = \tilde{H}^s(\Omega) = H_{\partial\Omega}^s(\Omega)$ and that for any $u \in H_0^s(\Omega)$ its extension by zero belongs to $H^s(\mathbb{R}^n)$.

Let $\Omega$ be a Lipschitz domain, let $l \in \mathbb{N}$ and let $(\Omega_i \subset \Omega)_{i \in \{1, \dots, l\}}$ be a family of Lipschitz domains satisfying $\bigcup_{i \in \{1, \dots, l\}} \bar{\Omega}_i = \bar{\Omega}$ and $\Omega_i \cap \Omega_j = \emptyset$ for all $i, j \in \{1, \dots, l\}$ with $i \ne j$. Let $(u_i \in H^1(\Omega_i))_{i \in \{1, \dots, l\}}$ be a set of functions such that for all $i, j \in \{1, \dots, l\}$ satisfying $\partial\Omega_i \cap \partial\Omega_j \ne \emptyset$, then $(\gamma_{\partial\Omega_i} u_i)|_{\partial\Omega_i \cap \partial\Omega_j} = (\gamma_{\partial\Omega_j} u_i)|_{\partial\Omega_i \cap \partial\Omega_j}$. For all $i \in \{1, \dots, l\}$ we denote by $\chi_{\Omega_i}$ the indicator function of the domain $\Omega_i$, it is classical that the function $u = \sum_{i \in \{1, \dots, l\}} u_i \chi_{\Omega_i}$ belongs to $H^1(\Omega)$.

The following lemma gives a similar result for functions in $H^s(\Omega)$, with $\frac{1}{2} < s < 1$ and $\Omega$ being the boundary of a polyhedral domain.

**Lemma 5.2** Let $K$ be a bounded polyhedral domain in $\mathbb{R}^3$, let $\mathcal{F}_K$ be the set of its polygonal faces and $\mathcal{E}_K$ the set of its edges. For all $\sigma \in \mathcal{F}_K$ we denote by $\mathcal{E}_\sigma$ the set of edges of $\sigma$ and for all $e \in \mathcal{E}_K$ by $\mathcal{F}_{K,e} \subset \mathcal{F}_K$ the two faces containing $e$. Let $\frac{1}{2} < s < 1$, for all $\sigma \in \mathcal{F}_K$ and $e \in \mathcal{E}_\sigma$ we denote by $\gamma_{\sigma,e}$ the trace operator from $H^s(\sigma)$ to $H^{s-\frac{1}{2}}(e)$.

Let $u \in L^2(\partial K)$ such that $u|_{\sigma} \in H^s(\sigma)$ for all $\sigma \in \mathcal{F}_K$ and such that for all $e \in \mathcal{E}_K$ with $\mathcal{F}_{K,e} = \{\sigma, \sigma'\}$ one has $\gamma_{\sigma,e}(u|_{\sigma}) = \gamma_{\sigma',e}(u|_{\sigma'})$. Then, $u \in H^s(\partial K)$.

**Proof:** Let $e \in \mathcal{E}_K$ and $\mathcal{F}_{K,e} = \{\sigma, \sigma'\}$, we associate with $e$ a couple of Lipschitz domains $(D_e^\sigma)_{\sigma \in \mathcal{F}_{K,e}}$ such that $D_e^\sigma \subset \sigma$ and $\partial D_e^\sigma \cap \partial \sigma = e$ for all $\sigma \in \mathcal{F}_{K,e}$, and such that $D_e^{\sigma'}$ is obtained by a rotation of $D_e^\sigma$ around $e$. We denote $\bar{D}_e = \bar{D}_e^\sigma \cup \bar{D}_e^{\sigma'}$ and $D_e = \bar{D}_e \setminus \partial \bar{D}_e$. The function $\gamma_{\sigma,e} u_e^\sigma$ is in $H^{s-\frac{1}{2}}(e)$, implying in view of statements (i) and (ii) of Lemma 5.1 that its extension by zero on $\partial D_e^\sigma$ belongs to $H^{s-\frac{1}{2}}(\partial D_e^\sigma)$. The trace operator $\gamma_{\partial D_e^\sigma}: H^s(D_e^\sigma) \to H^{s-\frac{1}{2}}(\partial D_e^\sigma)$ is surjective, therefore there exists $u_e^\sigma \in H^s(D_e^\sigma)$ such that

$$ \gamma_{\partial D_e^\sigma} u_e^\sigma = \begin{cases} \gamma_{\sigma,e} u_\sigma & \text{on } e, \\ 0 & \text{on } \partial D_e^\sigma \setminus e. \end{cases} $$

We denote by $u_e$ the extension by symmetry of $u_e^\sigma$ to $D_e$. One can show that $u_e \in H_0^s(D_e)$, implying that its extension by zero to the whole $\partial K$, denoted by $\bar{u}_e$, belongs to $H^s(\partial K)$. Let us consider the function

$$ v_\sigma = u|_{\sigma} - \sum_{e \in \mathcal{E}_\sigma} \bar{u}_e|_{\sigma}. $$
Since $v_\sigma \in H_0^s(\sigma)$ it can be extended by zero to the whole $\partial K$; we denote this extension by $\bar{v}_\sigma$. To complete the proof we remark that

$$u = \sum_{\sigma \in \mathcal{F}_K} \bar{v}_\sigma + \sum_{e \in \mathcal{E}_K} \bar{u}_e \in H^s(\partial K).$$

<figure><img src="image_7.png" /></figure>

**Lemma 5.3** The trace operator $\gamma$ is onto from $V$ to $H^1(\Gamma)$.

**Proof:** Given $u \in H^1(\Gamma)$ we prove that there exists $U \in V$ such that $\gamma U = u$. We focus on the case $d = 3$, the adaptation to the bidimensional setting is straightforward. The proof relies on the definition of a polyhedral mesh which is slightly different from Definition 4.1. More precisely in addition to Definition 4.1 it is assumed that all the mesh faces $\sigma \in \mathcal{F}$ are planar. On the other hand the existence of cell (and face) “centers” is not required. Such polyhedral partitioning of $\Omega \setminus \bar{\Gamma}$ always exists. For example we can define the set of cells $\mathcal{M}$ as the set of connected components of $\Omega \setminus \bigcup_{i \in I} \mathcal{P}_i$. In addition to the previous notations, we will denote by $\mathcal{F}_e$ the set of faces sharing a given edge $e \in \mathcal{E}$.

For all $e \in \mathcal{E}_\Gamma$ we denote by $\gamma_e$ the trace operator acting from $H^1(\Gamma)$ to $H^{\frac{1}{2}}(e)$, and for all $e \in \mathcal{E}$ we define

$$u_e = \begin{cases} \gamma_e u & \text{if } e \in \mathcal{E}_\Gamma, \\ 0 & \text{else.} \end{cases}$$

Let $\sigma \in \mathcal{F} \setminus \mathcal{F}_\Gamma$ and let us denote by $\chi_e$, $e \in \mathcal{E}_\sigma$ the indicator function of $e$ defined on $\partial \sigma$. In view of statements (i) and (ii) of Lemma 5.1 the function $u_{\partial \sigma} = \sum_{e \in \mathcal{E}_\sigma} \chi_e u_e$ belongs to $H^s(\partial \sigma)$ for any $0 < s < \frac{1}{2}$.

For all $\sigma \in \mathcal{F} \setminus \mathcal{F}_\Gamma$ there exists a function $u_\sigma \in H^{s+\frac{1}{2}}(\sigma)$ having the trace on $\partial \sigma$ equal to $u_{\partial \sigma}$. For $\sigma \in \mathcal{F}_\Gamma$ we denote by $u_\sigma \in H^1(\sigma)$ the restriction of $u$ on $\sigma$. Let $\chi_\sigma$ be the indicator function of $\sigma$ defined on $\bigcup_{\sigma \in \mathcal{F}} \sigma$, we set $u^* = \sum_{\sigma \in \mathcal{F}} \chi_\sigma u_\sigma$. Lemma 5.2 implies that for all $K \in \mathcal{M}$ the restriction of $u^*$ on $\partial K$ belongs $H^{s+\frac{1}{2}}(\partial K)$

Finally, for all $K \in \mathcal{M}$ we denote by $U_K$ some lifting of $u^*|_{\partial K}$ in $H^{s+1}(K)$ and we define $U = \sum_{K \in \mathcal{M}} \chi_K U_K$, where $\chi_K$ is the indicator function of a cell $K$. Then, it follows that $U \in H^1(\Omega)$ and that $\gamma U = u$. $\square$

Let us recall (see e.g. [7] Proposition 1.9) that the density of a linear subspace $M$ of $V$ is equivalent to the fact that any linear form of $V'$ vanishing on $M$ is vanishing on $V$.

The characterization of $V'$ can be obtained from the Riesz theorem, implying that any continuous linear form $\zeta$ on $V$ writes $\zeta = \xi + \gamma^t \tau$ where $\xi \in (H^1(\Omega))'$ and $\tau \in (H^1(\Gamma))'$. Then, assuming that $\langle \zeta, \varphi \rangle = 0$ for all $\varphi \in C^\infty(\bar{\Omega})$ it results from Lemma 5.4 stated and proved below that $\zeta = 0$. Therefore the space $C^\infty(\bar{\Omega})$ is dense in $V$.

**Lemma 5.4** Let $\zeta = \xi + \gamma^t \tau$ where $\xi \in (H^1(\Omega))'$ and $\tau \in (H^1(\Gamma))'$ be such that $\langle \zeta, \varphi \rangle = 0$ for all $\varphi \in C^\infty(\bar{\Omega})$, then $\zeta = 0$.

**Proof:** It is known that $C_c^\infty(\bar{\Omega} \setminus \bar{\Gamma})$, defined as the space of $C^\infty(\bar{\Omega})$ functions vanishing in a neighbourhood of $\Gamma$, is a dense subspace of $H_\Gamma^1(\Omega \setminus \bar{\Gamma})$ defined as the space of $H^1(\Omega \setminus \bar{\Gamma})$ functions vanishing on $\Gamma$. From the surjectivity and continuity of the trace operator $\gamma$ from $V$ to $H^1(\Gamma)$, there exists a continuous lifting operator denoted by $r_\Gamma$ from $H^1(\Gamma)$ to $V$. Using $\langle \xi, \varphi \rangle = 0$ for all $\varphi \in C_c^\infty(\bar{\Omega} \setminus \bar{\Gamma})$, and the density of $C_c^\infty(\bar{\Omega} \setminus \bar{\Gamma})$ in $H_\Gamma^1(\Omega \setminus \bar{\Gamma})$, we deduce that
< ξ, v > = < ξ, r_Γ(γ(v)) > for all v ∈ V. It results that ˜τ = r_Γ^t ξ ∈ (H^1(Γ))' is such that ξ = γ^t ˜τ. Hence, we can assume in the remaining of the proof that ξ = 0.

Let us set $E = \bigcup_{i \in I} \partial \Gamma_i$. We have $E = \bigcup_{e \in \mathcal{E}_\Gamma} \bar{e}$ where $\mathcal{E}_\Gamma$ is the set of edges of $\Gamma$ in the mesh defined above. Let $\gamma_E$ the trace operator from $H^1(\Gamma)$ to $L^2(E)$ and let us define the space $H^{\frac{1}{2}}(E)$ as $\gamma_E(H^1(\Gamma))$. We also define the space $H^{\frac{1}{2}}(\partial \Gamma_i)$ as the set of traces on $\partial \Gamma_i$ of functions in $H^1(\Gamma_i)$. Then, it is easy to show that a function $v \in L^2(E)$ belongs to $H^{\frac{1}{2}}(E)$ iff for all $i \in I$, the restriction $v_i$ of $v$ to $\partial \Gamma_i$ belongs to $H^{\frac{1}{2}}(\partial \Gamma_i)$. The function space $H^{\frac{1}{2}}(E)$ is endowed with the Hilbertian norm

$$
\|v\|_{H^{\frac{1}{2}}(E)} = \sum_{i \in I} \left( \|v_i\|_{H^{\frac{1}{2}}(\partial \Gamma_i)}^2 \right)^{\frac{1}{2}}.
$$

From the continuity and surjectivity of the trace operator $\gamma_E$ from $H^1(\Gamma)$ to $H^{\frac{1}{2}}(E)$, we deduce that there exists a continuous lifting operator denoted by $r_E$ from $H^{\frac{1}{2}}(E)$ to $H^1(\Gamma)$. Let us denote by $H_E^1(\Gamma \setminus E)$ the subspace of functions in $H^1(\Gamma)$ with a vanishing trace on $E$. From the known density of $C_c^\infty(\Gamma \setminus E)$ in $H^1(\Gamma \setminus E)$, we deduce as above that $l = r_E^t \tau \in (H^{\frac{1}{2}}(E))'$ is such that $\tau = \gamma_E^t l$

Let us denote by $\mathcal{V}_E$ the set of the vertices of $E$. For all $\phi \in C_c^\infty(E \setminus \mathcal{V}_E)$, there exists $\varphi \in C^\infty(\bar{\Omega})$ such that $\phi = \varphi|_E$. Hence, $l \in (H^{\frac{1}{2}}(E))'$ is such that $\langle l, \phi \rangle = 0$ for all $\phi \in C_c^\infty(E \setminus \mathcal{V}_E)$. Since $C_c^\infty(e)$ is dense in $H^{\frac{1}{2}}(e)$ for any edge $e$ ([18] Theorem 1.4.2.4), we can deduce that $C_c^\infty(E \setminus \mathcal{V}_E)$ is dense in $H^{\frac{1}{2}}(E)$ and hence that $l = 0$, and then that $\zeta = 0$. $\square$.

5.2 Proof of Proposition 2.3

We prove in this subsection that the space $C_{H\Sigma}^{\infty}(\Omega, \Gamma)$ is a dense subspace of $H_{\Sigma}(\Omega, \Gamma)$.

Let us start by the following technical lemma.

**Lemma 5.5** Let $v_m \in L^2(\Omega)$, $g_m \in L^2(\Omega)^d$, and $v_f \in L^2(\Gamma)$, $g_f \in L^2(\Gamma)^{d-1}$ be such that

$$
a_{\Sigma}((\mathbf{q}_m, \mathbf{q}_f), (v_m, v_f, \mathbf{g}_m, \mathbf{g}_f)) = 0 \quad \text{for all } (\mathbf{q}_m, \mathbf{q}_f) \in C_{H_{\Sigma}}^{\infty}(\Omega, \Gamma), \quad (36)
$$

then $v_m \in V^0$, $v_f = \gamma v_m$, $g_m = \nabla v_m$ and $g_f = \nabla_\tau \gamma v_m$.

**Proof:** It follows from (36) that $\sum_{\alpha \in A} \int_{\Omega_\alpha} (v_m \mathrm{div}(\mathbf{q}_{m,\alpha}) + \mathbf{q}_{m,\alpha} \cdot \mathbf{g}_m) d\mathbf{x} = 0$, for all $\mathbf{q}_m \in C_c^\infty(\Omega \setminus \bar{\Gamma})^d$, we classically deduce that $v_m \in H^1(\Omega \setminus \bar{\Gamma})$ and that $\mathbf{g}_m = \nabla v_m$. Next, let us denote by $\gamma_i^\pm : H^1(\Omega \setminus \bar{\Gamma}) \to L^2(\Gamma_i)$ the trace operators on $\Gamma_i$ from the sides $\pm$ of $\Omega \setminus \bar{\Gamma}$. For any given $i \in I$, let us denote by $\omega_{\Gamma_i}$ any open ball of $\Gamma_i$. For any $\varphi_i \in C_c^\infty(\omega_{\Gamma_i})$, one can build a function $\mathbf{q}_{m,\alpha_i^+} \in C_b^\infty(\Omega_{\alpha_i^+})$ such that $\mathbf{q}_{m,\alpha_i^+} \cdot \mathbf{n}_i^+ = \varphi_i$ on $\omega_{\Gamma_i}$ and $\mathbf{q}_{m,\alpha_i^-} \cdot \mathbf{n}_i^- = 0$ on $\omega_{\Gamma_i}$ if $\alpha_i^+ = \alpha_i^-$, $\mathbf{q}_{m,\alpha_i^+} \cdot \mathbf{n} = 0$ on $\partial\Omega_{\alpha_i^+} \cap \partial\Omega$, and $\mathbf{q}_{m,\alpha_j^+} \cdot \mathbf{n}_j^+ = 0$ on the sides $\pm$ of the fractures $i \neq j \in I$ such that $\alpha_j^+ = \alpha_i^+$. It follows from (36) that

$$
\int_{\omega_{\Gamma_i}} (\gamma_i^+ v_m - v_f) \varphi_i d\tau(\mathbf{x}) = 0.
$$

which implies that $\gamma_i^+ v_m = v_f|_{\Gamma_i}$ in $L^2(\Gamma_i)$ for all $i \in I$. Similarly, we can show that $\gamma_i^- v_m = v_f|_{\Gamma_i}$ in $L^2(\Gamma_i)$ for all $i \in I$, and $v_m = 0$ on $\partial\Omega$. Hence we deduce that $v_m \in H_0^1(\Omega)$ with $v_f = \gamma v_m$. Next, for all $i \in I$ and for all $\mathbf{q}_{f,i} \in C_c^\infty(\Gamma_i)^{d-1}$, one has

$$
\int_{\Gamma_i} (v_f \mathrm{div}_{\tau_i}(\mathbf{q}_{f,i}) + \mathbf{g}_f \cdot \mathbf{q}_{f,i}) d\tau(\mathbf{x}) = 0,
$$
which implies that $v_f|_{\Gamma_i} \in H^1(\Gamma_i)$ with $g_f|_{\Gamma_i} = \nabla_{\tau_i} v_f|_{\Gamma_i}$. Next, for all $(i, j) \in I \times I$, $i \neq j$, such that $\Sigma_{i,j} \setminus \Sigma_0$ is of codimension 2 non zero measure, let us consider any open segment $L_{i,j} \subset \Sigma_{i,j} \setminus \Sigma_0$, $r_{i,j} \in C_c^\infty(L_{i,j})$ and its extension $r_i$ (resp. $r_j$) by zero outside of $L_{i,j}$ on $\Sigma_i$ (resp. $\Sigma_j$). Then, one can build $(0, q_f) \in C_{H_\Sigma}^\infty(\Omega, \Gamma)$ such that $q_{f,i} \cdot n_{\Sigma_i} = r_i$, and $q_{f,j} \cdot n_{\Sigma_j} = -r_j$ and $q_{f,l} = 0$ for all $l \in I \setminus \{i, j\}$. It follows from (36) that

$$ \int_{L_{i,j}} \left( (\mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i}) \gamma_{\Sigma_i} v_f|_{\Gamma_i} + (\mathbf{q}_{f,j} \cdot \mathbf{n}_{\Sigma_j}) \gamma_{\Sigma_j} v_f|_{\Gamma_j} \right) dl(\mathbf{x}) \\ = \int_{L_{i,j}} r_{i,j} (\gamma_{\Sigma_i} v_f|_{\Gamma_i} - \gamma_{\Sigma_j} v_f|_{\Gamma_j}) dl(\mathbf{x}) = 0, $$

and hence that $\gamma_{\Sigma_i} v_f|_{\Gamma_i} = \gamma_{\Sigma_j} v_f|_{\Gamma_j}$ on $L_{i,j}$. Also, for all $i \in I$ such that $\Sigma_{i,0}$ has a non vanishing codimension 2 measure, let us consider any open segment $L_{i,0} \subset \Sigma_{i,0}$, $r_{i,0} \in C_c^\infty(L_{i,0})$ and its extension $r_i$ by zero outside of $L_{i,0}$ on $\Sigma_i$. Then, one can build $(0, q_f) \in C_{H_\Sigma}^\infty(\Omega, \Gamma)$ such that $q_{f,i} \cdot n_{\Sigma_i} = r_i$ and $q_{f,l} = 0$ for all $l \in I \setminus \{i\}$. It results from (36) that

$$ \int_{L_{i,0}} (\mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i}) \gamma_{\Sigma_i} v_f|_{\Gamma_i} dl(\mathbf{x}) = \int_{L_{i,0}} (\gamma_{\Sigma_i} v_f|_{\Gamma_i}) r_{i,0} dl(\mathbf{x}) = 0, $$

and hence that $\gamma_{\Sigma_i} v_f|_{\Gamma_i} = 0$ on $L_{i,0}$ which implies together with the previous properties that $v_f \in H_{\Sigma_0}^1(\Gamma)$ and concludes the proof. $\square$

Let us now prove the density of $C_{H_\Sigma}^\infty(\Omega, \Gamma)$ in $H_\Sigma(\Omega, \Gamma)$. Similarly as in the previous subsection, it is equivalent to prove that any linear form in $H_\Sigma(\Omega, \Gamma)'$ which vanishes on the subspace $C_{H_\Sigma}^\infty(\Omega, \Gamma)$, vanishes on the whole space $H_\Sigma(\Omega, \Gamma)$.

Since $H_\Sigma(\Omega, \Gamma)$ is a closed Hilbert subspace of $H(\Omega, \Gamma)$, any continuous linear form on $H_\Sigma(\Omega, \Gamma)$ can be continuously extended to $H(\Omega, \Gamma)$. From the Riez representation theorem, a continuous linear form $\xi$ on $H(\Omega, \Gamma)$ writes for all $(\mathbf{q}_m, \mathbf{q}_f) \in H(\Omega, \Gamma)$

$$ \langle \xi, (\mathbf{q}_m, \mathbf{q}_f) \rangle = a_\Sigma \left( (\mathbf{q}_m, \mathbf{q}_f), (v_m, v_f, \mathbf{g}_m, \mathbf{g}_f) \right) $$

with $v_m \in L^2(\Omega)$, $g_m \in L^2(\Omega)^d$, and $v_f \in L^2(\Gamma)$, $g_f \in L^2(\Gamma)^{d-1}$. It is assumed that $\xi$ vanishes on $C_{H_\Sigma}^\infty(\Omega, \Gamma)$. Then, it results from Lemma 5.5 that $v_m \in V^0$, $v_f = \gamma v_m$, $g_m = \nabla v_m$ and $g_f = \nabla_\tau \gamma v_m$. From the definition (2) of $H_\Sigma(\Omega, \Gamma)$, we conclude that $\xi$ vanishes on $H_\Sigma(\Omega, \Gamma)$ which proves the density of $C_{H_\Sigma}^\infty(\Omega, \Gamma)$ in $H_\Sigma(\Omega, \Gamma)$.

# 6 Numerical experiments

Let $\Omega$ denote a bounded domain in $\mathbb{R}^d$, $d = 3$, and let us consider four non immersed planar fractures splitting the domain $\Omega$ into four subdomains $\Omega_\alpha$, $\alpha = 1, \dots, 4$. Dirichlet boundary conditions are imposed on both the boundary of the domain $\partial\Omega$ and on the boundary of the fracture network $\partial\Gamma = \partial\Omega \cap \Gamma = \Sigma_0$. The permeability of the fractures is defined by $\Lambda_f(\mathbf{x}) = 100 \text{ Id}$ and their width by $d_f(\mathbf{x}) = 0.01$. In the matrix, the permeability tensor $\Lambda_m(\mathbf{x})$ is fixed to $\Lambda_{m,\alpha}$ on each subdomain $\Omega_\alpha$, $\alpha = 1, \dots, 4$ with two choices of the subdomain permeabilities. The first choice considers isotropic heterogeneous permeabilities setting $\Lambda_{m,\alpha} = \lambda_\alpha \text{ Id}$ with $\lambda_1 = 1$, $\lambda_2 = 0.1$, $\lambda_3 = 0.01$, $\lambda_4 = 10$. The second choice corresponds to the anisotropic heterogeneous permeabilities defined by

$$ \Lambda_{m,1} = \begin{pmatrix} a_1 & b_1 & 0 \\ b_1 & c_1 & 0 \\ 0 & 0 & \lambda \end{pmatrix}, \Lambda_{m,2} = \begin{pmatrix} a_2 & 0 & b_2 \\ 0 & \lambda & 0 \\ b_2 & 0 & c_2 \end{pmatrix}, \Lambda_{m,3} = \begin{pmatrix} a_3 & b_3 & 0 \\ b_3 & c_3 & 0 \\ 0 & 0 & \lambda \end{pmatrix}, \Lambda_{m,4} = \begin{pmatrix} \lambda & 0 & 0 \\ 0 & a_4 & b_4 \\ 0 & b_4 & c_4 \end{pmatrix}, $$
with $a_\alpha = \cos^2 \beta_\alpha + \omega \sin^2 \beta_\alpha$, $b_\alpha = (1-\omega)\cos \beta_\alpha \sin \beta_\alpha$, $c_\alpha = \omega \cos^2 \beta_\alpha + \sin^2 \beta_\alpha$, $\lambda = 0.01$, $\beta_1 = \frac{\pi}{6}$, $\beta_2 = -\frac{\pi}{6}$, $\beta_3 = 0$, $\beta_4 = \frac{\pi}{6}$ and $\omega = 0.01$.

Next, let us describe the two families of test cases that will be presented in this section.

*Test cases:* For the first test case $\Omega = (0,1)^3$, the fracture network is defined by the union of the two rectangles $\{(x,y,z) \in \Omega \mid x=0.5\}$ and $\{(x,y,z) \in \Omega \mid y=0.5\}$, and the four subdomains correspond to $\Omega_1 = \{(x,y,z) \in \Omega \mid x < 0.5, y < 0.5\}$, $\Omega_2 = \{(x,y,z) \in \Omega \mid x > 0.5, y < 0.5\}$, $\Omega_3 = \{(x,y,z) \in \Omega \mid x > 0.5, y > 0.5\}$ and $\Omega_4 = \{(x,y,z) \in \Omega \mid x < 0.5, y > 0.5\}$ (see the left picture of Figure 5). Let us define the functions $t_1(\mathbf{x}) = y-x+z$, $t_2(\mathbf{x}) = x+y+z-1$, $t_3(\mathbf{x}) = x-y+z$ and $t_4(\mathbf{x}) = 1-x-y+z$. One can check that the function $u(\mathbf{x}) = e^{\cos(t_\alpha(\mathbf{x}))}$, $\mathbf{x} \in \Omega_\alpha$, $\alpha = 1, \dots, 4$, belongs to $V$, and that the fluxes $(\mathbf{q}_m, \mathbf{q}_f) = (-\Lambda_m \nabla u, -d_f \Lambda_f \nabla_\tau \gamma u)$ belongs to $H_\Sigma(\Omega, \Gamma)$ since it satisfies $\sum_{i \in I} \mathbf{q}_{f,i} \cdot \mathbf{n}_{\Sigma_i} = 0$ on $\Sigma \setminus \Sigma_0$.

<figure><img src="image_8.png" /><figcaption>Figure 5: (Left): domain $\Omega = (0,1)^3$ and fracture network for the first test case. (Right): domain $\Omega = (-1.5, 1.5) \times (-2, 2) \times (0, 5)$ and fracture network for the second test case.</figcaption></figure>

For the numerical solutions using both the VAG and HFV schemes, three different families of meshes are considered. The first family is defined by uniform Cartesian grid of size $n \times n \times n$ with $n = 4, 8, 16, 32, 64, 128$ (see Table 1). The second family of meshes is obtained from the previous one by a perturbation of its nodes excluding the nodes on the boundary of $\Omega$ and on the boundary of each fracture $\Gamma_i$, $i \in I$. The perturbation is chosen randomly inside the ball of radius $\frac{1}{4n}$ and of center the Cartesian mesh node. The perturbation of a fracture node is done in the fracture plane. Note that it leads to hexahedral cells with non planar faces and hence the HFV scheme is no longer consistent on this family of meshes. Finally we consider a family of uniformly refined tetrahedral meshes generated by TetGen[29] (see Table 2). Tables 1 and 2 provide respectively for the Cartesian or randomly perturbed Cartesian meshes, and for the tetrahedral meshes, as well as for both schemes, the number of degrees of freedom (d.o.f.), the number of d.o.f. after elimination of the cell and Dirichlet unknowns (Reduced d.o.f.), and the number of nonzero element in the linear system after elimination without any fill-in of the cell and Dirichlet unknowns (Nonzero elts).
<table><caption>Table 1: For the first test case, the VAG and HFV schemes and the six Cartesian and randomly perturbed Cartesian meshes: mesh number, number of d.o.f., number of d.o.f. after elimination of the cell and Dirichlet unknowns, number of nonzero elements in the matrix after elimination.</caption><caption>Vertex Approximate Gradient Discretization</caption><thead><tr><th>Nb</th><th>d.o.f.</th><th>Reduced d.o.f.</th><th>Nonzero elts</th></tr></thead><tbody><tr><td>1</td><td>221</td><td>59</td><td>839</td></tr><tr><td>2</td><td>1 369</td><td>471</td><td>9 403</td></tr><tr><td>3</td><td>9 521</td><td>3 887</td><td>90 947</td></tr><tr><td>4</td><td>70 753</td><td>31 839</td><td>802 003</td></tr><tr><td>5</td><td>544 961</td><td>258 239</td><td>6 738 419</td></tr><tr><td>6</td><td>4 276 609</td><td>2 081 151</td><td>55 247 923</td></tr></tbody></table><table><caption>Hybrid Finite Volume Discretization</caption><thead><tr><th>Nb</th><th>d.o.f.</th><th>Reduced d.o.f.</th><th>Nonzero elts</th></tr></thead><tbody><tr><td>1</td><td>380</td><td>188</td><td>1 644</td></tr><tr><td>2</td><td>2 520</td><td>1 560</td><td>15 336</td></tr><tr><td>3</td><td>18 224</td><td>12 464</td><td>129 600</td></tr><tr><td>4</td><td>138 336</td><td>99 168</td><td>1 060 464</td></tr><tr><td>5</td><td>1 077 440</td><td>790 208</td><td>8 570 064</td></tr><tr><td>6</td><td>8 503 680</td><td>6 307 200</td><td>68 888 976</td></tr></tbody></table>

<table><caption>Table 2: For the first test case, the VAG and HFV schemes and the ten tetrahedral meshes: mesh number, number of d.o.f., number of d.o.f. after elimination of the cell and Dirichlet unknowns, number of nonzero elements in the matrix after elimination.</caption><caption>Vertex Approximate Gradient Discretization</caption><thead><tr><th>Nb</th><th>d.o.f.</th><th>Reduced d.o.f.</th><th>Nonzero elts</th></tr></thead><tbody><tr><td>1</td><td>1 888</td><td>294</td><td>2 660</td></tr><tr><td>2</td><td>13 593</td><td>1 924</td><td>23 148</td></tr><tr><td>3</td><td>121 818</td><td>16 780</td><td>233 978</td></tr><tr><td>4</td><td>263 391</td><td>36 214</td><td>519 694</td></tr><tr><td>5</td><td>509 038</td><td>69 762</td><td>1 021 940</td></tr><tr><td>6</td><td>939 007</td><td>128 324</td><td>1 904 390</td></tr><tr><td>7</td><td>1 386 833</td><td>189 300</td><td>2 830 880</td></tr><tr><td>8</td><td>1 874 186</td><td>255 370</td><td>3 840 778</td></tr><tr><td>9</td><td>2 383 038</td><td>324 682</td><td>4 901 360</td></tr><tr><td>10</td><td>4 813 285</td><td>654 670</td><td>9 979 004</td></tr></tbody></table><table><caption>Hybrid Finite Volume Discretization</caption><thead><tr><th>Nb</th><th>d.o.f.</th><th>Reduced d.o.f.</th><th>Nonzero elts</th></tr></thead><tbody><tr><td>1</td><td>4 569</td><td>2 661</td><td>17 677</td></tr><tr><td>2</td><td>34 150</td><td>21 409</td><td>146 147</td></tr><tr><td>3</td><td>311 261</td><td>201 748</td><td>1 395 908</td></tr><tr><td>4</td><td>675 298</td><td>440 798</td><td>3 058 868</td></tr><tr><td>5</td><td>1308 518</td><td>858 252</td><td>5 967 626</td></tr><tr><td>6</td><td>2 417 392</td><td>1 589 624</td><td>11 064 478</td></tr><tr><td>7</td><td>3 573 654</td><td>2 354 004</td><td>16 396 536</td></tr><tr><td>8</td><td>4 832 987</td><td>3 187 229</td><td>22 210 505</td></tr><tr><td>9</td><td>6 147 875</td><td>4 058 104</td><td>28 290 370</td></tr><tr><td>10</td><td>12 432 788</td><td>8 223 946</td><td>57 382 094</td></tr></tbody></table>

The second test case considers the domain $\Omega = (-1.5, 1.5) \times (-2, 2) \times (0, 5)$ and the fracture network defined by the union of two rectangles $\{(x, y, z) \in \Omega | y = mx\}$ and $\{(x, y, z) \in \Omega | y = -mx\}$ with $m = 8$ (see the right picture of Figure 5). The domain $\Omega$ is split into the following four subdomains: $\Omega_1 = \{(x, y, z) \in \Omega | mx < y, mx < -y\}$, $\Omega_2 = \{(x, y, z) \in \Omega | mx > y, mx < -y\}$, $\Omega_3 = \{(x, y, z) \in \Omega | mx > y, mx > -y\}$ and $\Omega_4 = \{(x, y, z) \in \Omega | mx < y, mx > -y\}$. In this test we set $t_1(\mathbf{x}) = 2y + z$, $t_2(\mathbf{x}) = 2mx + z$, $t_3(\mathbf{x}) = -2y + z$ and $t_4(\mathbf{x}) = -2mx + z$. It can be verified that the function $u(\mathbf{x}) = e^{\cos(t_\alpha(\mathbf{x}))}$, $\mathbf{x} \in \Omega_\alpha$, $\alpha = 1, \dots, 4$, belongs to $V$, and that the fluxes $(q_m, q_f) = (-\Lambda_m \nabla u, -d_f \Lambda_f \nabla \gamma u)$ belong to $H_\Sigma(\Omega, \Gamma)$. For the numerical solutions using the VAG and HFV schemes, a family of ten tetrahedral uniformly refined meshes is generated by TetGen [29]. Table 3 gives for both schemes, the number of degrees of freedom (d.o.f.), the number of d.o.f. after elimination of the cell and Dirichlet unknowns (Reduced d.o.f.), and the number of nonzero element in the linear system after elimination without any fill-in of the cell and Dirichlet unknowns (Nonzero elts).

<table><caption>Table 3: For the second test case, the VAG and HFV schemes and the ten tetrahedral meshes: mesh number, number of d.o.f., number of d.o.f. after elimination of the cell and Dirichlet unknowns, number of nonzero elements in the matrix after elimination.</caption><caption>Vertex Approximate Gradient Discretization</caption><thead><tr><th>Nb</th><th>d.o.f.</th><th>Reduced d.o.f.</th><th>Nonzero elts</th></tr></thead><tbody><tr><td>1</td><td>12 930</td><td>2 081</td><td>24 687</td></tr><tr><td>2</td><td>62 177</td><td>9 280</td><td>123 096</td></tr><tr><td>3</td><td>132 712</td><td>19 321</td><td>265 709</td></tr><tr><td>4</td><td>251 969</td><td>36 103</td><td>510 459</td></tr><tr><td>5</td><td>463 906</td><td>65 850</td><td>949 882</td></tr><tr><td>6</td><td>1 002 529</td><td>140 712</td><td>2 070 638</td></tr><tr><td>7</td><td>1 366 875</td><td>190 979</td><td>2 832 163</td></tr><tr><td>8</td><td>1 934 904</td><td>269 381</td><td>4 022 379</td></tr><tr><td>9</td><td>2 342 305</td><td>325 513</td><td>4 877 093</td></tr><tr><td>10</td><td>4 542 801</td><td>627 526</td><td>9 501 798</td></tr></tbody></table><table><caption>Hybrid Finite Volume Discretization</caption><thead><tr><th>Nb</th><th>d.o.f.</th><th>Reduced d.o.f.</th><th>Nonzero elts</th></tr></thead><tbody><tr><td>1</td><td>32 218</td><td>20 369</td><td>140 055</td></tr><tr><td>2</td><td>157 522</td><td>101 864</td><td>705 466</td></tr><tr><td>3</td><td>337 883</td><td>219 755</td><td>1 524 503</td></tr><tr><td>4</td><td>644 056</td><td>421 122</td><td>2 926 982</td></tr><tr><td>5</td><td>1 188 904</td><td>780 117</td><td>5 429 401</td></tr><tr><td>6</td><td>2 576 269</td><td>1 696 321</td><td>11 820 151</td></tr><tr><td>7</td><td>3 516 282</td><td>2 318 255</td><td>16 161 367</td></tr><tr><td>8</td><td>4 982 226</td><td>3 289 061</td><td>22 940 167</td></tr><tr><td>9</td><td>6 034 003</td><td>3 985 462</td><td>27 803 112</td></tr><tr><td>10</td><td>11 719 544</td><td>7 754 712</td><td>54 132 280</td></tr></tbody></table>
*Numerical Results:* All test cases are performed using the $\Pi_{D_m}$ and $\Pi_{D_f}$ function reconstruction operators obtained by setting $\omega_K = K$ for all $K \in \mathcal{M}$, and $\omega_\sigma = \sigma$ for all $\sigma \in \mathcal{F}_\Gamma$. To assess the error estimates of Proposition 3.3, we have computed the following relative errors

$$ \mathrm{Err}_u = \frac{\|u - \Pi_{D_m} u_D\|_{L^2(\Omega)}}{\|u\|_{L^2(\Omega)}} + \frac{\|\gamma u - \Pi_{D_f} u_D\|_{L^2(\Gamma)}}{\|\gamma u\|_{L^2(\Gamma)}}, $$

for the function reconstructions in the matrix and in the fractures, and

$$ \mathrm{Err}_g = \frac{\|\nabla u - \nabla_{D_m} u_D\|_{L^2(\Omega)^d}}{\|\nabla u\|_{L^2(\Omega)^d}} + \frac{\|\nabla \gamma u - \nabla_{D_f} u_D\|_{L^2(\Gamma)^{d-1}}}{\|\nabla \gamma u\|_{L^2(\Gamma)^{d-1}}}. $$

for the gradient reconstructions in the matrix and in the fractures. These errors are reported for both schemes in Figure 6 for the first test case, and in Figure 7 for the second test case as function of the number of d.o.f. after elimination of the cell and Dirichlet unknowns. The convergence rate between two successive meshes $k$ and $k+1$ is also provided and computed as follows:

$$ \mathrm{CR}_u^{k+1} = d \frac{\ln \left( \frac{\|\mathrm{Err}_u^k\|}{\|\mathrm{Err}_u^{k+1}\|} \right)}{\ln \left( \frac{\#(Nb_{cells}^{k+1})}{\#(Nb_{cells}^k)} \right)}, \quad \mathrm{CR}_g^{k+1} = d \frac{\ln \left( \frac{\|\mathrm{Err}_g^k\|}{\|\mathrm{Err}_g^{k+1}\|} \right)}{\ln \left( \frac{\#(Nb_{cells}^{k+1})}{\#(Nb_{cells}^k)} \right)}. $$

It is reported for both schemes in tables 4 and 5 for the first test case, and in table 6 for the second test case. A second order convergence rate is observed on Cartesian meshes for both the function and gradient reconstructions, which is a typical super convergence behavior on such meshes. On randomly perturbed Cartesian meshes, the VAG scheme exhibits a second order convergence rate for the function reconstructions and a first order convergence rate for the gradient reconstructions. Since on randomly perturbed Cartesian meshes the faces are no longer planar, the HFV scheme no longer converges as expected, at least for the gradient reconstructions. On tetrahedral meshes a second order of convergence is also obtained for the function reconstructions and a first order of convergence is noticed for the gradient reconstructions for both test cases. This observed second order of convergence for the function reconstructions is as usual better than the first order estimate given by Remark 4.1, while the observed first order of convergence for the gradient reconstructions confirms the estimate given by Remark 4.1. It is also clear on both test cases that the HFV scheme is much less robust w.r.t. anisotropy than the VAG scheme, especially on tetrahedral meshes.

In all test cases, the linear system obtained after elimination of the cell and Dirichlet unknowns is solved using the GMRes iterative solver with the stopping criteria $10^{-10}$ and a maximum Krylov subspace dimension fixed to 1000 (not attained in our tests). The GMRes solver is preconditioned by ILUT [26], [27] 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 4, 5 and 6, we report the number of GMRes iterations $It$, the fill-in factor $F$ of the ILUT factorization defined as the ratio between the number of nonzero elements of the ILUT factorization by the number of nonzero elements of the matrix. We also report the CPU time taking into account the elimination of the cell and Dirichlet unknowns, the ILUT factorization, the GMRes iterations, and the computation of the cell values. It can be noticed that, on topologically Cartesian meshes, the CPU time is roughly speaking 4 times larger for the HFV scheme than for the VAG scheme. This large difference is not due to the number of nonzero elements in the matrices which is only slightly larger for the HFV scheme than for the VAG scheme. As can be checked in table 4, this difference is due to a larger number of GMRes iterations and to a higher fill-in factor of the ILUT factorization for the HFV scheme
than for the VAG scheme. On tetrahedral meshes, the CPU time for the computation of the HFV solution is larger of a factor from 10 to 20 than the CPU time obtained with the VAG scheme. This is due to a larger number of GMRes iterations, a larger fill-in factor for ILUT combined with a 5 times larger number of nonzero elements in the matrix.

<figure><img src="image_9.png" /><figcaption>Figure 6: First test case. For the 3 families of meshes (top: Cartesian meshes, middle: randomly perturbed Cartesian meshes, and bottom: tetrahedral meshes), and for the isotropic (left) and anisotropic (right) subcases: sum of $L^2$ norm of the relative error in the matrix and in the fracture for the function and its gradients reconstruction both for VAG and HFV schemes as the function of the number of d.o.f. (after elimination of the cell and Dirichlet unknowns).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7: Second test case. For the 10 tetrahedral meshes and for the isotropic (left) and anisotropic (right) subcases: sum of $L^2$ norm of the relative error in the matrix and in the fracture for the function and its gradients reconstruction both for VAG and HFV schemes as the function of the number of d.o.f. (after elimination of the cell and Dirichlet unknowns).</figcaption></figure>

<table><thead><tr><th colspan="10">Isotropic case, Cartesian</th><th colspan="10">Anisotropic case, Cartesian</th></tr><tr><th colspan="10">Vertex Approximate Gradient Discretization</th><th colspan="10">Hybrid Finite Volume Discretization</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr></thead><tbody><tr><td>1</td><td>3</td><td>1.2</td><td>4.03·10<sup>-2</sup></td><td>0.11</td><td>n/a</td><td>n/a</td><td>9.6·10<sup>-4</sup></td><td>3</td><td>1.2</td><td>3.8·10<sup>-2</sup></td><td>8.5·10<sup>-2</sup></td><td>n/a</td><td>n/a</td><td>7.7·10<sup>-4</sup></td></tr><tr><td>2</td><td>5</td><td>2.1</td><td>1.08·10<sup>-2</sup></td><td>0.03</td><td>1.89</td><td>1.62</td><td>6.9·10<sup>-3</sup></td><td>5</td><td>1.9</td><td>8.7·10<sup>-3</sup></td><td>2.3·10<sup>-2</sup></td><td>2.12</td><td>1.87</td><td>6.1·10<sup>-3</sup></td></tr><tr><td>3</td><td>9</td><td>2.4</td><td>2.9·10<sup>-3</sup></td><td>1.04·10<sup>-2</sup></td><td>1.92</td><td>1.71</td><td>9.9·10<sup>-2</sup></td><td>9</td><td>2.2</td><td>2.1·10<sup>-3</sup></td><td>6.6·10<sup>-3</sup></td><td>2.06</td><td>1.83</td><td>8.1·10<sup>-2</sup></td></tr><tr><td>4</td><td>16</td><td>2.5</td><td>7.4·10<sup>-4</sup></td><td>3.03·10<sup>-3</sup></td><td>1.95</td><td>1.77</td><td>0.89</td><td>14</td><td>2.1</td><td>5.1·10<sup>-4</sup></td><td>1.9·10<sup>-3</sup></td><td>2.03</td><td>1.81</td><td>0.71</td></tr><tr><td>5</td><td>30</td><td>2.5</td><td>1.9·10<sup>-4</sup></td><td>8.5·10<sup>-4</sup></td><td>1.97</td><td>1.82</td><td>8.4</td><td>20</td><td>2.2</td><td>1.3·10<sup>-4</sup></td><td>5.3·10<sup>-4</sup></td><td>2.02</td><td>1.84</td><td>5.5</td></tr><tr><td>6</td><td>56</td><td>2.5</td><td>4.8·10<sup>-5</sup></td><td>2.3·10<sup>-4</sup></td><td>1.89</td><td>1.86</td><td>90</td><td>29</td><td>2.2</td><td>3.2·10<sup>-5</sup></td><td>1.4·10<sup>-4</sup></td><td>2.01</td><td>1.87</td><td>48</td></tr><tr><th colspan="10">Hybrid Finite Volume Discretization</th><th colspan="10">Isotropic case, Perturbated Cartesian</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr><tr><td>1</td><td>6</td><td>3.3</td><td>1.08·10<sup>-2</sup></td><td>3.8·10<sup>-2</sup></td><td>n/a</td><td>n/a</td><td>1.3·10<sup>-3</sup></td><td>4</td><td>2.5</td><td>1.3·10<sup>-2</sup></td><td>0.16</td><td>n/a</td><td>n/a</td><td>9.7·10<sup>-4</sup></td></tr><tr><td>2</td><td>10</td><td>3.6</td><td>2.7·10<sup>-3</sup></td><td>1.2·10<sup>-2</sup></td><td>1.98</td><td>1.64</td><td>1.2·10<sup>-2</sup></td><td>6</td><td>3.1</td><td>3.3·10<sup>-3</sup></td><td>5.8·10<sup>-2</sup></td><td>1.97</td><td>1.23</td><td>7.5·10<sup>-3</sup></td></tr><tr><td>3</td><td>17</td><td>3.6</td><td>6.9·10<sup>-4</sup></td><td>3.4·10<sup>-3</sup></td><td>1.99</td><td>1.71</td><td>0.14</td><td>10</td><td>3.6</td><td>8.2·10<sup>-4</sup></td><td>1.9·10<sup>-2</sup></td><td>1.99</td><td>1.49</td><td>9.4·10<sup>-2</sup></td></tr><tr><td>4</td><td>29</td><td>3.6</td><td>1.7·10<sup>-4</sup></td><td>9.6·10<sup>-4</sup></td><td>1.99</td><td>1.77</td><td>1.4</td><td>18</td><td>3.7</td><td>2.1·10<sup>-4</sup></td><td>6.1·10<sup>-3</sup></td><td>2.01</td><td>1.58</td><td>1.4</td></tr><tr><td>5</td><td>59</td><td>3.6</td><td>4.3·10<sup>-5</sup></td><td>2.6·10<sup>-4</sup></td><td>2</td><td>1.82</td><td>17.7</td><td>30</td><td>3.8</td><td>5.3·10<sup>-5</sup></td><td>1.8·10<sup>-3</sup></td><td>1.99</td><td>1.68</td><td>18</td></tr><tr><td>6</td><td>122</td><td>3.6</td><td>1.1·10<sup>-5</sup></td><td>7.1·10<sup>-5</sup></td><td>2</td><td>1.86</td><td>313</td><td>65</td><td>3.8</td><td>1.3·10<sup>-5</sup></td><td>5.1·10<sup>-4</sup></td><td>1.99</td><td>1.78</td><td>303</td></tr></tbody></table>

<table><caption>Table 4: First test case. For the isotropic (left) and anisotropic (right) subcases, the VAG and HFV schemes and the six Cartesian meshes (above) and the six perturbated Cartesian meshes (below): mesh number Nb, number IT of GMRES iterations preconditioned by ILUT, fill-in factor F, sum of the $L^2$ relative errors in the matrix and in the fractures for the function (Err$_u$) and for the gradient reconstruction (Err$_g$), convergence rates for the function (CR$_u$) and for the gradient (CR$_g$) reconstruction, CPU time in seconds.</caption><thead><tr><th colspan="10">Isotropic case, Perturbated Cartesian</th><th colspan="10">Anisotropic case, Perturbated Cartesian</th></tr><tr><th colspan="10">Vertex Approximate Gradient Discretization</th><th colspan="10">Hybrid Finite Volume Discretization</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr></thead><tbody><tr><td>1</td><td>4</td><td>1.2</td><td>4.03·10<sup>-2</sup></td><td>0.11</td><td>n/a</td><td>n/a</td><td>7.8·10<sup>-4</sup></td><td>3</td><td>1.2</td><td>3.8·10<sup>-2</sup></td><td>8.6·10<sup>-2</sup></td><td>n/a</td><td>n/a</td><td>7.8·10<sup>-4</sup></td></tr><tr><td>2</td><td>5</td><td>2.1</td><td>1.08·10<sup>-2</sup></td><td>3.4·10<sup>-2</sup></td><td>1.92</td><td>1.34</td><td>6.3·10<sup>-3</sup></td><td>5</td><td>1.9</td><td>8.7·10<sup>-3</sup></td><td>2.3·10<sup>-2</sup></td><td>2.14</td><td>1.44</td><td>5.8·10<sup>-3</sup></td></tr><tr><td>3</td><td>9</td><td>2.4</td><td>2.8·10<sup>-3</sup></td><td>1.03·10<sup>-2</sup></td><td>1.92</td><td>1.69</td><td>9.8·10<sup>-2</sup></td><td>8</td><td>2.2</td><td>2.1·10<sup>-3</sup></td><td>6.6·10<sup>-3</sup></td><td>2.02</td><td>1.61</td><td>7.6·10<sup>-2</sup></td></tr><tr><td>4</td><td>16</td><td>2.5</td><td>7.4·10<sup>-4</sup></td><td>3.02·10<sup>-3</sup></td><td>2.06</td><td>-0.11</td><td>0.908</td><td>14</td><td>2.1</td><td>5.1·10<sup>-4</sup></td><td>1.9·10<sup>-3</sup></td><td>2.16</td><td>-0.43</td><td>0.64</td></tr><tr><td>5</td><td>29</td><td>2.5</td><td>1.9·10<sup>-4</sup></td><td>8.5·10<sup>-4</sup></td><td>1.87</td><td>0.91</td><td>8.4</td><td>20</td><td>2.3</td><td>1.3·10<sup>-4</sup></td><td>5.2·10<sup>-4</sup></td><td>1.85</td><td>0.97</td><td>5.6</td></tr><tr><td>6</td><td>56</td><td>2.5</td><td>4.7·10<sup>-5</sup></td><td>2.3·10<sup>-4</sup></td><td>1.96</td><td>2.11</td><td>85</td><td>29</td><td>2.3</td><td>3.2·10<sup>-5</sup></td><td>1.4·10<sup>-4</sup></td><td>1.95</td><td>1.81</td><td>48</td></tr><tr><th colspan="10">Hybrid Finite Volume Discretization</th><th colspan="10">Isotropic case, Perturbated Cartesian</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>Err<sub>u</sub></th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr><tr><td>1</td><td>5</td><td>3.3</td><td>1.1·10<sup>-2</sup></td><td>7.4·10<sup>-2</sup></td><td>n/a</td><td>n/a</td><td>1.5·10<sup>-3</sup></td><td>5</td><td>2.9</td><td>1.5·10<sup>-2</sup></td><td>0.67</td><td>n/a</td><td>n/a</td><td>1.4·10<sup>-3</sup></td></tr><tr><td>2</td><td>9</td><td>3.6</td><td>2.8·10<sup>-3</sup></td><td>3.4·10<sup>-2</sup></td><td>1.99</td><td>1.12</td><td>1.6·10<sup>-2</sup></td><td>8</td><td>3.4</td><td>5.2·10<sup>-3</sup></td><td>0.48</td><td>1.52</td><td>0.48</td><td>1.7·10<sup>-2</sup></td></tr><tr><td>3</td><td>16</td><td>3.6</td><td>6.9·10<sup>-4</sup></td><td>1.1·10<sup>-2</sup></td><td>2.01</td><td>1.58</td><td>0.15</td><td>11</td><td>3.7</td><td>1.3·10<sup>-3</sup></td><td>0.17</td><td>2.01</td><td>1.47</td><td>0.17</td></tr><tr><td>4</td><td>27</td><td>3.6</td><td>2.7·10<sup>-4</sup></td><td>7.3·10<sup>-2</sup></td><td>1.35</td><td>-2.71</td><td>1.6</td><td>16</td><td>3.8</td><td>5.7·10<sup>-3</sup></td><td>1.38</td><td>-2.14</td><td>-3.01</td><td>1.9</td></tr><tr><td>5</td><td>50</td><td>3.6</td><td>7.8·10<sup>-5</sup></td><td>3.8·10<sup>-2</sup></td><td>1.79</td><td>0.96</td><td>20</td><td>31</td><td>3.9</td><td>2.2·10<sup>-3</sup></td><td>0.85</td><td>1.39</td><td>0.71</td><td>33</td></tr><tr><td>6</td><td>106</td><td>3.6</td><td>1.3·10<sup>-5</sup></td><td>9.5·10<sup>-3</sup></td><td>2.60</td><td>1.97</td><td>292</td><td>57</td><td>3.9</td><td>1.9·10<sup>-4</sup></td><td>0.21</td><td>3.53</td><td>2.03</td><td>356</td></tr></tbody></table>
<table><caption>Table 5: First test case. For the isotropic (left) and anisotropic (right) subcases, the VAG and HFV schemes and the ten tetrahedral meshes: mesh number Nb, number IT of GMRES iterations preconditioned by ILUT, fill-in factor F, sum of the L² relative errors in the matrix and in the fractures for the function (Err_u) and for the gradient reconstruction (Err_g), convergence rates for the function (CR_u) and for the gradient (CR_g) reconstruction, CPU time in seconds.</caption><thead><tr><th colspan="10">Isotropic case</th><th colspan="5">Anisotropic case</th></tr><tr><th colspan="10">Vertex Approximate Gradient Discretization</th><th colspan="5">Anisotropic case</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr></thead><tbody><tr><td>1</td><td>5</td><td>2.1</td><td>9.1·10<sup>-3</sup></td><td>0.22</td><td>n/a</td><td>n/a</td><td>3.6·10<sup>-3</sup></td><td>5</td><td>2.1</td><td>5.7·10<sup>-3</sup></td><td>0.21</td><td>n/a</td><td>n/a</td><td>3.8·10<sup>-3</sup></td></tr><tr><td>2</td><td>8</td><td>2.6</td><td>2.3·10<sup>-3</sup></td><td>0.11</td><td>1.96</td><td>1.08</td><td>4.9·10<sup>-2</sup></td><td>9</td><td>2.6</td><td>1.5·10<sup>-3</sup></td><td>0.11</td><td>1.94</td><td>1.04</td><td>4.5·10<sup>-2</sup></td></tr><tr><td>3</td><td>16</td><td>2.8</td><td>5.1·10<sup>-4</sup></td><td>4.6·10<sup>-2</sup></td><td>2.02</td><td>1.05</td><td>0.61</td><td>15</td><td>2.9</td><td>3.5·10<sup>-4</sup></td><td>4.8·10<sup>-2</sup></td><td>1.91</td><td>1.02</td><td>0.68</td></tr><tr><td>4</td><td>20</td><td>2.9</td><td>3.1·10<sup>-4</sup></td><td>3.6·10<sup>-2</sup></td><td>2.02</td><td>1.05</td><td>1.5</td><td>18</td><td>2.9</td><td>2.1·10<sup>-4</sup></td><td>3.7·10<sup>-2</sup></td><td>2.09</td><td>1.02</td><td>1.7</td></tr><tr><td>5</td><td>25</td><td>2.9</td><td>2.1·10<sup>-4</sup></td><td>2.8·10<sup>-2</sup></td><td>1.69</td><td>1.06</td><td>2.9</td><td>22</td><td>3</td><td>1.5·10<sup>-4</sup></td><td>2.9·10<sup>-2</sup></td><td>1.39</td><td>1.04</td><td>3.8</td></tr><tr><td>6</td><td>29</td><td>2.9</td><td>1.3·10<sup>-4</sup></td><td>2.3·10<sup>-2</sup></td><td>2.37</td><td>1.01</td><td>5.9</td><td>25</td><td>3</td><td>9.1·10<sup>-5</sup></td><td>2.4·10<sup>-2</sup></td><td>2.47</td><td>0.99</td><td>6.9</td></tr><tr><td>7</td><td>33</td><td>2.9</td><td>9.6·10<sup>-5</sup></td><td>2.1·10<sup>-2</sup></td><td>2.21</td><td>1.02</td><td>10</td><td>26</td><td>3</td><td>6.7·10<sup>-5</sup></td><td>2.1·10<sup>-2</sup></td><td>2.33</td><td>1.01</td><td>12</td></tr><tr><td>8</td><td>37</td><td>2.9</td><td>7.9·10<sup>-5</sup></td><td>1.8·10<sup>-2</sup></td><td>1.81</td><td>0.99</td><td>14</td><td>30</td><td>3</td><td>5.7·10<sup>-5</sup></td><td>1.9·10<sup>-2</sup></td><td>1.66</td><td>0.98</td><td>16</td></tr><tr><td>9</td><td>40</td><td>2.9</td><td>6.9·10<sup>-5</sup></td><td>1.7·10<sup>-2</sup></td><td>1.81</td><td>1.05</td><td>18</td><td>30</td><td>3</td><td>4.9·10<sup>-5</sup></td><td>1.7·10<sup>-2</sup></td><td>1.62</td><td>1.05</td><td>22</td></tr><tr><td>10</td><td>49</td><td>2.9</td><td>4.3·10<sup>-5</sup></td><td>1.3·10<sup>-2</sup></td><td>2.01</td><td>1.01</td><td>41</td><td>37</td><td>3</td><td>3.1·10<sup>-5</sup></td><td>1.4·10<sup>-2</sup></td><td>1.98</td><td>1.01</td><td>48</td></tr></tbody><tbody><tr><th colspan="10">Hybrid Finite Volume Discretization</th><th colspan="5">Anisotropic case</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr><tr><td>1</td><td>10</td><td>4.4</td><td>1.5·10<sup>-3</sup></td><td>0.12</td><td>n/a</td><td>n/a</td><td>1.7·10<sup>-2</sup></td><td>15</td><td>4.6</td><td>3.3·10<sup>-2</sup></td><td>1.7</td><td>n/a</td><td>n/a</td><td>2.7·10<sup>-2</sup></td></tr><tr><td>2</td><td>17</td><td>4.7</td><td>3.9·10<sup>-4</sup></td><td>6.4·10<sup>-2</sup></td><td>1.97</td><td>0.92</td><td>0.23</td><td>25</td><td>5.1</td><td>1.1·10<sup>-2</sup></td><td>1.1</td><td>1.72</td><td>0.72</td><td>0.49</td></tr><tr><td>3</td><td>33</td><td>4.8</td><td>9.8·10<sup>-5</sup></td><td>3.1·10<sup>-2</sup></td><td>1.87</td><td>0.98</td><td>3.3</td><td>39</td><td>5.2</td><td>2.6·10<sup>-3</sup></td><td>0.56</td><td>1.83</td><td>0.85</td><td>7.2</td></tr><tr><td>4</td><td>42</td><td>4.9</td><td>5.7·10<sup>-5</sup></td><td>2.4·10<sup>-2</sup></td><td>2.08</td><td>0.99</td><td>8.9</td><td>48</td><td>5.2</td><td>1.5·10<sup>-3</sup></td><td>0.44</td><td>1.90</td><td>0.92</td><td>18</td></tr><tr><td>5</td><td>52</td><td>4.9</td><td>4.3·10<sup>-5</sup></td><td>1.9·10<sup>-2</sup></td><td>1.25</td><td>0.97</td><td>22</td><td>62</td><td>5.2</td><td>1.1·10<sup>-3</sup></td><td>0.36</td><td>1.86</td><td>0.92</td><td>40</td></tr><tr><td>6</td><td>69</td><td>4.9</td><td>2.5·10<sup>-5</sup></td><td>1.5·10<sup>-2</sup></td><td>2.73</td><td>1.03</td><td>47</td><td>82</td><td>5.2</td><td>6.8·10<sup>-4</sup></td><td>0.29</td><td>2.03</td><td>0.98</td><td>87</td></tr><tr><td>7</td><td>78</td><td>4.9</td><td>1.8·10<sup>-5</sup></td><td>1.3·10<sup>-2</sup></td><td>2.42</td><td>1.05</td><td>74</td><td>94</td><td>5.2</td><td>5.3·10<sup>-4</sup></td><td>0.25</td><td>1.96</td><td>0.97</td><td>141</td></tr><tr><td>8</td><td>86</td><td>4.9</td><td>1.5·10<sup>-5</sup></td><td>1.2·10<sup>-2</sup></td><td>1.62</td><td>0.92</td><td>106</td><td>112</td><td>5.2</td><td>4.3·10<sup>-4</sup></td><td>0.23</td><td>1.94</td><td>0.94</td><td>203</td></tr><tr><td>9</td><td>94</td><td>4.9</td><td>1.3·10<sup>-5</sup></td><td>1.1·10<sup>-2</sup></td><td>1.60</td><td>0.98</td><td>154</td><td>125</td><td>5.2</td><td>3.7·10<sup>-4</sup></td><td>0.21</td><td>1.95</td><td>0.97</td><td>285</td></tr><tr><td>10</td><td>122</td><td>4.9</td><td>8.3·10<sup>-6</sup></td><td>8.9·10<sup>-2</sup></td><td>2.03</td><td>0.99</td><td>346</td><td>192</td><td>5.2</td><td>2.3·10<sup>-4</sup></td><td>0.17</td><td>1.97</td><td>0.98</td><td>725</td></tr></tbody></table>

<table><caption>Table 6: Second test case. For the isotropic (left) and the anisotropic (right) subcases, the VAG and HFV schemes and the ten tetrahedral meshes: mesh number Nb, number IT of GMRES iterations preconditioned by ILUT, fill-in factor F, sum of the L² relative errors in the matrix and in the fractures for the function (Err_u) and for the gradient reconstruction (Err_g), convergence rates for the function (CR_u) and for the gradient (CR_g) reconstruction, CPU time in seconds.</caption><thead><tr><th colspan="10">Isotropic case</th><th colspan="5">Anisotropic case</th></tr><tr><th colspan="10">Vertex Approximate Gradient Discretization</th><th colspan="5">Anisotropic case</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr></thead><tbody><tr><td>1</td><td>11</td><td>2.5</td><td>4.1·10<sup>-3</sup></td><td>0.61</td><td>n/a</td><td>n/a</td><td>4.3·10<sup>-2</sup></td><td>13</td><td>2.5</td><td>1.6·10<sup>-3</sup></td><td>0.42</td><td>n/a</td><td>n/a</td><td>4.8·10<sup>-2</sup></td></tr><tr><td>2</td><td>32</td><td>2.7</td><td>1.5·10<sup>-3</sup></td><td>0.38</td><td>1.96</td><td>0.91</td><td>0.34</td><td>20</td><td>2.8</td><td>7.2·10<sup>-4</sup></td><td>0.26</td><td>1.52</td><td>0.91</td><td>0.31</td></tr><tr><td>3</td><td>51</td><td>2.8</td><td>8.6·10<sup>-4</sup></td><td>0.29</td><td>2.28</td><td>1.06</td><td>0.77</td><td>52</td><td>2.8</td><td>4.7·10<sup>-4</sup></td><td>0.21</td><td>1.60</td><td>0.92</td><td>0.83</td></tr><tr><td>4</td><td>51</td><td>2.8</td><td>5.7·10<sup>-4</sup></td><td>0.23</td><td>1.94</td><td>1.05</td><td>1.5</td><td>51</td><td>2.9</td><td>3.5·10<sup>-4</sup></td><td>0.17</td><td>1.38</td><td>0.89</td><td>1.7</td></tr><tr><td>5</td><td>51</td><td>2.9</td><td>3.8·10<sup>-4</sup></td><td>0.18</td><td>2.07</td><td>1.04</td><td>3.1</td><td>51</td><td>2.9</td><td>2.7·10<sup>-4</sup></td><td>0.14</td><td>1.35</td><td>0.89</td><td>3.4</td></tr><tr><td>6</td><td>51</td><td>2.9</td><td>2.2·10<sup>-4</sup></td><td>0.13</td><td>2.08</td><td>1.06</td><td>7.9</td><td>51</td><td>2.9</td><td>1.8·10<sup>-4</sup></td><td>0.11</td><td>1.53</td><td>0.98</td><td>7.9</td></tr><tr><td>7</td><td>51</td><td>2.9</td><td>1.8·10<sup>-4</sup></td><td>0.12</td><td>1.88</td><td>0.97</td><td>9.8</td><td>51</td><td>2.9</td><td>1.6·10<sup>-4</sup></td><td>9.9·10<sup>-2</sup></td><td>1.37</td><td>0.88</td><td>11</td></tr><tr><td>8</td><td>52</td><td>2.9</td><td>1.5·10<sup>-4</sup></td><td>0.11</td><td>2.12</td><td>1.12</td><td>14</td><td>53</td><td>2.9</td><td>1.3·10<sup>-4</sup></td><td>8.8·10<sup>-2</sup></td><td>1.73</td><td>1.08</td><td>18</td></tr><tr><td>9</td><td>60</td><td>2.9</td><td>1.3·10<sup>-4</sup></td><td>9.7·10<sup>-2</sup></td><td>1.95</td><td>1.03</td><td>19</td><td>56</td><td>2.9</td><td>1.2·10<sup>-4</sup></td><td>8.3·10<sup>-2</sup></td><td>1.41</td><td>0.95</td><td>22</td></tr><tr><td>10</td><td>78</td><td>2.9</td><td>8.2·10<sup>-5</sup></td><td>7.4·10<sup>-2</sup></td><td>2.04</td><td>1.03</td><td>40</td><td>72</td><td>3</td><td>8.4·10<sup>-5</sup></td><td>6.7·10<sup>-2</sup></td><td>1.48</td><td>0.97</td><td>48</td></tr></tbody><tbody><tr><th colspan="10">Hybrid Finite Volume Discretization</th><th colspan="5">Anisotropic case</th></tr><tr><th>Nb</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th><th>It</th><th>F</th><th>ErrF</th><th>Err<sub>g</sub></th><th>CR<sub>u</sub></th><th>CR<sub>g</sub></th><th>CPU</th></tr><tr><td>1</td><td>51</td><td>4.7</td><td>1.08·10<sup>-3</sup></td><td>0.43</td><td>n/a</td><td>n/a</td><td>0.32</td><td>52</td><td>5.1</td><td>1.1·10<sup>-2</sup></td><td>8.6</td><td>n/a</td><td>n/a</td><td>0.54</td></tr><tr><td>2</td><td>51</td><td>4.8</td><td>3.7·10<sup>-4</sup></td><td>0.22</td><td>1.91</td><td>1.02</td><td>1.9</td><td>84</td><td>5.2</td><td>3.4·10<sup>-3</sup></td><td>4.4</td><td>2.09</td><td>1.09</td><td>3.9</td></tr><tr><td>3</td><td>53</td><td>4.8</td><td>2.3·10<sup>-4</sup></td><td>0.16</td><td>1.96</td><td>1.09</td><td>4.8</td><td>97</td><td>5.2</td><td>2.2·10<sup>-3</sup></td><td>3.7</td><td>1.73</td><td>0.89</td><td>9.8</td></tr><tr><td>4</td><td>71</td><td>4.9</td><td>1.4·10<sup>-4</sup></td><td>0.12</td><td>2.07</td><td>1.08</td><td>11</td><td>108</td><td>5.2</td><td>1.5·10<sup>-3</sup></td><td>3.1</td><td>1.68</td><td>0.86</td><td>22</td></tr><tr><td>5</td><td>88</td><td>4.9</td><td>9.1·10<sup>-5</sup></td><td>9.5·10<sup>-2</sup></td><td>2.07</td><td>1.07</td><td>25</td><td>146</td><td>5.2</td><td>1.1·10<sup>-3</sup></td><td>2.6</td><td>1.84</td><td>0.91</td><td>50</td></tr><tr><td>6</td><td>114</td><td>4.9</td><td>5.1·10<sup>-5</sup></td><td>6.9·10<sup>-2</sup></td><td>2.05</td><td>1.06</td><td>67</td><td>248</td><td>5.2</td><td>6.5·10<sup>-4</sup></td><td>2.1</td><td>1.80</td><td>0.89</td><td>155</td></tr><tr><td>7</td><td>132</td><td>4.9</td><td>4.1·10<sup>-5</sup></td><td>6.1·10<sup>-2</sup></td><td>1.98</td><td>1.01</td><td>104</td><td>532</td><td>5.2</td><td>5.4·10<sup>-4</sup></td><td>1.9</td><td>1.72</td><td>0.84</td><td>387</td></tr><tr><td>8</td><td>146</td><td>4.9</td><td>3.2·10<sup>-5</sup></td><td>5.3·10<sup>-2</sup></td><td>2.06</td><td>1.05</td><td>162</td><td>530</td><td>5.2</td><td>4.4·10<sup>-4</sup></td><td>1.7</td><td>1.83</td><td>0.88</td><td>565</td></tr><tr><td>9</td><td>165</td><td>4.9</td><td>2.8·10<sup>-5</sup></td><td>4.9·10<sup>-2</sup></td><td>1.99</td><td>1.01</td><td>216</td><td>312</td><td>5.2</td><td>3.9·10<sup>-4</sup></td><td>1.6</td><td>1.80</td><td>0.90</td><td>477</td></tr><tr><td>10</td><td>196</td><td>4.9</td><td>1.8·10<sup>-6</sup></td><td>3.8·10<sup>-2</sup></td><td>2.02</td><td>1.03</td><td>498</td><td>748</td><td>5.2</td><td>2.6·10<sup>-4</sup></td><td>1.3</td><td>1.89</td><td>0.94</td><td>906</td></tr></tbody></table>

# 7 Conclusion

In this paper, the gradient scheme framework [16] is extended to hybrid dimensional Darcy flow models in fractured porous media. Both the Vertex Approximate Gradient and the Hybrid Finite Volume schemes are shown to satisfy, whatever the choice of the control volumes, the coercivity, consistency, limit-conformity and compacity assumptions of the gradient scheme framework. These properties ensures in particular the convergence of the schemes to a weak solution of the model. One of the key ingredient to prove that both schemes satisfy this framework
is the density of smooth function subspaces for both the solution and flux Hilbert spaces. This result is obtained for a general network of polygonal fractures including intersecting, immersed or non immersed fractures. The numerical experiments carried out on Cartesian, hexahedral and tetrahedral families of meshes exhibit the convergence of both schemes except as expected for the HFV scheme with non planar faces. The results also clearly show that the VAG scheme is much cheaper in terms of CPU time than the HFV scheme on tetrahedral meshes and is also more robust regarding anisotropy.

**Acknowledgements** The authors would like to thank GDFSuez EP and Storengy for partially supporting this work, and Robert Eymard for fruitful discussions during the elaboration of this work.

## References

[1] Adams, R.A., Sobolev Spaces, Pure and Applied Mathematics, Academic Press, Vol. 140, New York, (1978).

[2] Ahmed, R., Edwards, M.G., Lamine, S., Huisman, B.A.H., Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model, J. Comp. Physics, 462-489, Vol. 284 (2015).

[3] Alboin, C., Jaffré, J., Roberts, J., Serres, C.: Modeling fractures as interfaces for flow and transport in porous media. Fluid flow and transport in porous media 295, 13-24 (2002).

[4] Angot, P., Boyer, F., Hubert, F.: Asymptotic and numerical modelling of flows in fractured porous media. ESAIM Mathematical Modelling and Numerical Analysis 23, 239-275 (2009).

[5] Brenner, K., Masson, R.: Convergence of a Vertex centered Discretization of Two-Phase Darcy flows on General Meshes. Int. Journal of Finite Volume Methods, June (2013).

[6] Brenner, K., Groza, M., Guichard, C., Masson, R.: Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media. ESAIM Mathematical Modelling and Numerical Analysis, 49, 303-330 (2015).

[7] Brezis, H.: Functional Analysis, Sobolev Spaces and Partial Differential Equations, Springer, Universitext, 2011.

[8] Brezzi, F., Lipnikov, K., Simoncini V.: A family of mimetic finite difference methods on polygonal and polyhedral meshes, Mathematical Models and Methods in Applied Sciences, vol. 15, 10, 2005, 1533-1552.

[9] Costabel, M.: Boundary integral operators on Lipschitz domains: elementary results. SIAM J. Math. Anal., 19, 613-626 (1988).

[10] D'Angelo, C., Scotti, A.: A mixed finite element method for Darcy flow in fractured porous media with non-matching grids. ESAIM Mathematical Modelling and Numerical Analysis 46,2, 465-489 (2012).
[11] Droniou, J., Eymard, R., Gallouët, T., Guichard, C., Herbin, R.: Gradient schemes for elliptic and parabolic problems. Personal Communication (2014).

[12] Droniou, J., Eymard, R., Gallouët, T., Herbin, R.: A Unified Approach to Mimetic Finite Difference, Hybrid Finite Volume and Mixed Finite Volume Methods. Math. Models and Methods in Appl. Sci. 20, 2, 265-295 (2010).

[13] Eymard, R., Gallouët, T., Herbin, R.: Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilisation and hybrid interfaces. IMA J. Numer. Anal. 30, 4, 1009-1043 (2010).

[14] Eymard, R., Guichard, C., Herbin, R.: Small-stencil 3D schemes for diffusive flows in porous media. Mathematical Modelling and Numerical Analysis 46, 265-290 (2012).

[15] Eymard, R., Herbin, R., Guichard, C., Masson, R.: Vertex Centered discretization of compositional Multiphase Darcy flows on general meshes. Comp. Geosciences 16, 4, 987-1005 (2012).

[16] Droniou, J., Eymard, R., Gallouët, T., Herbin, R.: Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations. Math. Models Methods Appl. Sci. 23, 13, 2395-2432 (2013).

[17] Flauraud, E., Nataf, F., Faille, I., Masson, R.: Domain Decomposition for an asymptotic geological fault modeling, Comptes Rendus à l'académie des Sciences de Mécanique 331, 849-855 (2003).

[18] Grisvard, P.: Elliptic Problems on non smooth domains. Pitman Publishing Inc., Marshfield, Massachusetts (1985).

[19] Karimi-Fard, M., Durlofsky, L.J., Aziz, K.: An efficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE journal, june (2004).

[20] Formaggia, L., Fumagalli, A., Scotti, A., Ruffo, P.: A reduced model for Darcy's problem in networks of fractures. ESAIM Mathematical Modelling and Numerical Analysis 48, 4, 1089-1116 (2014).

[21] Martin, V., Jaffré, J., Roberts, J.E.: Modeling fractures and barriers as interfaces for flow in porous media. SIAM J. Sci. Comput. 26, 5, 1667-1691 (2005).

[22] McLean, W.: Strongly Elliptic Systems and Boundary Integral Equations. Cambridge University Press, Cambridge, UK, (2000).

[23] Mikhailov, S. E.: Traces, extensions and co-normal derivatives for elliptic systems on Lipschitz domains. J. Math. Anal. Appl. 378, 1, 324-342 (2011)

[24] Monteagudu, J., Firoozabadi, A.: Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects. SPE Journal 12, 355-366 (2007).

[25] Reichenberger, V., Jakobs, H., Bastian, P., Helmig, R.: A mixed-dimensional finite volume method for multiphase flow in fractured porous media. Adv. Water Resources 29, 7, 1020-1036 (2006).
[26] Saad, Y.: Iterative Methods for Sparse Linear Systems. 2nd edition, SIAM, Philadelphia, PA, (2003)

[27] Saad, Y.: http://www-users.cs.umn.edu/ saad/software/

[28] Sandve, T.H., Berre, I., Nordbotten, J.M.: An efficient Multi-Point Flux Approximation Method for discrete Fracture-Matrix Simulations. J. Comp. Phys. 231, 3784-3800 (2012).

[29] Si, H.: http://tetgen.org

[30] Tunc, X., Faille, I., Gallouët, T., Cacas, M.C., Havé, P.: A model for conductive faults with non matching grids. Comp. Geosciences 16, 277-296 (2012).