# Mathematical and numerical analysis of the mixed formulation of single phase flow in three-dimensional fractured porous media

Michel Kern, Géraldine Pichot, Daniel Zegarra Vasquez

► To cite this version:

Michel Kern, Géraldine Pichot, Daniel Zegarra Vasquez. Mathematical and numerical analysis of the mixed formulation of single phase flow in three-dimensional fractured porous media. 2025. hal-05029638v2

HAL Id: hal-05029638

https://inria.hal.science/hal-05029638v2

Preprint submitted on 10 Oct 2025

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
# Mathematical and numerical analysis of the mixed formulation of single phase flow in three-dimensional fractured porous media

Michel Kern*†

Géraldine Pichot*†

Daniel Zegarra Vasquez*†

October 10, 2025

## Abstract

Three-dimensional fractured-porous media are typically modeled using a reduced approach known as the Discrete Fracture Matrix (DFM) model. In this model, the rock matrix is represented in three dimensions, while the fracture network is represented as a two-dimensional (co-dimension one) structure. The system is governed by a set of partial differential equations (PDEs) that couple the flow in the 3D rock matrix with the flow in the 2D fracture network. Mixed finite element methods are particularly well-suited for discretizing this system because they ensure local mass conservation. This paper presents an existence and uniqueness result for the PDEs in mixed form, assuming that the hydraulic head is continuous across fractures (i.e., the fractures are conductive). Additionally, error estimates for the mixed finite element discretization are derived under a regularity assumption, and a mixed-hybrid approximation for the coupled problem is developed. The theoretical findings are validated through comparison with an analytical solution.

## Keywords

single-phase flow, fractured porous media, discrete fracture-matrix model, $H(\text{div})$ spaces, conforming meshes, mixed-hybrid finite element method, numerical analysis.

## Introduction

Fractures are ubiquitous in the subsurface and play a key role in many physical and chemical phenomena. Simulating fluid flow in fractured rocks is essential in fields such as energy and environment, including oil and gas water extraction, geothermal energy production, nuclear waste storage, CO2 storage, groundwater resources as main examples. Imaging the subsurface is a challenge as measurements (boreholes, outcrops, lineament maps) provide data at some specific scales and there are a lot of uncertainties. The Discrete Fracture Matrix (DFM) model [9] provides an explicit representation of fractures by modeling them as lower-dimensional objects within a three-dimensional rock matrix. Fractures are considered two-dimensional flat surfaces (polygons) within this 3D space. A broader overview of fracture representation can be found in [9].

The model of flow in fractured porous media, valid for an arbitrary network of fractures, is by now well-understood, see for instance [9, 11, 15, 6]. Darcy-like averaged equations govern flow in both the rock matrix and the fracture network, with a coupling term that describes exchanges between the two media [57]. The fracture network is used to define a set of subdomains, and the coupling conditions are set on the interface between the subdomains (a fracture, or a passive interface), and on the fracture intersections. The reader is referred to [2, 55, 5, 13] for a more detailed discussion on the derivation of the mathematical model. More recently a more general derivation of mixed dimensional models from the initial full 3D equations has been carried out under different asymptotic regimes, see [50, 54, 42, 19].

In this work, the hydraulic head is assumed to be continuous across fractures, as in [2, 15]. We emphasize that this is a modeling choice, relevant when the fractures are considered highly

<footnote>*Inria, 48 rue Barrault, 75013 Paris, France.</footnote>

<footnote>†CERMICS, École nationale des ponts et chaussées, 8 Avenue Blaise Pascal, 77420 Champs-sur-Marne, France.</footnote>

<footnote>E-mails: michel.kern@inria.fr, geraldine.pichot@inria.fr, daniel.zegarra-vasquez@inria.fr.</footnote>
conductive. Another possible choice is the discontinuous hydraulic head model studied in [55, 3, 13], which is more appropriate when one considers impermeable fractures that act as barriers. The mathematical analysis for the discontinuous hydraulic head model cannot be used in the present setting, in large part because the hydraulic head jump leads to an additional unknown on each side of the fracture, and an additional equation for the normal fluxes on each side, granting control on the flux jumps, something that is missing in the model analyzed here.

A large variety of numerical methods have been applied to the simulation of flow in fracture networks, and in fractured porous media: several variants of finite volumes [5, 61, 56, 48, 17, 15, 16], Mixed and Mixed-Hybrid Finite Element [55, 4, 13, 29, 27, 44, 66, 59, 60, 35], Discontinuous Galerkin [6, 53, 52, 67], Mimetic Finite Difference [45, 58], Virtual Finite Element [7], Hybrid High-Order [20, 33], Extended Finite Element [26, 34, 37, 36] and others [47, 10, 64, 11, 41]. A comparison of several methods on test cases of varying difficulty can be found in [8] (the “Field” test case is used for validation in the companion paper [49]).

In this work, as well as in the two companion papers [49, 46], we use a Mixed-Hybrid Finite Element (MHFE) method on a simplicial conforming mesh. The Mixed Finite Element method has the advantage of ensuring mass conservation, and of being usable on unstructured simplicial meshes. Finite Volume methods, which are also conservative, are in some cases cheaper in terms on number of unknowns, but the simplest TPFA method is severely limited in the types of geometry it can accommodate. Therefore, more sophisticated finite volume methods (some of which can accommodate unstructured meshes) have been developed and successfully applied to fractured porous flow (see the above references). The Mixed-Hybrid variant is valuable as it leads to a symmetric and positive definite system, with one unknown per face of the tetrahedral mesh and one unknown per edge of the triangular mesh.

The details of the mesh generation are reported in [49]. We briefly recall that a triangular mesh of the fracture network is generated first using a dedicated surface mesh generation software [51] and then the volume mesh is generated with the fracture mesh as a constraint. As a result the meshes are matching at the fracture-rock matrix interface. We note that it is possible to work with non-matching grids, as is done in several of the works cited previously [26, 34, 37, 36, 35, 59, 60, 13].

Since the discretization method is based on a mixed formulation, it is natural to express the original system of PDEs in its mixed form and then check the well-posedness of the resulting system. The analysis of the mixed formulation, and its approximation, is carried out in [13] at a similar level of geometric description, but within a more general mixed-dimensional framework that encompasses flow across all structural dimensions, ranging from cross-points at fracture-fracture intersections (0D), to fracture intersections (1D), the fracture network (2D), and the surrounding rock matrix (3D). However, [13] uses the discontinuous hydraulic head model, and this turns out to make a large difference as far as proving existence and uniqueness is concerned. In this work, by using the continuous hydraulic head approximation, we lose a coupling condition on the fluxes. In [13], this condition is critical in order to control the traces on the boundary of the subdomains. The existence and uniqueness proof relies on the so-called “inf-sup” condition. We were able to establish the validity of the inf-sup condition by making use of the analysis for a continuous hydraulic head model in primal form presented in [15]. A significant difference between [13] and [15] is that the former builds the solution of the auxiliary problem piecewise by subdomain in each dimension, while the latter gives a global solution. We noticed that by being able to work globally, we did not need to assume that each 3D subdomain has Dirichlet conditions on part of its boundary (as is done in [13]), as this assumption is quite restrictive.

The coupling conditions linking the equations in the rock matrix (3D) to those in the fracture network (2D) lead to some technical difficulties when we define the functional setting. These coupling conditions involve mainly traces of vector fields on parts of the boundaries (faces in 3D, edges in 2D). We were able to handle these coupling conditions by expressing that the sum of the fluxes vanishes on fracture-fracture intersections thanks to a weak form of these conditions proposed in [15] (see also [66]).
The discrete problem uses the lowest order Raviart-Thomas-Nédélec element. As usual, checking the uniform inf-sup condition is the most delicate part of the numerical analysis. Once again, the fracture-fracture intersections need to be treated carefully. This time we built on [66], where a careful construction of the $H(\text{div})$ conforming interpolation operator is analyzed. The definition of the interpolation operators (both in 2D and 3D) require more smoothness on the solution than what is given from the existence theorem. Because of the complex geometry, in particular with immersed fractures, we do not know whether such smoothness is always attainable. Furthermore, the error estimates require even more smoothness. In the concluding section 5, we indicate some possible directions that could let us weaken the regularity requirements.

The numerical solution of the discrete problem is based on the Mixed-Hybrid formulation, which is by now well known to be equivalent to the Mixed formulation, and has the advantage of leading to a symmetric and positive definite linear system. We carried out this step in detail to show the differences entailed by the mixed-dimensional coupled problem. Thus, we were able to show how one proceeds from the general Mixed-Hybrid formulation to the actual block linear system that has to be assembled and solved. The derivation clearly shows that one may first work with the 3D rock matrix part and the 2D fracture network part separately, provided one has identified the faces in the tetrahedral mesh that belong to the fractures and those that do not. In a second step, the blocks that express the coupling are added together.

An important additional step was to adapt to our purpose the alternative formulation of the Mixed-Hybrid system as a mesh-dependent variational formulation, with a bilinear form that can be computed from local information on each mesh element. We have followed the presentation in [23, 24], but note that this work was already present in [18] or [21]. This result plays an important role in the companion paper [46] for showing that the linear system arising from the Mixed-Hybrid formulation can be solved by using the recently developed domain decomposition method named GenEO [63, 28].

The main contributions of this paper are:

* a proof of existence and uniqueness for the Mixed formulation of the coupled problem for single-phase flow in fractured porous media;

* the numerical analysis of the discrete problem using Mixed Finite Element, as well as the derivation of the linear system from the associated Mixed-Hybrid formulation;

* a numerical convergence study that illustrates the theoretical result for an example from [15] with two intersecting fractures, where an analytical solution is available.

The outline of the paper is as follows. Section 1 presents the geometrical setting and the physical model. Section 2 presents an appropriate functional setting for the continuous problem and proves its well-posedness. Section 3 contains the numerical analysis of the Mixed Finite Element approximation. A convergence result is proved, and the Mixed-Hybrid approximation for the coupled problem is developed. Section 4 illustrates the theoretical results by studying the numerical convergence of the MHFE method thanks to an exact solution in a simple DFM geometry. Conclusions and perspectives for further work are discussed in Section 5.

# 1 Physical problem description

## 1.1 Geometrical setting

We denote by $\Omega$ the domain representing a permeable fractured porous rock. The fracture network $\Omega^f$ embedded in $\bar{\Omega}$ is modeled as a collection of $N^f$ fractures represented as polygons: $\Omega^f = \cup_{\ell=1}^{N^f} \Omega_\ell^f$. The rock matrix domain is defined as $\Omega^m = \Omega \setminus \bar{\Omega}^f$. We can find an illustration in Figure 1a.
For the purpose of the mathematical description, we divide each polygon $\Omega_{\ell}^{\mathfrak{f}}$, $\ell = 1, \dots, N^{\mathfrak{f}}$, into sub-polygons $\gamma_{\ell}$ along each intersection line, such that $\Omega^{\mathfrak{f}} = \cup_{\ell \in \mathcal{I}} \gamma_{\ell}$, where $\mathcal{I}$ is the index set of the sub-polygons. The sub-polygons are built so as to satisfy: if $\ell \neq k$, then $\gamma_{\ell} \cap \gamma_{k} = \emptyset$; if $\ell \neq k$, then $\bar{\gamma}_{\ell} \cap \bar{\gamma}_{k}$ is either an edge or a point or an empty set, as described in [66, 6]. An example of decomposition in sub-polygons, based on a mesh of the fractures, is shown on Figure 1b.

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) A cubic fractured porous medium $\Omega$ with $N^{\mathfrak{f}} = 3$ fractures; (b) Example of a decomposition of the fracture network in sub-polygons (one color per sub-polygons).</figcaption></figure>

The sub-polygons $\gamma_{\ell}$, $\ell \in \mathcal{I}$, further partition the rock matrix domain $\Omega^{\mathfrak{m}}$ into connected components $\omega_{\alpha}$ such that $\Omega^{\mathfrak{m}} = \cup_{\alpha \in \mathcal{C}} \omega_{\alpha}$, where $\mathcal{C}$ is the index set of connected components. The subdomains $\omega_{\alpha}$ may also be further subdivided to satisfy: if $\alpha \neq \beta$, then $\omega_{\alpha} \cap \omega_{\beta} = \emptyset$; if $\alpha \neq \beta$, then $\bar{\omega}_{\alpha} \cap \bar{\omega}_{\beta}$ is either a face or an edge or an empty set. Let $\mathcal{I}_{\alpha} = \{\ell, \gamma_{\ell} \subset \partial\omega_{\alpha}\}$ the set of sub-polygons on $\partial\omega_{\alpha}$.

We further assume that both the sub-polygons and the subdomains are Lipschitz (see also [6]). This is always achievable in practice. For example, after generating a mesh for both $\Omega^{\mathfrak{f}}$ and $\Omega^{\mathfrak{m}}$, the sub-polygons and the subdomains may be built by agglomeration of mesh elements.

Dirichlet and Neumann boundary conditions are set on the rock matrix domain on $\Gamma^D$ and $\Gamma^N$ respectively, with $\bar{\partial}\Omega^{\mathfrak{m}} = \Gamma_0$ where $\Gamma_0 = \Gamma^D \cup \Gamma^N$. Let $\Gamma_{\alpha} = \partial\omega_{\alpha} = \Gamma_{\alpha}^D \cup \Gamma_{\alpha}^N \cup \Gamma_{\alpha}^{\gamma} \cup \Gamma_{\alpha}^{\text{fict}}$, with $\Gamma_{\alpha}^D = \partial\omega_{\alpha} \cap \Gamma^D$, $\Gamma_{\alpha}^N = \partial\omega_{\alpha} \cap \Gamma^N$, $\Gamma_{\alpha}^{\gamma} = \cup_{\ell \in \mathcal{I}_{\alpha}} \gamma_{\ell}$ and $\Gamma_{\alpha}^{\text{fict}}$ which gathers the faces introduced by the subdivision process (they will only play a role in Sections 1.2 and A). Let $\Gamma^{\text{fict}}$ be the set of all such “fictitious” faces (each one belongs to two subdomains).

We assume that Dirichlet and Neumann boundary conditions also hold for the fracture network on the cube sides: $\bar{\partial}\Omega^{\mathfrak{f}} \cap \bar{\partial}\Omega = \Sigma^D \cup \Sigma^N$. For each $\ell \in \mathcal{I}$, $\Sigma_{\ell}^{D} = \partial\gamma_{\ell} \cap \Sigma^D$ and $\Sigma_{\ell}^{N} = \partial\gamma_{\ell} \cap \Sigma^N$. We denote by $\Sigma_{\ell,0} = \Sigma_{\ell}^{D} \cup \Sigma_{\ell}^{N}$ the segments shared between the subdomain $\gamma_{\ell}$ and the domain $\Omega$ and $\Sigma_0 = \cup_{\ell \in \mathcal{I}} \Sigma_{\ell,0}$.

The intersections are segments shared by two or more sub-polygons. Let $\Sigma$ be the set of intersections between the sub-polygons, and $\mathcal{I}_{\sigma} = \{\ell \in \mathcal{I}, \sigma \subset \Sigma \cap \partial\gamma_{\ell}\}$ the index set of polygons that contain $\sigma$ as intersection. For each $\ell \in \mathcal{I}$, we let $\Sigma_{\ell}^{\text{inter}}$ be the set of intersections of the sub-polygon $\gamma_{\ell}$.

In this study, homogeneous Neumann boundary conditions will be applied at the fracture tips, defined as the borders of the sub-polygons $\gamma_{\ell}$, $\ell \in \mathcal{I}$ that are not intersections nor on the border of the cube. The fracture tips are part of $\Sigma^N$. For each $\ell \in \mathcal{I}$, we define $\Sigma_{\ell} = \Sigma_{\ell}^{D} \cup \Sigma_{\ell}^{N} \cup \Sigma_{\ell}^{\text{inter}}$.

## 1.2 Physical model

We consider a single-phase, incompressible, and time-independent flow problem, as detailed for instance in [55, 15, 16, 8, 9, 3, 11]. The unknowns are:

* the fluid flow velocity in the rock matrix domain $u^{\mathfrak{m}} = (u_{\alpha}^{\mathfrak{m}})_{\alpha \in \mathcal{C}}$, obtained from the fluid velocity $u_{\alpha}^{\mathfrak{m}}$ computed in each subdomain $\omega_{\alpha}$, in m/s,
* the fluid flow velocity in the fracture network $u^f = (u_\ell^f)_{\ell \in \mathcal{I}}$, obtained from the fluid velocity $u_\ell^f$ computed in each sub-polygon $\gamma_\ell$, in m$^2$/s,

* the hydraulic head in the rock matrix $p^m = (p_\alpha^m)_{\alpha \in \mathcal{C}}$ obtained from the hydraulic head $p_\alpha^m$ computed in each subdomain $\omega_\alpha$, in m,

* the hydraulic head in the fracture network $p^f = (p_\ell^f)_{\ell \in \mathcal{I}}$ obtained from the hydraulic head $p_\ell^f$ computed in each subdomain $\gamma_\ell$, in m.

In the rock matrix, the hydraulic conductivity is a tensor denoted by $K^m = (K_\alpha^m)_{\alpha \in \mathcal{C}}$. In the fracture network, the tangential hydraulic conductivity is also a tensor denoted by $\tilde{K}^f = (\tilde{K}_\ell^f)_{\ell \in \mathcal{I}}$. Their physical physical unit is m/s. Furthermore, in the fracture network, the aperture is a scalar denoted by $d_\ell^f = (d_\ell^f)_{\ell \in \mathcal{I}}$, whose physical unit is m. Thus, in the fracture network, one can define the transmissivity as a tensor: $K^f = (K_\ell^f)_{\ell \in \mathcal{I}} = (d_\ell^f \tilde{K}_\ell^f)_{\ell \in \mathcal{I}}$, whose unit is m$^2$/s. In the sequel, we only work with $K^f$.

The source term in the rock matrix (resp. in the fracture network) is denoted by $f^m = (f_\alpha^m)_{\alpha \in \mathcal{C}}$ (resp. $f^f = (f_\ell^f)_{\ell \in \mathcal{I}}$), whose unit is 1/s (resp. m/s).

The Neumann condition in the rock matrix (resp. in the fracture network) is denoted by $q^{m,N} = (q_\alpha^{m,N})_{\alpha \in \mathcal{C}}$ (resp. $q^{f,N} = (q_\ell^{f,N})_{\ell \in \mathcal{I}}$), whose unit is m/s (resp. m$^2$/s).

The Dirichlet condition in the rock matrix (resp. in the fracture network) is denoted by $g^{m,D} = (g_\alpha^{m,D})_{\alpha \in \mathcal{C}}$ (resp. $g^{f,D} = (g_\ell^{f,D})_{\ell \in \mathcal{I}}$), whose unit is m (resp. m).

For a smooth function $p^m$ (resp. $p_\alpha^m$) defined in $\Omega^m$ (resp. a subdomain $\omega_\alpha$), we let $\gamma^{m,\partial\Omega^m} p^m$ (resp. $\gamma^{m,\partial\omega_\alpha} p_\alpha^m$) denote the restriction of $p^m$ to the boundary of $\Omega^m$ (resp. $\omega_\alpha$).

For a smooth vector field $u^m$ defined in $\Omega^m$, we let $\gamma_n^{m,\partial\Omega^m} u^m$ denote the normal trace of $u^m$ along (a part of) the boundary of $\Omega^m$. Informally, $\gamma_n^{m,\partial\Omega^m} u^m = u^m \cdot n|_{\partial\Omega^m}$, where $n$ is the unit outward normal to $\partial\Omega^m$. If we want to specify the subdomain $\omega_\alpha$ over which $u_\alpha^m$ is defined, we will write $\gamma_n^{m,\partial\omega_\alpha} u_\alpha^m$. If a part of the boundary is an immersed sub-polygon $\gamma_\ell$, with $\ell \in \mathcal{I}$, then $\gamma_\ell$ is shared by two subdomains, namely $\omega_\alpha^+$ and $\omega_\alpha^-$. We denote by $n_\ell^\pm$ the unit outward normal to $\partial\omega_\alpha^\pm$, and by $u_\alpha^\pm$ the solution inside the subdomain $\omega_\alpha^\pm$. The normal trace of $u_\alpha^\pm$ along $\gamma_\ell$ will be written as $\gamma_{n_\ell^\pm}^{m,\gamma_\ell} u_\alpha^\pm$, depending on the side ($\pm$) of the sub-polygon $\gamma_\ell$. Now, we define the jump across $\gamma_\ell$ as:

$$ [[u^m]]|_{\gamma_\ell} = \gamma_{n_\ell^+}^{m,\gamma_\ell} u_{\alpha_\ell^+}^m + \gamma_{n_\ell^-}^{m,\gamma_\ell} u_{\alpha_\ell^-}^m. \quad (1.1) $$

Similarly, for a scalar function $p^f$ or a (tangential) vector field $u^f$ defined on $\Omega^f$ (resp. $p_\ell^f$ or $u_\ell^f$ defined over a sub-polygon $\gamma_\ell$), we denote by $\gamma^{f,\partial\Omega^f} p^f$ and $\gamma_n^{f,\partial\Omega^f} u^f$ (resp. $\gamma^{f,\partial\gamma_\ell} p_\ell^f$ and $\gamma_n^{f,\partial\gamma_\ell} u_\ell^f$) its (normal) trace along the boundary of $\Omega^f$ (resp. of $\gamma_\ell$). If we are looking at the trace along a particular intersection segment $\sigma$, we will denote the (normal) trace by $\gamma^{f,\sigma} p_\ell^f$ and $\gamma_n^{f,\sigma} u_\ell^f$.

In each rock matrix subdomain $\omega_\alpha$, $\alpha \in \mathcal{C}$, the flow problem is written as:

$$ \nabla \cdot u_{\alpha}^{m} = f_{\alpha}^{m} \quad \text{in } \omega_{\alpha}, \qquad (1.2a) $$

$$ u_{\alpha}^{m} = -K_{\alpha}^{m} \nabla p_{\alpha}^{m} \quad \text{in } \omega_{\alpha}, \qquad (1.2b) $$

$$ \gamma^{m,\Gamma_{\alpha}^{D}} p_{\alpha}^{m} = g_{\alpha}^{m,D} \quad \text{on } \Gamma_{\alpha}^{D}, \qquad (1.2c) $$

$$ \gamma_{n}^{m,\Gamma_{\alpha}^{N}} u_{\alpha}^{m} = q_{\alpha}^{m,N} \quad \text{on } \Gamma_{\alpha}^{N}. \qquad (1.2d) $$

Additionally, the hydraulic head is continuous and the jump of the flux is zero across a fictitious face shared by two subdomains.
In each sub-polygon $\gamma_\ell$, $\ell \in \mathcal{I}$, the flow problem is written as:

$$ \nabla_{\tau} \cdot u_{\ell}^{\mathfrak{f}} = [[u^{\mathfrak{m}}]]_{|\gamma_{\ell}} + f_{\ell}^{\mathfrak{f}} \quad \text{in } \gamma_{\ell}, \qquad (1.3a) $$

$$ u_{\ell}^{\mathfrak{f}} = -\mathcal{K}_{\ell}^{\mathfrak{f}} \nabla_{\tau} p_{\ell}^{\mathfrak{f}} \quad \text{in } \gamma_{\ell}, \qquad (1.3b) $$

$$ \gamma_{\ell}^{\mathfrak{f}, \Sigma_{\ell}^{D}} p_{\ell}^{\mathfrak{f}} = g_{\ell}^{\mathfrak{f}, D} \quad \text{on } \Sigma_{\ell}^{D}, \qquad (1.3c) $$

$$ \gamma_{n}^{\mathfrak{f}, \Sigma_{\ell}^{N}} u_{\ell}^{\mathfrak{f}} = q_{\ell}^{\mathfrak{f}, N} \quad \text{on } \Sigma_{\ell}^{N}, \qquad (1.3d) $$

where we have denoted the surface divergence by $\nabla_{\tau} \cdot$. Since the fractures are assumed to be planar, $\nabla_{\tau} \cdot$ is simply the two-dimensional divergence, once a coordinate system has been chosen for each fracture. The exchange term $[[u^{\mathfrak{m}}]]_{|\gamma_{\ell}}$ that couples the flow between the rock matrix subdomains $\omega_{\alpha+}$ and $\omega_{\alpha-}$ and the fracture $\gamma_{\ell}$ acts as a source term for the fracture flow problem.

We also assume the continuity of the hydraulic head between the subdomains and the sub-polygons (relevant for modeling conductive fractures):

$$ \gamma^{\mathfrak{m}, \gamma_{\ell}} p_{\alpha_{\ell}}^{\mathfrak{m}} = \gamma^{\mathfrak{m}, \gamma_{\ell}} p_{\alpha_{\ell}}^{\mathfrak{m}} = p_{\ell}^{\mathfrak{f}}. \qquad (1.4) $$

We also add conditions at each fracture-fracture intersection $\sigma \in \Sigma$ [66]:

$$ \gamma^{\mathfrak{f}, \sigma} p_{\ell}^{\mathfrak{f}}|_{\sigma} = \gamma^{\mathfrak{f}, \sigma} p_{k}^{\mathfrak{f}}|_{\sigma} \quad \text{for } (\ell, k) \in \mathcal{I}_{\sigma}^{2}, \qquad (1.5a) $$

$$ \sum_{\ell \in \mathcal{I}_{\sigma}} \gamma_{n}^{\mathfrak{f}, \sigma} u_{\ell}^{\mathfrak{f}} = 0, \quad \text{for each } \sigma \subset \Sigma. \qquad (1.5b) $$

Condition (1.5a) represents the continuity of the hydraulic head $p^{\mathfrak{f}}$ for each intersection $\sigma \subset \Sigma$. Condition (1.5b) represents the mass balance of the fracture fluid's flux per unit length across each intersection $\sigma \subset \Sigma$ between fractures.

We assume that $\mathcal{K}^{\mathfrak{m}}(x)$ (resp. $\mathcal{K}^{\mathfrak{f}}(x)$) is a symmetric and positive definite tensor at each point $x \in \Omega^{\mathfrak{m}}$ (resp. $x \in \Omega^{\mathfrak{f}}$) and that there exists constants $0 < K_{\min}^{\mathfrak{m}}, K_{\max}^{\mathfrak{m}} < \infty$ (resp. $0 < K_{\min}^{\mathfrak{f}}, K_{\max}^{\mathfrak{f}} < \infty$) such that:

$$ K_{\min}^{\mathfrak{m}} \|\xi\|_{2}^{2} \le \xi^{T} \mathcal{K}^{\mathfrak{m}}(x) \xi \le K_{\max}^{\mathfrak{m}} \|\xi\|_{2}^{2}, \quad \forall \xi \in \mathbb{R}^{3}, \forall x \in \Omega^{\mathfrak{m}}, $$

$$ K_{\min}^{\mathfrak{f}} \|\xi\|_{2}^{2} \le \xi^{T} \mathcal{K}^{\mathfrak{f}}(x) \xi \le K_{\max}^{\mathfrak{f}} \|\xi\|_{2}^{2}, \quad \forall \xi \in \mathbb{R}^{3}, \forall x \in \Omega^{\mathfrak{f}}. $$

We further assume that $f^{\mathfrak{m}} \in L^2(\Omega^{\mathfrak{m}})$, that $f^{\mathfrak{f}} \in L^2(\Omega^{\mathfrak{f}})$, that $g_{\alpha}^{\mathfrak{m}, D} \in H^{1/2}(\Gamma_{\alpha}^{D})$, and that $g_{\ell}^{\mathfrak{f}, D} \in H^{1/2}(\Sigma_{\ell}^{D})$.

For the sake of simplicity, we assume that $\forall \alpha \in \mathcal{C}, q_{\alpha}^{\mathfrak{m}} = 0$ on $\Gamma_{\alpha}^{N}$ and that $\forall \ell \in \mathcal{I}, q_{\ell}^{\mathfrak{f}} = 0$ on $\Sigma_{\ell}^{N}$. Otherwise, we could perform a lifting of the Neumann boundary conditions.

# 2 Continuous setting : existence and uniqueness

## 2.1 Complements on the trace of $H(\text{div})$ functions

As shown in the informal derivation of the weak derivation in Appendix A, we need to work with traces of functions over part of the boundary, for instance when looking at $[[u^{\mathfrak{m}}]]_{|\gamma_{\ell}}$ (the jump of $u^{\mathfrak{m}}$ across a fracture $\gamma_{\ell}$), or if we need to consider mixed boundary conditions. When the functions involved are only in an $H(\text{div})$ space, this presents some technical difficulties. In this section, we remind the reader how to tackle those difficulties thanks to a duality technique, following [5, 38] (see also [40]).

We now consider a domain $\mathcal{O} \subset \mathbb{R}^d$ ($d = 2$ or 3), assumed to be polygonal in 2D (resp. polyhedral in 3D), with a boundary denoted by $\partial \mathcal{O}$. We fix an edge (resp. face) denoted by $\Sigma$. We let $\gamma^{\Sigma}: H^1(\mathcal{O}) \to H_{00}^{1/2}(\Sigma)$ be the trace mapping on $\Sigma$. Then one may define a continuous
lifting $R_Σ : H_{00}^{1/2}(Σ) \to H^1(\mathcal{O})$ that is a right inverse of the trace mapping (see for instance [5] or [40]):

$$ \gamma^{\Sigma}(R_{\Sigma}g) = g, \quad \forall g \in H_{00}^{1/2}(\Sigma). $$

We may now define the trace of a vector field $v$ in $H(\mathrm{div}, \mathcal{O})$ by duality. The mapping $\gamma_n^\Sigma$ defined by:

$$ \langle \gamma_n^\Sigma v, g \rangle_\Sigma = \int_{\mathcal{O}} v \cdot \nabla R_\Sigma g + \mathrm{div} v R_\Sigma g, \quad \forall g \in H_{00}^{1/2}(\Sigma) $$

is continuous over $H_{00}^{1/2}(\Sigma)$, and thus defines an element of $H_{00}^{-1/2}(\Sigma)$, its dual.

## 2.2 Definition of the function spaces and forms

We denote $u = (u^m, u^f)$, and $v = (v^m, v^f)$, and similarly $p = (p^m, p^f)$ and $s = (s^m, s^f)$.

The definition of the functional spaces $V$ for the velocity and $S$ for the hydraulic head needs to satisfy two main constraints:

* $V$ must be a Hilbert space compatible with the application of the Ladyzhenskaya-Babuška-Brezzi theory;

* the “essential” boundary conditions ((1.2d) for the rock matrix, (1.3d) for the fracture network and (1.5b) at the fracture-fracture intersections) must be satisfied.

In order to define the jump of the flux across a fracture $\gamma_ℓ$, we apply the results from Section 2.1. As noted previously, the fracture is a boundary face of two subdomains $\omega_{\alpha_ℓ^+}$ and $\omega_{\alpha_ℓ^-}$. For a function $v^m \in L^2(\Omega^m)$, with $v_\alpha^m \in H(\mathrm{div}, \omega_\alpha)$, for all $\alpha \in \mathcal{C}$, the jump, defined as in (1.1)

$$ [[v^m]]|_{\gamma_\ell} = \gamma_{n_\ell^+}^{m, \gamma_\ell} v_{\alpha_\ell^+}^m + \gamma_{n_\ell^-}^{m, \gamma_\ell} v_{\alpha_\ell^-}^m, $$

is an element of $H_{00}^{-1/2}(\gamma_\ell)$.

For the hydraulic head, we choose the space $S = L^2(\Omega^m) \times L^2(\Omega^f)$.

For the velocity, we proceed as in [15] in order to properly handle the continuity condition (1.5b) on the intersections. We first define an auxiliary space $\tilde{V}$:

$$ \begin{gathered} \tilde{V} = \{v = (v^m, v^f) \in L^2(\Omega^m) \times L^2(\Omega^f) \mid v_\alpha^m \in H(\mathrm{div}, \omega_\alpha), \forall \alpha \in \mathcal{C}; \\ \gamma_n^{m, \Gamma_\alpha^N} v_\alpha^m = 0 \text{ on } \Gamma_\alpha^N, \forall \alpha \in \mathcal{C}; \\ (\nabla_\tau \cdot v_\ell^f - [[v^m]]|_{\gamma_\ell}) \in L^2(\gamma_\ell), \forall \ell \in \mathcal{I} \}. \end{gathered} \tag{2.1} $$

We endow $\tilde{V}$ with the following norm:

$$ \|v\|_{V}^{2} = \sum_{\alpha \in \mathcal{C}} \|v_{\alpha}^{m}\|_{H(\mathrm{div}, \omega_{\alpha})}^{2} + \|v^{f}\|_{L^{2}(\Omega^{f})}^{2} + \sum_{\ell \in \mathcal{I}} \|\nabla_{\tau} \cdot v_{\ell}^{f} - [[v^{m}]]|_{\gamma_{\ell}}\|_{L^{2}(\gamma_{\ell})}^{2}. \tag{2.2} $$

**Remark 2.1** (Traces of functions in $\tilde{V}$). In the definition of $\tilde{V}$ in (2.1), the “combined divergence” term $\nabla_\tau \cdot v_\ell^f - [[v^m]]|_{\gamma_\ell}$ is only square integrable as a whole. We have no information on the two terms $\nabla_\tau \cdot v_\ell^f$ and $[[v^m]]|_{\gamma_\ell}$ separately. In fact, the vector $v_\ell^f$ is not in $H(\mathrm{div}, \gamma_\ell)$, thus, it is not possible to define directly its trace on the boundary of $\gamma_\ell$. This is why we proceed as in [15] and will define the subspace $V$ below as a subset of $\tilde{V}$.

**Lemma 2.2.** $\tilde{V}$ is a Hilbert space for the norm defined in (2.2).

*Proof.* The proof follows [13, Lemma 2.2], see also [1].

We consider a Cauchy sequence $(v_k^m, v_k^f)_{k \in \mathbb{N}}$. Each one of the terms in the sum defined by (2.2) gives rise to a corresponding Cauchy sequence in the appropriate norm. Thus, for each $\alpha \in \mathcal{C}$, $(v_k^m)_\alpha$ is a Cauchy sequence in $H(\mathrm{div}, \omega_\alpha)$, so that there exists $v^m = (v_\alpha^m)_{\alpha \in \mathcal{C}}$ such that for all $\alpha \in \mathcal{C}$,
$\lim_{k \to +\infty} (v_k^m)_\alpha = v_\alpha^m \in H(\text{div}, \omega_\alpha)$. Also, by the continuity of the normal trace operators, the jump $[[v_k^m]]_{\gamma_\ell}$ converges in $H_{00}^{-1/2}(\gamma_\ell)$. Similarly, there exists $v_\ell^f = (v_\ell^f)_{\ell \in \mathcal{I}} \in L^2(\gamma_\ell)$ such that for all $\ell \in \mathcal{I}$, $\lim_{k \to +\infty} v_{k,\ell}^f = v_\ell^f$, but this time the convergence only takes place in $L^2(\gamma_\ell)$. Finally, the sequence $\eta_k = \nabla_\tau \cdot v_{k,\ell}^f - [[v_k^m]]_{|\gamma_\ell}$ converges in $L^2(\gamma_\ell)$ to a limit, say $\eta \in L^2(\gamma_\ell)$.

If we prove that $\eta = \nabla_\tau \cdot v_\ell^f - [[v^m]]_{|\gamma_\ell}$, we will have proved that $\tilde{V}$ is complete for the given norm (2.2). To do this, we take a function $\varphi \in C_0^\infty(\gamma_\ell)$, and we compute:

$$ \begin{aligned} \langle \eta, \varphi \rangle &= \lim_{k \to +\infty} \langle \eta_k, \varphi \rangle = \lim_{k \to +\infty} \langle \nabla_\tau \cdot v_{k,\ell}^f - [[v_k^m]]_{|\gamma_\ell}, \varphi \rangle = \lim_{k \to +\infty} \left( -\langle v_{k,\ell}^f, \nabla_\tau \varphi \rangle - \langle [[v_k^m]]_{|\gamma_\ell}, \varphi \rangle \right) \\ &= -\langle v_\ell^f, \nabla_\tau \varphi \rangle - \langle [[v^m]]_{|\gamma_\ell}, \varphi \rangle = \langle \nabla_\tau \cdot v_\ell^f - [[v^m]]_{|\gamma_\ell}, \varphi \rangle. \end{aligned} $$

As noted above, $\tilde{V}$ does not account for the boundary conditions for $v^f$. In order to properly take into account those boundary conditions, we follow [15] and impose the boundary conditions (1.3d) and (1.5b) in a weak form. We need to define a space for test functions, which will turn out to also be critical in the proof of the inf-sup condition later on (see Lemma 2.7 from Section 2.3). We define the following Hilbert space, where $\gamma$ denotes the trace operator in $H^1(\Omega)$:

$$ \tilde{H} = \{ s \in H^1(\Omega) \mid s^f := \gamma^{\Omega^f} s \in H^1(\Omega^f) \}. $$

We also define its closed subspace of functions with traces that match at fracture-fracture intersections:

$$ H = \{ s \in \tilde{H} \mid \gamma^{\mathcal{F},\sigma} s_j^f = \gamma^{\mathcal{F},\sigma} s_k^f, \forall \sigma \in \Sigma, \forall (j,k) \in \mathcal{I}_{\sigma}^2 \}. $$

We finally define:

$$ H_D = \{ s \in H \mid \gamma^{\mathfrak{m}, \Gamma_\alpha^D} s_\alpha^\mathfrak{m} = 0 \text{ on } \Gamma_\alpha^D, \forall \alpha \in \mathcal{C}, \gamma^{\mathfrak{f}, \Sigma_\ell^D} s_\ell^\mathfrak{f} = 0, \text{ on } \Sigma_\ell^D, \forall \ell \in \mathcal{I} \}. $$

The space $H$ is well defined in the case of immersed fractures, again as a consequence of [5]. More precisely, because we work with functions that are continuous across the sub-polygons, the traces of a function $s^\mathfrak{m} \in H^1(\Omega^\mathfrak{m})$ on the two sides of the sub-polygon $\gamma_\ell$ coincide, and we denote their common value by $s_\ell^\mathfrak{f}$ (see continuity condition (1.4) from Section 1.2).

We next define the bilinear form $d$ on $\tilde{V} \times H_D$:

$$ d(v, s) = \int_{\Omega^{\mathfrak{m}}} v^{\mathfrak{m}} \cdot \nabla s^{\mathfrak{m}} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} v_{\ell}^{\mathfrak{f}} \cdot \nabla_{\tau} s_{\ell}^{\mathfrak{f}} + \sum_{\alpha \in \mathcal{C}} \int_{J_{\omega_{\alpha}}} s_{\alpha}^{\mathfrak{m}} \nabla \cdot v_{\alpha}^{\mathfrak{m}} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} s_{\ell}^{\mathfrak{f}} \left( \nabla_{\tau} \cdot v_{\ell}^{\mathfrak{f}} - [[v^{\mathfrak{m}}]]_{|\gamma_{\ell}} \right). $$

We finally define the space for test functions:

$$ V = \{ v \in \tilde{V} \mid d(v, s) = 0, \forall s \in H_D \}. \quad (2.3) $$

**Remark 2.3.** The motivation for defining $V$ this way is to impose the boundary conditions (1.3d) and (1.5b) in a weak form. Indeed, we verify that, for functions $(v^m, v^f)$ with sufficient regularity (e.g., $v \in H(\text{div}, \Omega^m)$ with $\gamma_n^m(v_\alpha^m) \in L^2(\partial\omega_\alpha)$, $\forall \alpha \in \mathcal{C}$, and $v^f \in H(\text{div}, \gamma_\ell^f)$, $\forall \ell \in \mathcal{I}$), we have:

$$ \begin{aligned} d(v, s) &= \sum_{\alpha \in \mathcal{C}} \int_{\omega_\alpha} v^m \cdot \nabla s_\alpha^m + s_\alpha^m \nabla \cdot v_\alpha^m - \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} s_\ell^f \cdot [[v^m]]_{|\gamma_\ell} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} s_\ell^f \nabla_\tau \cdot v_\ell^f + v_\ell^f \cdot \nabla_\tau s_\ell^f \\ &= \sum_{\alpha \in \mathcal{C}} \int_{\omega_\alpha} \nabla \cdot (s_\alpha^m v_\alpha^m) - \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} s_\ell^f \cdot [[v^m]]_{|\gamma_\ell} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} \nabla_\tau \cdot (s_\ell^f v_\ell^f) \\ &= \sum_{\alpha \in \mathcal{C}} \int_{\partial\omega_\alpha} s_\ell^f \gamma_n^{\mathfrak{m}, \partial\omega_\alpha} v_\alpha^m - \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} s_\ell^f \cdot [[v^m]]_{|\gamma_\ell} + \sum_{\ell \in \mathcal{I}} \int_{\partial\gamma_\ell} \gamma^{\mathcal{F}, \partial\gamma_\ell} s_\ell^f \gamma_n^{\mathcal{F}, \partial\gamma_\ell} v_\ell^f \\ &= \sum_{\sigma \in \Sigma} \sum_{\ell \in \mathcal{I}_\sigma} \int_{\sigma} \gamma^{\mathcal{F}, \sigma} s_\ell^f \gamma_n^{\mathcal{F}, \sigma} v_\ell^f + \sum_{\ell \in \mathcal{I}} \int_{\Sigma_\ell^N} \gamma^{\mathcal{F}, \Sigma_\ell^N} s_\ell^f \gamma_n^{\mathcal{F}, \Sigma_\ell^N} v_\ell^f. \end{aligned} $$
**Lemma 2.4.** V is a Hilbert space for the norm defined in (2.2).

*Proof.* As a closed subspace of the Hilbert space $\tilde{V}$, V is itself a Hilbert space.

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

**Remark 2.5.** The definitions of the spaces H and V come from the paper [15]. In that paper, density results for (appropriately defined) smooth functions in the above space were given. We refer to [15] for details.

## 2.3 Weak mixed solution

We obtained the (formal) weak formulation for the flow problem in Section A. We now define the following bilinear forms:

$$a(u, v) = \int_{\Omega^m} (K^m)^{-1} u^m \cdot v^m + \int_{\Omega^f} (K^f)^{-1} u^f \cdot v^f, \quad \forall (u, v) \in V^2, \quad (2.4a)$$

$$b(u, s) = \sum_{\alpha \in \mathcal{C}} \int_{\omega_{\alpha}} s_{\alpha}^{\mathfrak{m}} \nabla \cdot u_{\alpha}^{\mathfrak{m}} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} s_{\ell}^{\mathfrak{f}} (\nabla_{\tau} \cdot u_{\ell}^{\mathfrak{f}} - [[u^{\mathfrak{m}}]]|_{\gamma_{\ell}}), \quad \forall u \in V, \forall s \in S. \quad (2.4b)$$

We also define the following linear forms:

$$L_a(v) = -\langle g^{\mathfrak{m},D}, \gamma_n^{\mathfrak{m}, \Gamma^D} v^m \rangle - \sum_{\ell \in \mathcal{I}} \langle g_{\ell}^{\mathfrak{f},D}, \gamma_n^{\mathfrak{f}, \Sigma_{\ell}^D} v_{\ell}^{\mathfrak{f}} \rangle, \quad \forall v \in V, \quad (2.5a)$$

$$L_b(s) = \int_{\Omega^m} f^{\mathfrak{m}} s^{\mathfrak{m}} + \int_{\Omega^f} f^{\mathfrak{f}} s^{\mathfrak{f}}, \quad \forall s \in S. \quad (2.5b)$$

A weak solution of problem (1.2)– (1.5) is a couple $(u, p) \in V \times S$ satisfying:

$$a(u, v) - b(v, p) = L_a(v), \quad \forall v \in V, \quad (2.6a)$$

$$b(u, s) = L_b(s), \quad \forall s \in S. \quad (2.6b)$$

We prove the existence and uniqueness for the mixed problem (2.6) by the usual Ladyzhenskaya-Babuška-Brezzi theory. On the one hand, we need to prove that the bilinear form $a$ is coercive on the kernel of $b$ (Lemma 2.6). On the other hand, we need to prove that $b$ satisfies the inf-sup condition (Lemma 2.7).

**Lemma 2.6 (Coercivity).** The bilinear form $a$ is continuous on $V$ and coercive on the space $W = \{v \in V | b(v, s) = 0, \forall s \in S\}$:

$$\exists \xi > 0 : a(v, v) \ge \xi ||v||_V^2, \forall v \in W. \quad (2.7)$$

*Proof.* The continuity of $a$ follows easily from the Cauchy-Schwarz inequality and the hypotheses on the permeability fields.

We first identify the subspace $W$, the kernel of $b$. From the definition of $b$ in (2.4b), we see that:

$$W = \{v = (v^{\mathfrak{m}}, v^{\mathfrak{f}}) \in V \mid \nabla \cdot v_{\alpha}^{\mathfrak{m}} = 0, \forall \alpha \in \mathcal{C}, \nabla_{\tau} \cdot v_{\ell}^{\mathfrak{f}} - [[v^{\mathfrak{m}}]]|_{\gamma_{\ell}} = 0, \forall \ell \in \mathcal{I}\}. \quad (2.8)$$

It is now clear that property (2.7) holds, thanks again to the hypotheses on the hydraulic conductivity fields.

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

**Lemma 2.7 (Inf-sup condition).** The bilinear form $b$ satisfies the inf-sup condition:

$$\exists \beta > 0 : \inf_{s \in S, s \neq 0} \sup_{v \in V, v \neq 0} \frac{b(v, s)}{||v||_V ||s||_S} \ge \beta. \quad (2.9)$$
*Proof.* We will prove the result by showing that there exists $C > 0$ such that:

Given $s = (s^m, s^\mathfrak{f}) \in S$, there exists $v = (v^m, v^\mathfrak{f})$ and a constant $C > 0$ such that $b(v, s) = \|s\|_S^2$, and that $\|v\|_V \le C \|s\|_S$. Indeed, assuming we have proven this claim, we would have:

$$ \frac{b(v, s)}{\|v\|_V \|s\|_S} \ge \frac{\|s\|_S^2}{C \|s\|_S^2} = \frac{1}{C}, $$

and the Lemma would follow.

Let $s = (s^m, s^\mathfrak{f}) \in S$. We obtain $v$ as the result of solving an auxiliary problem: Find $z \in H_D$ such that

$$ -\nabla \cdot (\nabla z) = s^m \qquad \text{in } \Omega^m, \qquad (2.10a) $$

$$ -\nabla_{\tau} \cdot (\nabla_{\tau}(\gamma^{m, \gamma_{\ell}} z)) + [[\nabla z]]|_{\gamma_{\ell}} = s^{\mathfrak{f}} \qquad \text{in } \gamma_{\ell}, \text{ for } \ell \in \mathcal{I}. \qquad (2.10b) $$

To obtain the weak formulation of (2.10), we first multiply (2.10a) by a test function $\varphi \in H_D$ and apply Green's formula. The contributions from the exterior boundaries vanish, and we regroup the terms over the sub-polygons:

$$ \int_{\Omega^m} \nabla z \cdot \nabla \varphi - \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} [[z]]|_{\gamma_\ell} \gamma^{m, \gamma_\ell} \varphi = \int_{\Omega^m} s^m \varphi. \qquad (2.11) $$

Then we turn to (2.10b) by multiplying by $\gamma^m \varphi$, using Green's formula over each sub-polygon $\gamma_\ell$ and adding the results. Again, the contributions from the exterior boundaries and the tips vanish, and we obtain:

$$ \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} \nabla_{\tau} \gamma^{m, \gamma_{\ell}} z \cdot \nabla_{\tau} \gamma^{m, \gamma_{\ell}} \varphi - \sum_{\sigma \in \Sigma} \int_{\sigma} \nabla_{\tau} \gamma^{m, \sigma} z \cdot \nu_{\sigma} \gamma^{\mathfrak{f}, \sigma} (\gamma^{\mathfrak{f}, \sigma} \varphi) + \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} [[z]]|_{\gamma_{\ell}} \gamma^{m, \gamma_{\ell}} \varphi = \int_{\Omega^{\mathfrak{f}}} s^{\mathfrak{f}} \gamma^{m, \partial \Omega^{\mathfrak{f}}} \varphi. \qquad (2.12) $$

Finally we add equations (2.11) and (2.12), so that the terms involving the trace jumps across the fractures cancel out and we obtain the (primal) weak formulation of (2.10):

$$ \int_{\Omega^m} \nabla z \cdot \nabla \varphi + \sum_{\ell \in \mathcal{I}} \int_{\gamma_\ell} \nabla_\tau \gamma^{m, \gamma_\ell} z \cdot \nabla_\tau \gamma^{m, \gamma_\ell} \varphi = \int_{\Omega^m} s^m \varphi + \int_{\Omega^\mathfrak{f}} s^\mathfrak{f} \gamma^{m, \partial \Omega^\mathfrak{f}} \varphi, \quad \forall \varphi \in H_D, \qquad (2.13) $$

Proposition 4 in [15] states that problem (2.13) has a unique solution $z \in H_D$ such that:

$$ \|z\|_{H_D} \le C \|s\|_S^2, $$

for some constant $C > 0$.

Thus, we set $v = (-\nabla z, -\nabla_\tau \gamma^m z)$. The already quoted result from [15] implies furthermore that the chosen $v$ belongs to $V$ (see its definition in (2.3)). In particular, $v^\mathfrak{f}$ satisfies the continuity condition (in the weak sense of (2.3)). We then see that the chosen $v$ satisfies:

$$ \begin{aligned} \nabla \cdot v^m &= s^m \quad \text{in } \Omega^m, \\ \nabla_\tau \cdot v^\mathfrak{f} - [[v^m]]|_{\gamma_\ell} &= s_\ell^\mathfrak{f} \quad \text{in } \gamma_\ell, \text{ for } \ell \in \mathcal{I}, \end{aligned} $$

so that $b(v, s) = \|s\|_S^2$.

Finally, we compute:

$$ \begin{aligned} \|v\|_V^2 &= \|v^m\|_{L^2(\Omega^m)}^2 + \sum_{\alpha \in \mathcal{C}} \|\nabla \cdot v_\alpha^m\|_{L^2(\omega_\alpha)}^2 + \sum_{\ell \in \mathcal{I}} \left( \|v_\ell^\mathfrak{f}\|_{L^2(\gamma_\ell)}^2 + \|\nabla_\tau \cdot v^\mathfrak{f} - [[v^m]]|_{\gamma_\ell}\|_{L^2(\gamma_\ell)}^2 \right) \\ &\le \|z\|_{H_D}^2 + \|s^m\|_{L^2(\Omega^m)}^2 + \|s^\mathfrak{f}\|_{L^2(\Omega^\mathfrak{f})}^2 \le (1+C) \|s\|_S^2, \end{aligned} $$

which shows that $v$ satisfies the required estimate, and the Lemma now follows as indicated at the beginning of the proof. ■
We now have all the results needed to state and prove the existence and uniqueness result.

**Theorem 2.8.** *Problem (2.6) has a unique solution $(u, p) \in V \times S$.*

*Proof.* The proof follows from Lemmas 2.6 and 2.7 and from the general existence theory for mixed problems. We refer the reader for instance to [12, 38, 32]. ■

**Remark 2.9.** *We note that the proof that the inf-sup condition (2.9) holds is different from those given in [13] or [3]. Those proofs worked piece-wise by dimension (increasing for [13], decreasing for [3]). This was not possible here because of the different transmission conditions for the hydraulic head. Because we use a continuous hydraulic head, it is not possible to control separately the divergence of the velocity on the fractures and the jump of the normal traces. Thus, we had to use a more global construction such as the one that was used in [15]. The difference between this work and the result proven in [15] (that we used in our proof) is that here we use a mixed formulation, whereas [15] studies a primal formulation. Of course, in both cases the formulation is linked to the numerical scheme. This work uses a (hybrid-)mixed finite element method, as we will see in Section 3, while [15] uses the Vertex Approximate Gradient finite volume scheme.*

# 3 Discrete setting : existence and uniqueness

In this section we present the discrete setting for approximating the problem (2.6).

## 3.1 Mixed finite element approximation

### 3.1.1 Meshes and discrete approximation spaces

We first describe the meshes used for the approximation. We saw in Section 1.1 that we have two distinct types of geometric entities: the rock matrix and the fracture network. Accordingly, we also need two separate meshes: one for the fracture network and another one for the rock matrix. In this work, we only use conforming meshes, in the sense that the matrix mesh is supported by the fracture mesh.

Accordingly, we describe the fracture mesh first. Let $K_h$ be a mesh of the fracture network, composed of $N(K_h)$ triangles, and let $\mathcal{E}_h$ be the set of edges of $K_h$, with the number of edges denoted by $N(\mathcal{E}_h)$. We define $\mathcal{E}_h^{\text{in}}$ to be the set of inner edges, (not intersection and not on the boundary), $\mathcal{E}_h^{\text{inter}}$ to be the set of intersection edges, $\mathcal{E}_h^{\text{D}}$ the set of Dirichlet edges and $\mathcal{E}_h^{\text{N}}$ the set of Neumann edges. The number of elements in each of these sets is denoted respectively by $N(\mathcal{E}_h^{\text{in}})$, $N(\mathcal{E}_h^{\text{inter}})$, $N(\mathcal{E}_h^{\text{D}})$ and $N(\mathcal{E}_h^{\text{N}})$. Let $\mathcal{E}_h \setminus \mathcal{E}_h^{\text{D}}$ be the set of non-Dirichlet edges of cardinal $N(\mathcal{E}_h \setminus \mathcal{E}_h^{\text{D}})$. We define $h^\mathcal{f} = \max_{K \in \mathcal{K}_h} h_K$, where $h_K$ is the diameter of $K$.

In the same way, we let $\mathcal{T}_h$ be a mesh of the rock matrix domain composed of $N(\mathcal{T}_h)$ tetrahedral elements. Let $\mathcal{F}_h$ be the set of faces of $\mathcal{T}_h$, whose number is denoted by $N(\mathcal{F}_h)$. We define $\mathcal{F}_h^{\text{in}}$ as the set of inner faces not on the boundary and not fractures (rock matrix faces), $\mathcal{F}_h^{\text{frac}}$ the set of tetrahedral faces that are fracture faces (triangles, see below), $\mathcal{F}_h^{\text{D}}$ the set of tetrahedral faces that are Dirichlet faces and $\mathcal{F}_h^{\text{N}}$ the set of Neumann faces. We have $\mathcal{F}_h = \mathcal{F}_h^{\text{in}} \cup \mathcal{F}_h^{\text{frac}} \cup \mathcal{F}_h^{\text{D}} \cup \mathcal{F}_h^{\text{N}}$. The number of inner faces is denoted by $N(\mathcal{F}_h^{\text{in}})$, the number of fracture faces is denoted by $N(\mathcal{F}_h^{\text{frac}})$, the number of Dirichlet and Neumann faces by $N(\mathcal{F}_h^{\text{D}})$ and $N(\mathcal{F}_h^{\text{N}})$ respectively and $N(\mathcal{F}_h \setminus \mathcal{F}_h^{\text{D}})$ is the number of non-Dirichlet faces. Let $\mathcal{F}_h^{\text{in,N}} = \mathcal{F}_h^{\text{in}} \cup \mathcal{F}_h^{\text{N}}$ be the set of faces that are inner or Neumann (neither Dirichlet nor fractures). We define $h^\mathcal{m} = \max_{T \in \mathcal{T}_h} h_T$, where $h_T$ is the diameter of $T$. Finally we let $h = \max(h^\mathcal{m}, h^\mathcal{f})$.

We assume that $K_h$ and $\mathcal{T}_h$ both satisfy the usual properties of a (simplicial) finite element mesh (see for instance [31, Sec. 8.2]). Additionally, we assume that both the rock matrix tetrahedral mesh $\mathcal{T}_h$ and the fracture network triangular mesh $K_h$ are regular, meaning that:

for $\mathcal{T}_h$ (resp. $K_h$) there exists $c^\mathcal{m} > 0$ (resp. $c^\mathcal{f} > 0$) such that $\frac{h_T}{\rho_T} \le c^\mathcal{m}$ (resp. $\frac{h_K}{\rho_K} \le c^\mathcal{f}$) for all $h > 0$ and for all $T \in \mathcal{T}_h$ (resp. $K \in K_h$), where $\rho_T$ (resp. $\rho_K$) is the diameter of the largest ball contained in $T$ (resp. $K$).
Finally, we assume that both meshes are built so as to accommodate the boundary conditions.

The approximation will use the lowest order Raviart-Thomas(-Nédélec) mixed finite element space for the velocity. We recall their definition:

### For the rock matrix

$$RTN^0(T) = \left\{ \varphi : \mathbb{R}^3 \to \mathbb{R}^3 \mid \forall x \in T, \varphi(x) = ax + (b \ c \ d)^T, (a, b, c, d) \in \mathbb{R}^4 \right\}, \forall T \in \mathcal{T}_h;$$

### For the fracture network

$$RT^0(K) = \left\{ \varphi : \mathbb{R}^2 \to \mathbb{R}^2 \mid \forall x \in K, \varphi(x) = ax + (b \ c)^T, (a, b, c) \in \mathbb{R}^3 \right\}, \forall K \in \mathcal{K}_h.$$

We recall two important properties of the Raviart-Thomas-Nédélec spaces:

* For the rock matrix, $\text{div}(RTN^0(T)) \subset P_0(\mathbb{R}^3)$;

* For the fracture network, $\text{div}(RT^0(K)) \subset P_0(\mathbb{R}^2)$;

* For $\varphi \in RTN^0(T)$, the normal trace $\gamma_n^{\text{m},\partial T}\varphi$ belongs to $P_0(\mathbb{R}^2)$.

We may now define the approximation spaces.

For the approximation of the hydraulic heads, we use piecewise constant functions:

$$S_h = \left\{ s_h = (s_h^{\text{m}}, s_h^{\text{f}}) \in L^2(\Omega^{\text{m}}) \times L^2(\Omega^{\text{f}}) \mid s_h^{\text{m}}|_T \in P_0(\mathbb{R}^3), \forall T \in \mathcal{T}_h; s_h^{\text{f}}|_K \in P_0(\mathbb{R}^2), \forall K \in \mathcal{K}_h \right\}. \quad (3.1)$$

For the velocity space, we first define two auxiliary spaces

$$V_h^{\text{m}} = \left\{ v_h^{\text{m}} \in L^2(\Omega^{\text{m}}) \mid \begin{array}{l} \forall T \in \mathcal{T}_h, v_h^{\text{m}}|_T \in \text{RTN}^0(T); \forall \alpha \in \mathcal{C}, v_h^{\text{m}}|_{\omega_\alpha} \in H(\text{div}, \omega_\alpha); \\ \forall F \in \mathcal{F}_h^{\text{N}}, \gamma_n^{\text{m},F} v_h^{\text{m}} = 0 \end{array} \right\},$$

and

$$V_h^{\text{f}} = \left\{ v_h^{\text{f}} \in L^2(\Omega^{\text{f}}) \mid \begin{array}{l} \forall K \in \mathcal{K}_h, v_h^{\text{f}}|_K \in \text{RT}^0(K); \forall \ell \in \mathcal{I}, v_h^{\text{f}}|_{\gamma_\ell} \in H(\text{div}, \gamma_\ell); \\ \forall \sigma \in \Sigma, \sum_{E \subset \sigma} \sum_{K, E \subset \partial K} \gamma_n^{\text{f},E} v_h^{\text{f}}|_K = 0, \text{ and } \forall E \in \mathcal{E}_h^{\text{N}}, \gamma_n^{\text{f},E} v_h^{\text{f}} = 0 \end{array} \right\}.$$

In the discrete setting we work with functions that are smoother than in the continuous setting. This means that we no longer have to deal with the rather cumbersome condition $(\text{div}(v_h^{\text{f}}) - [[u^\text{m}]]|_{\gamma_\ell}) \in L^2(\gamma_\ell), \ell \in \mathcal{I}$ (see Section 2.2 for details). Instead, we can simply consider functions for which the divergence and the normal trace are both in $L^2$. We can study the two spaces $V_h^{\text{m}}$ and $V_h^{\text{f}}$ separately and simply set:

$$V_h = V_h^{\text{m}} \times V_h^{\text{f}}. \quad (3.2)$$

Also, because the discrete functions are more regular than their continuous counterpart, we can impose the boundary conditions, and the conditions on the intersections, in a strong way. As a consequence, we have $V_h \subset V$.

We note the following alternative characterization of $V_h^{\text{m}}$. A function $v_h^{\text{m}} \in L^2(\Omega^{\text{m}})$ belongs to $V_h^{\text{m}}$ if and only if the following conditions hold:

$$v_h^{\text{m}}|_T \in \text{RTN}^0(T), \quad \forall T \in \mathcal{T}_h, \quad (3.3a)$$

$$[[v_h^{\text{m}}]]|_F = 0, \quad \forall F \in \mathcal{F}_h^{\text{in}}, \quad (3.3b)$$

$$\gamma_n^{\text{m},F} v_h^{\text{m}} = 0, \quad \forall F \in \mathcal{F}_h^{\text{N}}. \quad (3.3c)$$
### 3.1.2 The mixed approximation

We can now define the approximate problem for (2.6). We look for a couple $(u_h, p_h) \in V_h \times S_h$ satisfying:

$$ 
\begin{aligned}
a(u_h, v_h) - b(v_h, p_h) &= L_a(v_h), \quad \forall v_h \in V_h, \\
b(u_h, s_h) &= L_b(s_h), \quad \forall s_h \in S_h.
\end{aligned} 
\quad (3.4) $$

In the rest of this section, we prove that this discrete problem is well posed, and give an error estimate.

We first note that we are in the favorable situation, where the discrete kernel of the bilinear form $b$, defined in (2.4b), is a subspace of the continuous kernel characterized in (2.8). Let $W_h$ be the discrete kernel of the bilinear form $b$ over $V_h$:

$$ W_h = \{v_h \in V_h \mid b(v_h, s_h) = 0, \forall s_h \in S_h\}. $$

**Lemma 3.1.** The discrete kernel $W_h$ is a subspace of $W$.

*Proof.* Let $v_h \in W_h$. Then by definition:

$$ \sum_{T \in \mathcal{T}_h} \int_T (\nabla \cdot v_h^m) s_h^m + \sum_{K \in \mathcal{K}_h} \int_K (\nabla_\tau \cdot v_h^\mathfrak{f} - [[v_h^m]]_K) s_h^\mathfrak{f} = 0, \quad \forall (s_h^m, s_h^\mathfrak{f}) \in S_h. $$

This implies that, on each tetrahedral element $T \in \mathcal{T}_h$, $\nabla \cdot (v_h^m)|_T = 0$, so that $\nabla \cdot v_h^m = 0$ almost everywhere in $\omega_\alpha$, for all $\alpha \in \mathcal{C}$. Similarly, one sees that $\nabla_\tau \cdot v_h^\mathfrak{f} - [[v_h^m]]|_{\gamma_\ell} = 0$ almost everywhere in $\gamma_\ell$, for all $\ell \in \mathcal{I}$. Together, these two facts imply that $b(v_h, s) = 0, \forall s \in S$. ■

We now prove that the discrete problem is well-posed.

**Proposition 3.2.** The discrete problem (3.4) has a unique solution $(u_h, p_h) \in V_h \times S_h$.

*Proof.* Because the discrete problem (3.4) is equivalent to a square linear system, all we need to check is that the only solution for the homogeneous problem (i.e., with $L_a(v_h) = 0$ and $L_b(s_h) = 0$) is the trivial solution.

We start by taking $v_h = u_h$ and $s_h = p_h$, then adding the two equations in (3.4), to obtain:

$$ a(u_h, u_h) = 0. $$

The previous lemma, together with the coercivity of $a$ over $W$ (Lemma 2.6), implies that $u_h = 0$. Then we are left with $b(v_h, p_h) = 0, \forall v_h \in V_h$. This boils down to showing that if:

$$ \sum_{T \in \mathcal{T}_h} \int_T p_h^m \nabla \cdot v_h^m + \sum_{K \in \mathcal{K}_h} \int_K p_h^\mathfrak{f} (\nabla \cdot v_h^\mathfrak{f} - [[v_h^m]]_K) = 0, \quad \forall (v_h^m, v_h^\mathfrak{f}) \in V_h^m \times V_h^\mathfrak{f}, $$

then $p_h = 0$.

We take $v_h^m = 0$. Because the discrete divergence is surjective from $RT_0(K)$ onto $P_0(K)$, we see that $p_h^\mathfrak{f}|_K = 0, \forall K \in \mathcal{K}_h$. The same argument over the tetrahedral elements shows that $p_h^m|_T = 0, \forall T \in \mathcal{T}_h$. ■

We now proceed to check the discrete inf-sup condition. This will be done by constructing a Fortin operator (see [32, Lemma 26.9] or [12, Sec. 5.4.3]) that is related to the interpolation operator for the discrete spaces. We first define the operator at the element level, and then glue together the local definitions to obtain the global operator. This second step requires more smoothness on the solution than what is available from the existence theorem.

The definition of the local interpolation operator requires a more precise description of the degrees of freedom and the associated basis functions for the space $V_h$.
For the rock matrix domain, we can simply use the usual basis functions for the Raviart-Thomas-Nédélec space, with the only difference that each fracture face will carry two basis functions. Actually, fracture faces are boundary faces for the subdomains $\omega_\alpha$, $\alpha \in \mathcal{C}$. This means that $\dim(V_h^m) = \mathcal{N}(\mathcal{E}_h^{\text{in}}) + 2\mathcal{N}(\mathcal{F}_h^{\text{frac}}) + \mathcal{N}(\mathcal{F}_h^D)$.

For the fracture network, we will rely on the work by Vohralík et al. [66], where it is shown how to define the degrees of freedom and the basis functions so as to take into account the condition at the fracture intersections. The main point is that each intersection removes one basis function on the corresponding edges of the fracture mesh. We refer the reader to [66] for the details. As a result, we have $\dim(V_h^f) = \mathcal{N}(\mathcal{E}_h^{\text{inter}}) + \mathcal{N}(\mathcal{E}_h^D)$.

We begin by recalling the known interpolation operators for the local Raviart-Thomas(-Nédélec) spaces (see [12, 38]). As usual, the interpolation operator can only be defined for functions that are smoother than those in $V$. Following [12, Sec 2.5.1], we define the spaces (for a given $\mathfrak{s} > 2$):

$$ 
\begin{aligned}
W^{\mathfrak{m}}(T) &= \{u^{\mathfrak{m}} \in (L^{\mathfrak{s}}(T))^3 \mid \nabla \cdot u^{\mathfrak{m}} \in L^2(T)\}, \quad \forall T \in \mathcal{T}_h \\
W^{\mathfrak{f}}(K) &= \{u^{\mathfrak{f}} \in (L^{\mathfrak{s}}(K))^2 \mid \nabla_{\tau} \cdot u^{\mathfrak{f}} \in L^2(K)\}, \quad \forall K \in \mathcal{K}_h.
\end{aligned}
$$

We denote by $\pi_T^{\mathfrak{m}}: W^{\mathfrak{m}}(T) \to RTN^0(T)$ (resp. $\pi_K^{\mathfrak{f}}: W^{\mathfrak{f}}(K) \to RT^0(K)$) the local interpolation operators and we recall that they satisfy the following properties.

* In the rock matrix, for all tetrahedral element $T \in \mathcal{T}_h$:

$$ \int_F \gamma_n^{\mathfrak{m},F} \pi_T^{\mathfrak{m}} u^{\mathfrak{m}} = \int_F \gamma_n^{\mathfrak{m},F} u^{\mathfrak{m}}, \quad \forall F \in \mathcal{F}_h, F \subset \partial T \quad \text{and} \quad \int_T \nabla \cdot (\pi_T^{\mathfrak{m}} u^{\mathfrak{m}}) = \int_T \nabla \cdot u^{\mathfrak{m}}. \quad (3.5a) $$

* In the fracture network, for all triangles $K \in \mathcal{K}_h$:

$$ \int_E \gamma_n^{\mathfrak{f},E} \pi_K^{\mathfrak{f}} u^{\mathfrak{f}} = \int_E \gamma_n^{\mathfrak{f},E} u^{\mathfrak{f}}, \quad \forall E \in \mathcal{E}_h, E \subset \partial K, \quad \text{and} \quad \int_K \nabla_{\tau} \cdot (\pi_K^{\mathfrak{f}} u^{\mathfrak{f}}) = \int_K \nabla_{\tau} \cdot u^{\mathfrak{f}}. \quad (3.5b) $$

The next step is to build an interpolation operator for the spaces $V$ and $V_h$.

For the remaining of this section, we assume that the solution belongs to a space with more regularity that what is just afforded by $V$, denoted by $Z$ and defined as follows:

$$ Z = Z^{\mathfrak{m}} \times Z^{\mathfrak{f}}, $$

with

$$ 
\begin{aligned}
Z^{\mathfrak{m}} &= \{u^{\mathfrak{m}} \in L^2(\Omega^{\mathfrak{m}})^3 \mid u_{\alpha}^{\mathfrak{m}} \in H^{\mathfrak{s}}(\omega_{\alpha}), \forall \alpha \in \mathcal{C}\} \cap V^{\mathfrak{m}}, \\
Z^{\mathfrak{f}} &= \{u^{\mathfrak{f}} \in L^2(\Omega^{\mathfrak{f}})^2 \mid u_{\ell}^{\mathfrak{f}} \in H^{\mathfrak{s}}(\gamma_{\ell}), \forall \ell \in \mathcal{I}\} \cap V^{\mathfrak{f}},
\end{aligned}
$$

for some $\mathfrak{s} > 0$.

We also assume, for the remaining of this section, that the velocity $u$ belongs to the space $Z$. This is at least true when all fractures end on the boundaries of the domain, or on other fractures, as was assumed in [3, 41]. It is not clear whether or not the additional regularity still holds when there are immersed fractures.

The interpolation operator will be constructed separately for the rock matrix and for the fracture network, in each case by using already existing operators.

For the rock matrix domain, we can simply choose the usual interpolation operator for the global Raviart-Thomas-Nédélec space:

$$ \pi_h^{\mathfrak{m}}: Z^{\mathfrak{m}} \to V_h^{\mathfrak{m}}, \quad \pi_h^{\mathfrak{m}}(u^{\mathfrak{m}})|_T = \pi_T^{\mathfrak{m}}(u|_T^{\mathfrak{m}}). \quad (3.6) $$

Let $T \in \mathcal{T}_h$. By introducing the $L^2(T)$ orthogonal projection onto $P_0(\mathbb{R}^3)$, denoted by $P_T^{\mathfrak{m}}$, and its global counterpart from $L^2(\Omega^{\mathfrak{m}})$, denoted by $P_h^{\mathfrak{m}}$, we can summarize the main property of $\pi_h^{\mathfrak{m}}$ in
the following commutative diagram:

$$ \begin{array}{ccc} Z^{\mathfrak{m}} & \xrightarrow{\text{div}} & L^2(\Omega^{\mathfrak{m}}) \\ \pi_h^{\mathfrak{m}} & \downarrow & \downarrow \\ V_h^{\mathfrak{m}} & \xrightarrow{\text{div}} & P_0(\mathbb{R}^2) \end{array} \qquad (3.7) $$

For the fracture network part, we need to be more careful, as already mentioned. The difficulty is to take into account the condition of zero total flux at fracture-fracture intersections, where there can be more than two triangles meeting. We will use the interpolation operator proposed in [66], which will be denoted by $\pi_h^f$. We omit the details, and refer the reader to the cited paper, where it is shown how to eliminate one degree of freedom for each intersection, and how to subsequently combine the basis functions in order to obtain a basis for the space $V_h^f$. The end result is that the interpolation operator $\pi_h^f$ satisfies the same properties as $\pi_h^m$, in particular a commutative diagram similar to (3.7). Let $K \in \mathcal{K}_h$. We will also need the $L^2(K)$ orthogonal projection onto $P_0(\mathbb{R}^2)$, denoted by $P_K^f$, and its global counterpart from $L^2(\Omega^f)$, denoted by $P_h^f$.

Finally, we define the operator $\pi_h : Z \to V_h$ by:

$$ \pi_h(u) = (\pi_h^m(u^m), \pi_h^f(u^f)), \quad \forall u \in Z. $$

**Lemma 3.3.** The operator $\pi_h$ is a Fortin operator, i.e. it satisfies:

$$ 1. b(v - \pi_h v, s_h) = 0, \forall v \in Z, \forall s_h \in S_h; $$

$$ 2. \text{ There exists a constant } C > 0 \text{ such that } \|\pi_h v\|_V \le C \|v\|_V. $$

*Proof.* 1. We split the integrals defining $b$ into a sum over the elements in the matrix and the fracture meshes as:

$$ b(v - \pi_h v, s_h) = \sum_{T \in \mathcal{T}_h} I_T(v - \pi_h v, s_h) + \sum_{K \in \mathcal{K}_h} I_K(v - \pi_h v, s_h), \quad \forall v \in Z, \forall s_h \in S_h, $$

with

$$ I_T(v - \pi_h v, s_h) = \int_T s_h^{\mathfrak{m}} \nabla \cdot (v^{\mathfrak{m}} - \pi_h^{\mathfrak{m}} v^{\mathfrak{m}}), $$

and

$$ I_K(v - \pi_h v, s_h) = \int_K s_h^f \left[ \nabla_{\tau} \cdot (v^f - \pi_h^f v^f) - [[v^{\mathfrak{m}} - \pi_h^{\mathfrak{m}} v^{\mathfrak{m}}]]|_K \right]. $$

The result for $I_T$ is clear from the commutative diagram property (3.7). For $I_K$, because $v$ is in the more regular space $Z$, we may split the integral in two:

$$ I_K(v - \pi_h v, s_h) = \int_K (s_h^f)|_K \nabla_{\tau} \cdot (v^f|_K - \pi_K^f(v^f|_K)) - \int_K (s_h^f)|_K [[(v^{\mathfrak{m}} - \pi_h^{\mathfrak{m}} v^{\mathfrak{m}})]|_K. $$

The first integral is zeros again by the commutative diagram property of the interpolation operator $\pi_h^f$. The second one is also zero thanks to the definition of the degrees of freedom for $V_h^f$ (see (3.5a)).

2. This is a consequence of the continuity of the individual interpolation operators. We must prove that there exists $C > 0$ such that:

$$ \begin{aligned} & \sum_{\alpha \in C} \int_{\omega_{\alpha}} |\pi_h^m v_{\alpha}^m|^2 + |\nabla \cdot \pi_h^m v_{\alpha}^m|^2 + \sum_{\ell \in I} \int_{\gamma_{\ell}} |\pi_h^f v_{\ell}^f|^2 + |\nabla_{\tau} \cdot \pi_h^f v_{\ell}^f - [[\pi_h^m v^m]]|_{\gamma_{\ell}}|^2 \\ & \le C \left[ \sum_{\alpha \in C} \int_{\omega_{\alpha}} |v_{\alpha}^m|^2 + |\nabla \cdot v_{\alpha}^m|^2 + \sum_{\ell \in I} \int_{\gamma_{\ell}} |v_{\ell}^f|^2 + |\nabla_{\tau} \cdot v_{\ell}^f - [[v^m]]|_{\gamma_{\ell}}|^2 \right]. \end{aligned} $$
The known continuity properties of the Fortin operators $\pi_h^m$ and $\pi_h^f$ (see for instance [38, 12])
mean that it is sufficient to check that there holds:

$$
\int_K |\nabla_\tau \cdot \pi_K^f(v_\ell^f) - [[\pi_h^m v^m]]|_K|^2 \le C \int_K |\nabla_\tau \cdot v_\ell^f - [[v^m]]|_K|^2, \quad \forall K \in \mathcal{K}_h. \quad (3.8)
$$

Recall that for each triangular element $K \in \mathcal{K}_h$ there exists two tetrahe-
dra $(T^+(K), T^-(K)) \in \mathcal{T}_h^2$ that share $K$ as a common face. We rewrite (3.5a) as:

$$
\gamma_{n}^{m,K} \pi_{T^{\pm}(K)} v^{m} = P_{K}^{\mathcal{F}} (\gamma_{n}^{m,K} v_{|T^{\pm}(K)}^{m}),
$$

so that

$$
[[\pi_h^m v^m]]|_K = P_K^f ([[v^m]]|_K),
$$

We also have:

$$
\nabla_{\tau} \cdot \pi_{K}^{\mathcal{F}}(v_{\ell}^{\mathcal{F}}) = P_{K}^{\mathcal{F}}(\nabla_{\tau} \cdot v_{\ell}^{\mathcal{F}}).
$$

Inequality (3.8) is thus a consequence of

$$
\int_K |P_K^\mathcal{F} (\nabla_\tau \cdot v_\ell^\mathcal{F} - [[v^\mathcal{m}]]|_K)|^2 \le \int_K |\nabla_\tau \cdot v_\ell^\mathcal{F} - [[v^\mathcal{m}]]|_K|^2,
$$

which is true because $P_K^\mathcal{F}$ is an orthogonal projector.

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

The proof of the next result now follows from the continuous inf-sup condition, Lemma 2.7,
and [12, Proposition 5.4.3].

**Corollary 3.4.** *The bilinear form b satisfies the discrete uniform inf-sup condition:*

$$
\exists \tilde{\beta} > 0, \forall h > 0, \quad \inf_{s_h \in S_h, s_h \neq 0} \sup_{v_h \in V_h, v_h \neq 0} \frac{b(v_h, s_h)}{\|v_h\|_V \|s_h\|_S} \ge \tilde{\beta} > 0. \quad (3.9)
$$

We now state the main result of this section.

**Theorem 3.5.** Let $(u,p) \in V \times S$ be the solution of the continuous problem (2.6) and $(u_h, p_h) \in V_h \times S_h$ be the solution of the discrete problem (3.4). We recall that the meshes $\mathcal{T}_h$ and $\mathcal{K}_h$ are both assumed to satisfy the mesh regularity condition. Finally, assume that $u \in Z$.

1. The following estimates hold:

$$
\|u - u_h\|_V \le C \inf_{v_h \in V_h} \|u - v_h\|_V, \qquad (3.10a)
$$

$$
\|p - p_h\|_S \le C \left( \inf_{s_h \in S_h} \|p - s_h\|_S + \inf_{v_h \in V_h} \|u - v_h\|_V \right). \quad (3.10\text{b})
$$

2. If we further assume that the solution has additional regularity, namely:

$$
\begin{array}{l}
u_{\alpha}^{\mathcal{m}} \in H^1(\omega_{\alpha})^3, \text{ and } \nabla \cdot u_{\alpha}^{\mathcal{m}} \in H^1(\omega_{\alpha}), \forall \alpha \in \mathcal{C} \\
u_{\ell}^{\mathcal{f}} \in H^1(\gamma_{\ell})^2 \text{ and } \nabla_{\tau} \cdot u_{\ell}^{\mathcal{f}} \in H^1(\gamma_{\ell}), \forall \ell \in \mathcal{I}, \\
p \in H_D,
\end{array}
$$

then the following error estimates hold:

$$
\|u - u_h\|_V \le h \|u\|_1, \qquad (3.11a)
$$

$$
\|p - p_h\|_S \le h (\|u\|_1 + \|p\|_{H^1(\Omega)}), \qquad (3.11\text{b})
$$

where

$$
\|u\|_1^2 = \sum_{\alpha \in \mathcal{C}} \left( |u_{\alpha}^{\mathcal{m}}|_{H^1(\omega_{\alpha})}^2 + |u_{\alpha}^{\mathcal{m}}|_{H^1(\omega_{\alpha})}^2 \right) + \sum_{\ell \in \mathcal{I}} \left( |u_{\ell}^{\mathcal{f}}|_{H^1(\gamma_{\ell})}^2 + |u_{\ell}^{\mathcal{f}}|_{H^1(\gamma_{\ell})}^2 \right).
$$
Proof. The proof of point 1 follows from the general theory of mixed approximation [12, 38, 32]. We note that the estimates here simplify because we are in the favorable situation where $W_h \subset W$.

Here, we detail the proof of point 2.

Since $u \in Z$, we have:

$$
\begin{align*}
\|u - \pi_h u\|_V^2 &= \sum_{\alpha \in C} \sum_{T \in \mathcal{T}_h^\alpha} \int_T |u_\alpha^\mathfrak{m} - \pi_h^\mathfrak{m} u_\alpha^\mathfrak{m}|^2 + |\nabla \cdot u_\alpha^\mathfrak{m} - \nabla \cdot (\pi_h^\mathfrak{m} u_\alpha^\mathfrak{m})|^2 \\
&\quad + \sum_{\ell \in \mathcal{I}} \sum_{K \in \mathcal{K}_h^\ell} \left[ \int_K |u_\ell^\mathfrak{f} - \pi_h^\mathfrak{f} u_\ell^\mathfrak{f}|^2 + \int_K |\nabla_\tau \cdot u_\ell^\mathfrak{f} - \nabla_\tau \cdot (\pi_h^\mathfrak{f} u_\ell^\mathfrak{f}) - [[u_\alpha^\mathfrak{m} - \pi_h^\mathfrak{m} u_\alpha^\mathfrak{m}]]|_K|^2 \right] \\
&\le \sum_{\alpha \in C} \sum_{T \in \mathcal{T}_h^\alpha} \underbrace{\int_T |u_\alpha^\mathfrak{m} - \pi_h^\mathfrak{m} u_\alpha^\mathfrak{m}|^2 + |\nabla \cdot u_\alpha^\mathfrak{m} - \nabla \cdot (\pi_h^\mathfrak{m} u_\alpha^\mathfrak{m})|^2}_{=: I_{T,\alpha}} \\
&\quad + 2 \sum_{\ell \in \mathcal{I}} \sum_{K \in \mathcal{K}_h^\ell} \left[ \underbrace{\int_K |u_\ell^\mathfrak{f} - \pi_h^\mathfrak{f} u_\ell^\mathfrak{f}|^2 + |\nabla_\tau \cdot u_\ell^\mathfrak{f} - \nabla_\tau \cdot (\pi_h^\mathfrak{f} u_\ell^\mathfrak{f})|^2}_{=: I_{K,\ell,1}} + \underbrace{\int_K |[[u^\mathfrak{m} - \pi_h^\mathfrak{m} u^\mathfrak{m}]]|_K|^2}_{=: I_{K,\ell,2}} \right] \\
&= \sum_{\alpha \in C} \sum_{T \in \mathcal{T}_h^\alpha} I_{T,\alpha} + 2 \sum_{\ell \in \mathcal{I}} \sum_{K \in \mathcal{K}_h^\ell} (I_{K,\ell,1} + I_{K,\ell,2}).
\end{align*}
$$

The estimates for $I_{T,\alpha}$ and $I_{K,\ell,1}$ are obtained separately for each subdomain $\omega_\alpha$ and each fracture $\gamma_\ell$ by using standard arguments (e.g. [38, Thm. 3.6], used here in the lowest order case $k=0$). Then there exists a constant $C > 0$ independent of $h^\mathfrak{m}$ and $h^\mathfrak{f}$ (but possibly depending on the mesh regularity parameters $c^\mathfrak{m}$ and $c^\mathfrak{f}$) such that:

$$ \sum_{\alpha \in C} \sum_{T \in \mathcal{T}_h^\alpha} I_{T,\alpha} \le C(h^\mathfrak{m})^2 \sum_{\alpha \in C} \left( |u_{\alpha}^{\mathfrak{m}}|_{H^1(\omega_\alpha)}^2 + |\nabla \cdot u_{\alpha}^{\mathfrak{m}}|_{H^1(\omega_\alpha)}^2 \right), \quad (3.12) $$

$$ \sum_{\ell \in \mathcal{I}} \sum_{K \in \mathcal{K}_h^\ell} I_{K, \ell, 1} \le C(h^\mathfrak{f})^2 \sum_{\ell \in \mathcal{I}} \left( |u_\ell^\mathfrak{f}|_{H^1(\gamma_\ell)}^2 + |\nabla_\tau \cdot u_\ell^\mathfrak{f}|_{H^1(\gamma_\ell)}^2 \right). \quad (3.13) $$

In order to obtain a bound for $I_{K,\ell,2}$, we first recall that each triangular element $K \in \mathcal{K}_h$ is the boundary of two tetrahedra $T^+(K)$ and $T^-(K)$ belonging to two different subdomains $\omega_{\alpha\pm}$.
Thus:

$$ I_{K,\ell,2} \le 2 \int_K |\gamma_n^\mathfrak{m}(u_{\alpha^+}^\mathfrak{m} - \pi_h^\mathfrak{m} u_{\alpha^+}^\mathfrak{m})|^2 + 2 \int_K |\gamma_n^\mathfrak{m}(u_{\alpha^-}^\mathfrak{m} - \pi_h^\mathfrak{m} u_{\alpha^-}^\mathfrak{m})|^2. $$

We treat only the first integral and we omit the + sign. By using first the definition of the interpolation operator $\pi_h^\mathfrak{m}$ and then [38, Lemma 3.18] to bound the interpolation error for the normal trace, there holds:

$$ \int_K |\gamma_n^\mathfrak{m}(u_\alpha^\mathfrak{m} - \pi_h^\mathfrak{m} u_\alpha^\mathfrak{m})|^2 = \int_K |\gamma_n^\mathfrak{m} u_\alpha^\mathfrak{m} - \gamma_n^\mathfrak{m} \pi_T^\mathfrak{m}(u_\alpha^\mathfrak{m})|_T^2 \le C|K| |u_\alpha^\mathfrak{m}|_{H^1(\omega_\alpha)}^2 \le C h_K^2 |u_\alpha^\mathfrak{m}|_{H^1(\omega_\alpha)}^2. $$

By summing over all triangles and all fractures, we obtain:

$$ \sum_{\ell \in \mathcal{I}} \sum_{K \in \mathcal{K}_h^\ell} I_{K, \ell, 2} \le C(h^\mathfrak{f})^2 \sum_{\alpha \in C} |u_\alpha^\mathfrak{m}|_{H^1(\omega_\alpha)}^2. \quad (3.14) $$

Estimate (3.11a) results from summing equations (3.12), (3.13) and (3.14). Estimate (3.11b) follows in the same way from standard approximation results for scalar functions (see for instance [22]) and from the one for the velocity.
### 3.1.3 Mixed-hybrid formulation

In this section and the next we will use the following notation to relate tetrahedral elements to their faces as well as triangles to their edges:

* for a given tetrahedron $T \in \mathcal{T}_h$, we let $\mathcal{F}(T) = \{F \in \mathcal{F}_h | F \text{ is a face of } T\}$;

* for a given triangle $K \in \mathcal{K}_h$, we let $\mathcal{E}(K) = \{E \in \mathcal{E}_h | E \text{ is an edge of } K\}$;

Furthermore we let $\mathcal{F}(T)^{\text{frac}} = \mathcal{F}(T) \cap \mathcal{F}_h^{\text{frac}}$ and $\mathcal{F}(T)^{\text{in,N}} = \mathcal{F}(T) \cap \mathcal{F}_h^{\text{in,N}}$.

In order to obtain the mixed-hybrid formulation, we relax the continuity of the normal trace of the velocity across the interior faces. We need the following discontinuous spaces:

$$V_h^{\text{hyb,m}} = \{v_h^m \in L^2(\Omega^m) \mid v_h^m|_T \in \text{RTN}^0(T), \forall T \in \mathcal{T}_h\}, \quad (3.15a)$$

$$V_h^{\text{hyb,f}} = \{v_h^f \in L^2(\Omega^f) \mid v_h^f|_K \in \text{RT}^0(K), \forall K \in \mathcal{K}_h\}, \quad (3.15b)$$

and finally

$$V_h^{\text{hyb}} = V_h^{\text{hyb,m}} \times V_h^{\text{hyb,f}}. \quad (3.16)$$

We note the following simple but crucial characterizations of functions in $V_h^{\text{hyb}}$. For $v_h = (v_h^m, v_h^f) \in V_h^{\text{hyb}}$, we have the following equivalences:

$$v_h^m \in V_h^m \iff [[v_h^m]]|_F = 0, \forall F \in \mathcal{F}_h^{\text{in}}, \text{ and } \gamma_n^{m,F} v_h^m = 0, \forall F \in \mathcal{F}_h^N, \quad (3.17a)$$

$$v_h^f \in V_h^f \iff \left[ \sum_{K \in \mathcal{K}(E)} \gamma_n^{f,E} v_h^f|_K = 0, \forall E \in \mathcal{E}_h^{\text{in}} \cup \mathcal{E}_h^{\text{inter}}, \text{ and } \gamma_n^{f,E} v_h^f = 0, \forall E \in \mathcal{E}_h^D \right]. \quad (3.17b)$$

We note that the characterization (3.17a) of $V_h^m$ excludes fracture faces. We also note that the characterization (3.17b) of $V_h^f$ actually takes into account all edges in the mesh: if an edge is not an intersection, the sum is over two triangles and expresses the usual flux continuity. The proof of the equivalences is well known, and follows exactly as in [12].

We formulate the characterization in (3.17) by means of Lagrange multipliers that are piecewise constant functions on the interior faces and all the edges. We define the following spaces:

$$M_{h,0,D}^{m} = \{ \mu_h^m \in L^2(\Omega^m) \mid \mu_h^m|_F \in P_0(F), \forall F \in \mathcal{F}_h \setminus \mathcal{F}_h^{\text{frac}}, \text{ and } \mu_h^m|_F = 0, \forall F \in \mathcal{F}_h^D \}, \quad (3.18a)$$

$$M_{h,0,D}^{f} = \{ \mu_h^f \in L^2(\Omega^f) \mid \mu_h^f|_E \in P_0(E), \forall E \in \mathcal{E}_h, \text{ and } \mu_h^f|_E = 0, \forall E \in \mathcal{E}_h^D \}, \quad (3.18b)$$

and finally:

$$M_{h,0,D} = M_{h,0,D}^{m} \times M_{h,0,D}^{f}. \quad (3.19)$$

We then define the following bilinear form over $V_h^{\text{hyb}} \times M_{h,0,D}$:

$$c(v_h, \mu_h) = \sum_{F \in \mathcal{F}_h^{\text{in}}} \int_F [[v_h^m]]|_F \mu_h^m|_F + \sum_{F \in \mathcal{F}_h^N} \int_F \gamma_n^{m,F} v_h^m \mu_h^m|_F + \sum_{E \in \mathcal{E}_h} \int_E \sum_{K \in \mathcal{K}(E)} \gamma_n^{f,E} v_h^f|_K \mu_h^f|_E. \quad (3.20)$$

Then, the conditions in (3.17) can be equivalently rewritten as:

$$\text{for } v_h \in V_h^{\text{hyb}}, \quad (v_h \in V_h \iff c(v_h, \mu_h) = 0, \forall \mu_h \in M_{h,0,D}).$$

Because functions in $V_h^{\text{hyb}}$ are not necessarily in $V_h$, we need to introduce a mesh dependent bilinear form for all $u_h \in V_h^{\text{hyb}}$ and for all $s_h \in S_h$:

$$b_T(u_h, s_h) = \sum_{T \in \mathcal{T}_h} \int_T \nabla \cdot u_h^m|_T s_h^m|_T + \sum_{K \in \mathcal{K}_h} \int_K (\nabla_{\tau} \cdot u_h^f|_K - [[u_h^m]]|_K) s_h^f|_K. \quad (3.21)$$
We can then write the hybridized formulation of the discrete mixed problem (3.4). We look for a triple $(u_h, p_h, \lambda_h) \in V_h^{\text{hyb}} \times S_h \times M_{h,0,D}$ such that:

$$
\begin{align}
a(u_h, v_h) - b_T(v_h, p_h) + c(\lambda_h, v_h) &= L_a(v_h), && \forall v_h \in V_h^{\text{hyb}}, \nonumber \\
b_T(u_h, s_h) &= L_b(s_h), && \forall s_h \in S_h, \tag{3.22} \\
c(u_h, \mu_h) &= 0, && \forall \mu_h \in M_{h,0,D}. \nonumber
\end{align}
$$

We have given the same name to the vector unknown in both problems (3.4) and (3.22) because of the well-known equivalence theorem (see [12, Thm 7.2.1] or [32, Prop. 52.1]) that states that the solution $u_h$ of (3.22) is actually in $V_h$ and that $(u_h, p_h)$ is the solution of (3.4).

As is well known (see the references above), the main advantage of the mixed-hybrid finite element formulation (3.22) with respect to the mixed finite element formulation (3.4) is that the former leads to a linear system that is both of smaller dimension and with a symmetric and positive definite matrix. In the context of this work, the choice of the mixed-hybrid formulation has the further advantage that the hydraulic head degrees of freedom in the rock matrix elements are located “at the same place” as the Lagrange multipliers (trace of hydraulic head) for the rock matrix element. At the element level, this means that we can work “as if” there were Lagrange multipliers on all element faces, and that this additional degree of freedom coincides with the hydraulic head on the fracture face. This is the point of view taken in [44].

## 3.2 Mixed-hybrid finite element algebraic formulation

In this section we describe the linear system that is obtained after basis functions and degrees of freedom are chosen for the spaces involved in (3.22).

The global vectors, degrees of freedom and associated functions are given in Table 1. The discrete spaces for the fluxes $V_h^{\text{hyb},\mathfrak{m}}$ and $V_h^{\text{hyb},\mathfrak{f}}$, and thus $V_h^{\text{hyb}}$, are provided with the natural bases obtained from the product of the local Raviart-Thomas(-Nédélec) spaces. The Raviart-Thomas(-Nédélec) basis functions for triangles (resp. tetrahedra) are denoted by $w_{K,e}^{\mathfrak{f}}$ (resp. $w_{T,f}^{\mathfrak{m}}$), see [12, 38, 62] for the definitions. The spaces for the hydraulic head are piecewise constant functions over each tetrahedron or triangle respectively. Finally, the spaces for the Lagrange multipliers are piecewise constant on the inner faces of the tetrahedral mesh, and on the edges of the triangular mesh respectively.

<table><caption>Table 1: Definition of discrete functions, degrees of freedom, and associated global vectors.</caption><thead><tr><th>Discrete functions</th><th>Degrees of freedom</th><th>Global vector</th></tr></thead><tbody><tr><td>$u_h^{\mathfrak{m}} = \sum_{T \in \mathcal{T}_h} \sum_{f=1}^4 q_{T,f}^{\mathfrak{m}} w_{T,f}^{\mathfrak{m}}$<br>$u_h^{\mathfrak{f}} = \sum_{K \in \mathcal{K}_h} \sum_{e=1}^3 q_{K,e}^{\mathfrak{f}} w_{K,e}^{\mathfrak{f}}$</td><td>$Q^{\mathfrak{m}} = (q_{T,f}^{\mathfrak{m}})_{T \in \mathcal{T}_h, f=1,...,4} \in \mathbb{R}^{4\mathcal{N}(\mathcal{T}_h)}$<br>$Q^{\mathfrak{f}} = (q_{K,e}^{\mathfrak{f}})_{K \in \mathcal{K}_h, e=1,...,3} \in \mathbb{R}^{3\mathcal{N}(\mathcal{K}_h)}$</td><td>$Q = \begin{pmatrix} Q^{\mathfrak{m}} \\ Q^{\mathfrak{f}} \end{pmatrix}$</td></tr><tr><td>$p_h^{\mathfrak{m}} = \sum_{T \in \mathcal{T}_h} p_T^{\mathfrak{m}} \mathbb{1}_T$<br>$p_h^{\mathfrak{f}} = \sum_{K \in \mathcal{K}_h} p_K^{\mathfrak{f}} \mathbb{1}_K$</td><td>$P^{\mathfrak{m}} = (p_T^{\mathfrak{m}})_{T \in \mathcal{T}_h} \in \mathbb{R}^{\mathcal{N}(\mathcal{T}_h)}$<br>$P^{\mathfrak{f}} = (p_K^{\mathfrak{f}})_{K \in \mathcal{K}_h} \in \mathbb{R}^{\mathcal{N}(\mathcal{K}_h)}$</td><td>$P = \begin{pmatrix} P^{\mathfrak{m}} \\ P^{\mathfrak{f}} \end{pmatrix}$</td></tr><tr><td>$\lambda_h^{\mathfrak{m}} = \sum_{F \in \mathcal{F}_h^{\text{in},N}} \lambda_F^{\mathfrak{m}} \mathbb{1}_F$<br>$\lambda_h^{\mathfrak{f}} = \sum_{E \in \mathcal{E}_h} \lambda_E^{\mathfrak{f}} \mathbb{1}_E$</td><td>$\Lambda^{\mathfrak{m}} = (\lambda_F^{\mathfrak{m}})_{F \in \mathcal{F}_h^{\text{in},N}} \in \mathbb{R}^{\mathcal{N}(\mathcal{F}_h^{\text{in},N})}$<br>$\Lambda^{\mathfrak{f}} = (\lambda_E^{\mathfrak{f}})_{E \in \mathcal{E}_h \setminus \mathcal{E}_h^D} \in \mathbb{R}^{\mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)}$</td><td>$\Lambda = \begin{pmatrix} \Lambda^{\mathfrak{m}} \\ \Lambda^{\mathfrak{f}} \end{pmatrix}$</td></tr></tbody></table>

An illustration with two tetrahedral elements $T_1$ and $T_2$ (rock) and one triangle $K$ (fracture) is displayed on Figure 2:
We now define the right hand sides as $\mathbf{g} = \begin{pmatrix} g^m \\ g^\mathcal{F} \end{pmatrix}$ and $\mathbf{f} = \begin{pmatrix} f^m \\ f^\mathcal{F} \end{pmatrix}$.

The vector $f^m = (f_T^m)_{T \in \mathcal{T}_h}$ (resp. $f^\mathcal{F} = (f_K^\mathcal{F})_{K \in \mathcal{K}_h}$) is of size $\mathcal{N}(\mathcal{T}_h)$ (resp. $\mathcal{N}(\mathcal{K}_h)$), such that for each $T \in \mathcal{T}_h$ (resp. each $K \in \mathcal{K}_h$): $f_T^m = \int_T f^m$, (resp. $f_K^\mathcal{F} = \int_K f^\mathcal{F}$).

The vector $g^m$ (resp. $g^\mathcal{F}$) is of size $4\mathcal{N}(\mathcal{T}_h)$ (resp. $3\mathcal{N}(\mathcal{K}_h)$), such that

$$g_{T,f}^m = \begin{cases} -\int_F g^{m,D} & \text{if } F \subset \Gamma^D, \\ 0 & \text{otherwise,} \end{cases} \quad \text{where face } F \in \mathcal{F}(T) \text{ has local number } f \text{ in } T, f \in \{1, 2, 3, 4\}.$$

$$g_{K,e}^\mathcal{F} = \begin{cases} -\int_E g^{\mathcal{F},D} & \text{if } E \subset \Sigma^D, \\ 0 & \text{otherwise,} \end{cases} \quad \text{where edge } E \in \mathcal{E}(K) \text{ has local number } e \text{ in } K, e \in \{1, 2, 3\}.$$

When a hybrid formulation is used, the matrices $\mathcal{A}$ and also $\mathcal{B}$ are block diagonal. It is usual to take advantage of this feature to reduce the linear system (3.23) by locally eliminating first the fluxes and then the hydraulic heads to obtain a linear system that only involves the Lagrange multipliers [12, Sec. 7.2]. Furthermore, the matrix for this linear system can be shown to be symmetric and positive definite. After the Lagrange multipliers have been computed, the hydraulic head, and then the fluxes, can be recovered by solving local problems.

In our application, this elimination procedure must be adapted to take into account the coupling between the rock matrix and the fractures (this is due to the presence of the coupling matrix $B^{m,\mathcal{F}}$ so that $\mathcal{B}$ is no longer block diagonal). Because the process is quite technical, the details are presented in Appendix B. The end result is the linear system obtained in (B.8):

$$\underbrace{\begin{pmatrix} A^{m,m} & A^{m,\mathcal{F}} & 0 \\ (A^{m,\mathcal{F}})^T & A^{\mathcal{F},\mathcal{F}} + D^\mathcal{F} & -R^\mathcal{F} \\ 0 & -(R^\mathcal{F})^T & M^\mathcal{F} \end{pmatrix}}_{\mathbf{A}} \begin{pmatrix} \mathbf{\Lambda}^m \\ \mathbf{P}^\mathcal{F} \\ \mathbf{\Lambda}^\mathcal{F} \end{pmatrix} = \underbrace{\begin{pmatrix} \tilde{\mathcal{F}}^{m} + \mathcal{V}^m \\ \tilde{\mathcal{F}}^{m,\mathcal{F}} + V^{m,\mathcal{F}} + \mathcal{F}^\mathcal{F} \\ \mathcal{V}^\mathcal{F} \end{pmatrix}}_{\mathbf{F}},$$

with the block matrices and vectors defined in Appendix B.3.

This linear system will be denoted by:

$$A\Lambda = F, \qquad (3.25)$$

$$\text{with the unknown vector } \Lambda = \begin{pmatrix} \mathbf{\Lambda}^m \\ \mathbf{P}^\mathcal{F} \\ \mathbf{\Lambda}^\mathcal{F} \end{pmatrix}.$$

It can be observed that the matrix $\mathbf{A}$ can be assembled from a rock matrix “part”: $\mathbf{A}^m = \begin{pmatrix} A^{m,m} & A^{m,\mathcal{F}} \\ (A^{m,\mathcal{F}})^T & A^{\mathcal{F},\mathcal{F}} \end{pmatrix}$, and a fracture network “part”: $\mathbf{A}^\mathcal{F} = \begin{pmatrix} D^\mathcal{F} & -R^\mathcal{F} \\ -(R^\mathcal{F})^T & M^\mathcal{F} \end{pmatrix}$, provided that the triangular faces from the fracture network mesh $\mathcal{K}_h$ are identified with the non-inner triangular faces in the rock matrix mesh $\mathcal{T}_h$, and that the corresponding blocks are added. It is shown in Appendix B that the matrix $\mathbf{A}$ is symmetric and positive definite.

Once the system (3.25) has been solved, the mean hydraulic head $\mathbf{P}^m$ can be recovered by solving (B.3) for each tetrahedron, and the fluxes can be recovered by solving (B.1a) for each tetrahedron and (B.2a) for each triangle.

## 3.3 A discrete bilinear form

In this section, we follow [24, 23] to obtain a mesh-dependent bilinear form for the discrete problem. This will prove to be useful in [46] for showing that the discrete problem fits the GenEO
framework [63, 28]. The structure of the discrete problem, in particular the fact that the mean hydraulic heads on fracture network triangles are not eliminated, means that the details are slightly different from the computations in the above references.

Before we state the result, we follow [24] and introduce two auxiliary problems. First, for a given vector $\Lambda = (\Lambda^m, P^f, \Lambda^g)^T$, we define $(Q_\Lambda^m, P_\Lambda^m, Q_\Lambda^g)^T$ as the solution of:

$$A^m Q_\Lambda^m - B^m P_\Lambda^m = -B^{m,f} P^f - C^m \Lambda^m, \quad (3.26a)$$

$$(B^m)^T Q_\Lambda^m = 0, \quad (3.26b)$$

$$A^f Q_\Lambda^f = B^f P^f - C^f \Lambda^f. \quad (3.26c)$$

Second, for given right-hand sides $(g^m, f^m, g^f)^T$, we define $(\overline{Q}^m, \overline{P}^m, \overline{Q}^f)^T$ as a solution of the problem:

$$A^m \overline{Q}^m - B^m \overline{P}^m = g^m, \quad (3.27a)$$

$$(B^m)^T \overline{Q}^m = f^m, \quad (3.27b)$$

$$A^f \overline{Q}^f = g^f. \quad (3.27c)$$

We emphasize that both problems can be solved locally, that is, element by element. We then define:

$$u_h^m(\Lambda) = \sum_{T \in \mathcal{T}_h} (Q_\Lambda^m)_{T,f} w_{T,f}^m, \quad u_h^f(\Lambda) = \sum_{K \in \mathcal{K}_h} (Q_\Lambda^f)_{K,e} w_{K,e}^f, \quad p_h^m(\Lambda) = \sum_{T \in \mathcal{T}_h} (P_\Lambda^m)_{T \perp T}$$

where $(w_{T,f}^m)$ (resp. $w_{K,e}^f$) are the basis functions of $V_h^{\text{hyb},m}$ (resp. $V_h^{\text{hyb},f}$), for $T \in \mathcal{T}_h$, $f \in \{1, 2, 3, 4\}$ (resp. $K \in \mathcal{K}_h$, $e \in \{1, 2, 3\}$), as well as $u_h(\Lambda) = (u_h^m(\Lambda), u_h^f(\Lambda))$. We say that $u_h(\Lambda)$ and $p_h^m(\Lambda)$ are the solutions of problem (3.26) with data $\Lambda$.

**Proposition 3.6.** A vector $(\begin{pmatrix} Q \\ P \\ \Lambda \end{pmatrix})$ is the solution of the mixed-hybrid linear system (3.23) if and only if

$$1. \quad \begin{pmatrix} Q^m \\ P^m \\ Q^f \end{pmatrix} = \begin{pmatrix} Q_\Lambda^m \\ P_\Lambda^m \\ Q_\Lambda^f \end{pmatrix} + \begin{pmatrix} \overline{Q}^m \\ \overline{P}^m \\ \overline{Q}^f \end{pmatrix};$$

2. the vector $\Lambda$ is the solution of the discrete variational equation

$$a_T(\Lambda, \hat{\Lambda}) = L(\hat{\Lambda}), \quad \forall \hat{\Lambda}, \quad (3.28)$$

where $a_T(\Lambda, \hat{\Lambda})$ is a bilinear form given by an integral:

$$a_T(\Lambda, \hat{\Lambda}) = a(u_h(\Lambda), u_h(\hat{\Lambda})) = \int_{\Omega^m} (\mathcal{K}^m)^{-1} u_h^m(\Lambda) \cdot u_h^m(\hat{\Lambda}) + \int_{\Omega^f} (\mathcal{K}^f)^{-1} u_h^f(\Lambda) \cdot u_h^f(\hat{\Lambda}), \quad (3.29)$$

and $u_h(\Lambda)$ (resp. $u_h(\hat{\Lambda})$) is the solution of the auxiliary problem (3.26) with data vector $\Lambda$ (resp. $\hat{\Lambda}$). Similarly, the linear form $L(\hat{\Lambda})$ is given by:

$$L(\hat{\Lambda}) = \int_{\Omega^m} p_h^m(\hat{\Lambda}) f^m + \int_{\Gamma^D} \gamma_n^{m, \Gamma^D} u_h^m(\hat{\Lambda}) g^m + \int_{\Omega^f} p_h^f(\hat{\Lambda}) f^f + \int_{\Sigma^D} \gamma_n^{f, \Sigma^D} u_h^f(\hat{\Lambda}) g^f, \quad (3.30)$$

$p_h^m(\hat{\Lambda})$ is the solution of (3.26) with data vector $\hat{\Lambda}$, and $f$ and $g$ are the right hand side data for the mixed-hybrid problem (3.23).

**Remark 3.7.** The bilinear form (3.29) and the linear form (3.30) can be computed locally.
Proof. We will apply Theorem A.1 in [24], as it is written at the right level of generality for our purpose. To distinguish between the abstract quantities used in that reference and those in the present work, we will use a superscript CG for all quantities defined in [24]. In order to fit that paper’s framework, we need a different decomposition of the vectors and the matrices. Thus we let

$$ \boldsymbol{p}^{\text{CG}} = \begin{pmatrix} \boldsymbol{Q}^{\text{m}} \\ \boldsymbol{Q}^{\text{f}} \end{pmatrix}, \quad v^{\text{CG}} = \boldsymbol{P}^{\text{m}}, \quad \lambda^{\text{CG}} = \begin{pmatrix} \boldsymbol{P}^{\text{f}} \\ \boldsymbol{\Lambda}^{\text{m}} \\ \boldsymbol{\Lambda}^{\text{f}} \end{pmatrix}, \quad \boldsymbol{\alpha}^{\text{CG}} = \begin{pmatrix} g^{\text{m}} \\ g^{\text{f}} \end{pmatrix}, \quad \beta^{\text{CG}} = f^{\text{m}}, \quad \gamma^{\text{CG}} = \begin{pmatrix} f^{\text{f}} \\ 0 \\ 0 \end{pmatrix}, $$

and

$$ A^{\text{CG}} = \begin{pmatrix} A^{\text{m}} & 0 \\ 0 & A^{\text{f}} \end{pmatrix}, \quad B^{\text{CG}} = \begin{pmatrix} B^{\text{m}} \\ 0 \end{pmatrix}^T, \quad C^{\text{CG}} = \begin{pmatrix} -(B^{\text{m,f}})^T & -(C^{\text{m}})^T & 0 \\ (B^{\text{f}})^T & 0 & -(C^{\text{f}})^T \end{pmatrix}. $$

Note that we have had to adjust the transposes and the signs because our conventions are different than those in [24] (compare (3.23) above and (A.4) in [24]). Also the matrix $D^{\text{CG}}$ is zero.

The hypotheses of Theorem A.1 are clearly satisfied ($A^{\text{CG}}$ is symmetric and positive definite, and $B^{\text{CG}}$ is surjective, as follows from the known results on Raviart-Thomas-Nédélec spaces). Then we see that the auxiliary problem defined in (3.26) is the same as the auxiliary operator defined in (A.1) of [24]. For the right hand sides, we have gathered the two separate auxiliary problems of [24] into a single one, defined in (3.27a).

Everything is now in place so that we can apply Theorem A.1 in [24]. As explained above, we can extract the solution $\Lambda$ to the variational equation (3.28) from the solution of the mixed-hybrid equation (3.23). Moreover, the discrete bilinear form $a_{\tau}$ is defined as

$$ a_{\tau}(\Lambda, \hat{\Lambda}) = \begin{pmatrix} Q_{\Lambda}^{\text{m}} \\ Q_{\Lambda}^{\text{f}} \end{pmatrix}^T \begin{pmatrix} A^{\text{m}} & 0 \\ 0 & A^{\text{f}} \end{pmatrix} \begin{pmatrix} Q_{\hat{\Lambda}}^{\text{m}} \\ Q_{\hat{\Lambda}}^{\text{f}} \end{pmatrix} = (Q_{\Lambda}^{\text{m}})^T A^{\text{m}} Q_{\hat{\Lambda}}^{\text{m}} + (Q_{\Lambda}^{\text{f}})^T A^{\text{f}} Q_{\hat{\Lambda}}^{\text{f}}. \quad (3.31) $$

and the linear form $L$ is obtained in the same way:

$$ L(\hat{\Lambda}) = (\boldsymbol{P}_{\hat{\Lambda}}^{\text{m}})^T f^{\text{m}} + (\boldsymbol{Q}_{\hat{\Lambda}}^{\text{m}})^T g^{\text{m}} + (\boldsymbol{P}_{\hat{\Lambda}}^{\text{f}})^T f^{\text{f}} + (\boldsymbol{Q}_{\hat{\Lambda}}^{\text{f}})^T g^{\text{f}}. \quad (3.32) $$

Finally, (3.29) and (3.30) are simply translations from vectors to the discrete functions they represent. It is sufficient to observe that (3.31) and (3.29) are the same statement (similarly for the linear form in (3.30) and (3.32)).

# 4 Numerical convergence

We consider the test case 1 proposed by [15].

Let $\Omega$ be a cube of size 1. The domain $\Omega$ contains two orthogonal non-immersed fractures, denoted by $\Omega_1^f$ and $\Omega_2^f$. The first fracture $\Omega_1^f$ corresponds to the rectangle $\{(x, y, z) \in \Omega \mid x = 0.5\}$ and the second fracture $\Omega_2^f$ to the rectangle $\{(x, y, z) \in \Omega \mid y = 0.5\}$.

Dirichlet boundary conditions are imposed on the boundary of the domain $\Gamma_0$ and on the boundary of the fracture network $\Sigma_0$. The two fractures split the rock matrix domain into four subdomains $\omega_\alpha$, $\alpha \in \{1, 2, 3, 4\}$:

$$ \omega_1 = \{(x, y, z) \in \Omega \mid x < 0.5, y < 0.5\}, \quad \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\}, \quad \omega_4 = \{(x, y, z) \in \Omega \mid x < 0.5, y > 0.5\}. $$

Each fracture is divided into two sub-polygons along the intersection line, resulting in a network containing four sub-polygons $\gamma_\ell$, $\ell \in \{1, 2, 3, 4\}$. The transmissivity of the fracture network is
equal to $K_{\ell}^{\mathcal{F}} = \text{Id}, \ell \in \{1, 2, 3, 4\}$. The permeability of the rock matrix is equal to $K_{\alpha}^{\mathcal{M}} = \lambda_{\alpha} \text{Id}, \alpha \in \{1, 2, 3, 4\}$, where: $\lambda_1 = 1, \lambda_2 = 0.1, \lambda_3 = 0.01$ and $\lambda_4 = 10$.

We consider the function $p^{\mathcal{M}}(\mathbf{x}) = \exp(\cos(t_{\alpha}(\mathbf{x})))$, $\mathbf{x} \in \omega_{\alpha}, \alpha \in \{1, 2, 3, 4\}$, with:

$$t_1(\mathbf{x}) = y - x + z, \quad t_2(\mathbf{x}) = x + y + z - 1, \quad t_3(\mathbf{x}) = x - y + z, \quad t_4(\mathbf{x}) = 1 - x - y + z,$$

as the rock matrix hydraulic head analytical solution of this test case.

The fracture network surface mesh is generated by the software MODFRAC [14, 51] and the rock matrix volume mesh is performed with the software GHS3D [39]. We recall that we use conforming meshes, in the sense that the volume mesh is supported by the surface mesh.

The analytical solution is displayed on Figure 3.

<figure><img src="image_6.png" /><figcaption>Figure 3: Analytical solution on the rock matrix and fracture network meshes.</figcaption></figure>

Now let us consider a family of uniformly refined tetrahedral meshes, as described in Table 3.

<table><caption>Table 3: Mesh properties of the fractured porous medium presented in this test case.</caption><thead><tr><th>Mesh name</th><th>N(T<sub>h</sub>)</th><th>N(F<sub>h</sub>)</th><th>N(K<sub>h</sub>)</th><th>N(E<sub>h</sub>)</th><th>Linear system size</th><th>Non-zero elements</th></tr></thead><tbody><tr><td>1</td><td>96</td><td>216</td><td>16</td><td>30</td><td>182</td><td>1,114</td></tr><tr><td>2</td><td>7,850</td><td>16,345</td><td>444</td><td>696</td><td>15,671</td><td>107,128</td></tr><tr><td>3</td><td>54,190</td><td>111,053</td><td>1,782</td><td>2,733</td><td>108,280</td><td>747,263</td></tr><tr><td>4</td><td>351,738</td><td>714,383</td><td>7,282</td><td>11,043</td><td>703,292</td><td>4,879,505</td></tr><tr><td>5</td><td>813,284</td><td>1,643,678</td><td>11,438</td><td>17,307</td><td>1,626,365</td><td>11,316,300</td></tr></tbody></table>

The linear system (3.25) is solved with a direct solver, here LU factorization, as it has shown more robustness with respect to floating point arithmetic than the Cholesky decomposition in the MATLAB implementation.

Let us now compute the following relative errors on the mean hydraulic heads in the rock matrix and in the fracture network respectively (computed using quadrature rules on the tetrahedral mesh $T_h$ and on the triangular mesh $K_h$):

$$ \text{Err}_{p}^{\mathcal{M}} = \frac{\|p^{\mathcal{M}} - p_{h}^{\mathcal{M}}\|_{L^2(T_h)}}{\|p^{\mathcal{M}}\|_{L^2(T_h)}}, \qquad (4.1a) $$

$$ \text{Err}_{p}^{\mathcal{F}} = \frac{\|p^{\mathcal{F}} - p_{h}^{\mathcal{F}}\|_{L^2(K_h)}}{\|p^{\mathcal{F}}\|_{L^2(K_h)}}, \qquad (4.1b) $$
We also compute the following relative errors on the velocity in the rock matrix and on the velocity in the fracture network respectively (again, computed using quadrature rules on the tetrahedral mesh $T_h$ and on the triangular mesh $K_h$):

$$ \mathrm{Err}_u^m = \frac{\|u^m - u_h^m\|_{L^2(\mathcal{T}_h)}}{\|u^m\|_{L^2(\mathcal{T}_h)}}, \qquad (4.2a) $$

$$ \mathrm{Err}_u^f = \frac{\|u^f - u_h^f\|_{L^2(\mathcal{K}_h)}}{\|u^f\|_{L^2(\mathcal{K}_h)}}, \qquad (4.2b) $$

These errors are displayed on Figure 4.

<figure><img src="image_7.png" /><figcaption>Figure 4: Relative errors $Err_p$ and $Err_u$ computed on the four meshes presented in Table 3.</figcaption></figure>

The convergence rate between two successive meshes $k$ and $k+1$ are computed as in [15], both in the porous medium and in the fractures:

$$ \mathrm{CR}_p^{m,k+1} = 3 \frac{\ln \left( \frac{\| \mathrm{Err}_p^{m,k} \|}{\| \mathrm{Err}_p^{m,k+1} \|} \right)}{\ln \left( \frac{\#(\mathcal{N}(\mathcal{T}_h)^{k+1})}{\#(\mathcal{N}(\mathcal{T}_h)^k)} \right)}, $$

$$ \mathrm{CR}_p^{\mathcal{F},k+1} = 2 \frac{\ln \left( \frac{\| \mathrm{Err}_p^{\mathcal{F},k} \|}{\| \mathrm{Err}_p^{\mathcal{F},k+1} \|} \right)}{\ln \left( \frac{\#(\mathcal{N}(\mathcal{K}_h)^{k+1})}{\#(\mathcal{N}(\mathcal{K}_h)^k)} \right)}, $$

$$ \mathrm{CR}_u^{m,k+1} = 3 \frac{\ln \left( \frac{\| \mathrm{Err}_u^{m,k} \|}{\| \mathrm{Err}_u^{m,k+1} \|} \right)}{\ln \left( \frac{\#(\mathcal{N}(\mathcal{T}_h)^{k+1})}{\#(\mathcal{N}(\mathcal{T}_h)^k)} \right)}, $$

$$ \mathrm{CR}_u^{\mathcal{F},k+1} = 2 \frac{\ln \left( \frac{\| \mathrm{Err}_u^{\mathcal{F},k} \|}{\| \mathrm{Err}_u^{\mathcal{F},k+1} \|} \right)}{\ln \left( \frac{\#(\mathcal{N}(\mathcal{K}_h)^{k+1})}{\#(\mathcal{N}(\mathcal{K}_h)^k)} \right)}. $$

The relative errors and convergence order in the rock matrix medium (resp. the fracture network medium) are reported in Table 4 (resp. in Table 5).

<table><caption>Table 4: Relative errors and orders of convergence in the rock matrix domain.</caption><thead><tr><th>Mesh name</th><th>N(T<sub>h</sub>)</th><th>Err<sub>p</sub><sup>m</sup></th><th>Err<sub>u</sub><sup>m</sup></th><th>CR<sub>p</sub><sup>m,k+1</sup></th><th>CR<sub>u</sub><sup>m,k+1</sup></th></tr></thead><tbody><tr><td>1</td><td>96</td><td>6.50 &middot; 10<sup>-2</sup></td><td>2.52 &middot; 10<sup>-1</sup></td><td>n/a</td><td>n/a</td></tr><tr><td>2</td><td>7,850</td><td>1.46 &middot; 10<sup>-2</sup></td><td>5.48 &middot; 10<sup>-2</sup></td><td>1.02</td><td>1.04</td></tr><tr><td>3</td><td>54,190</td><td>7.89 &middot; 10<sup>-3</sup></td><td>2.92 &middot; 10<sup>-2</sup></td><td>0.96</td><td>0.98</td></tr><tr><td>4</td><td>351,738</td><td>4.47 &middot; 10<sup>-3</sup></td><td>1.63 &middot; 10<sup>-2</sup></td><td>0.91</td><td>0.93</td></tr><tr><td>5</td><td>813,284</td><td>3.19 &middot; 10<sup>-3</sup></td><td>1.17 &middot; 10<sup>-2</sup></td><td>1.20</td><td>1.19</td></tr></tbody></table>
<table><caption>Table 5: Relative errors and orders of convergence in the fracture network.</caption><thead><tr><th>Mesh name</th><th>N(K<sub>h</sub>)</th><th>Err<sub>p</sub><sup>f</sup></th><th>Err<sub>u</sub><sup>f</sup></th><th>CR<sub>p</sub><sup>f,k+1</sup></th><th>CR<sub>u</sub><sup>f,k+1</sup></th></tr></thead><tbody><tr><td>1</td><td>16</td><td>7.20 &middot; 10<sup>-2</sup></td><td>2.32 &middot; 10<sup>-1</sup></td><td>n/a</td><td>n/a</td></tr><tr><td>2</td><td>444</td><td>1.42 &middot; 10<sup>-2</sup></td><td>4.86 &middot; 10<sup>-2</sup></td><td>0.98</td><td>0.94</td></tr><tr><td>3</td><td>1,782</td><td>6.98 &middot; 10<sup>-3</sup></td><td>2.40 &middot; 10<sup>-2</sup></td><td>1.02</td><td>1.02</td></tr><tr><td>4</td><td>7,282</td><td>3.43 &middot; 10<sup>-3</sup></td><td>1.18 &middot; 10<sup>-2</sup></td><td>1.01</td><td>1.01</td></tr><tr><td>5</td><td>11,438</td><td>2.73 &middot; 10<sup>-3</sup></td><td>9.39 &middot; 10<sup>-3</sup></td><td>1.01</td><td>1.00</td></tr></tbody></table>

Figure 4, along with Tables 4 and 5, illustrates a first-order convergence rate, in agreement with Theorem 3.5. It is important to note that this example exclusively features fractures that span the entire domain, thereby satisfying the smoothness assumptions of Theorem 3.5. Future research will investigate configurations that include immersed fractures. Based on the examples presented in [13], we anticipate observing a reduction in the convergence rate in these more complex scenarios.

# 5 Conclusion

In this work, we investigated various aspects of simulating single-phase flow in three-dimensional fractured porous media using the 2D-3D Discrete Fracture Matrix (DFM) model. We defined the geometrical setting, where each 2D fracture is represented as a union of 2D polygons, inducing a partition of the 3D rock matrix domain into connected components. The physical problem involves coupled Darcy-type PDEs with coupling conditions across fractures and at their intersections. The model assumes the continuity of the hydraulic head, making it relevant for conductive fractures.

We derived a weak formulation and proposed a functional setting that makes this formulation well-posed. Then we went through the discrete setting thanks to the mixed finite element approximation in conforming 2D-3D meshes, by defining the finite element spaces and proving first the existence and uniqueness of the discrete solution, and second by showing first order error estimates between the continuous and the discrete solutions.

By leveraging the equivalence with a mixed-hybrid formulation, we derived the algebraic formulation and then the global linear system. The linear system has as unknowns the traces of the hydraulic head (Lagrange multipliers) on element faces in both the rock matrix and the fracture network, and the mean hydraulic head on each element in the fracture network. The resulting linear system is square, sparse, symmetric and positive definite. The remaining unknowns, such as the hydraulic head in the rock matrix and the velocities in both media, are reconstructed in a post-processing phase. We also showed that the mixed-hybrid formulation is equivalent to a mesh-dependent variational formulation, with a discrete bilinear form computable element by element.

We concluded with a numerical convergence test, comparing the computed solution with an analytical solution. The test confirmed that the method converges to first order both in the rock matrix and in the fractures, at least in the absence of immersed fractures.

The numerical scheme analyzed in the present paper has been applied to the simulation of flow in large fractured porous media [49, 46] (containing up to 697,000 fractures in [46]).

Future work will aim at relaxing the smoothness requirements on the solution in order to derive an improved error estimate. In this direction, we intend to exploit recent advances on commuting projections that require minimal smoothness [31, Chapter 23] to replace the interpolation operators. We also expect that recent work [65, 30] may allow us to replace the global norms of the solution in the error estimates by a sum of norms over the individual rock matrix subdomains and the fractures, where the solution may exhibit more regularity.
Further research should also address the influence of fracture tips on the convergence rate and investigate the actual regularity of the continuous solution, both locally (in the rock matrix subdomains and in the fracture network) and globally.

## Acknowledgements

The authors thank Martin Vohralík for helpful discussions.

# Appendices

## A Derivation of the weak formulation

We proceed to establish the weak formulation of the problem described in equations (1.2a) to (1.5b). The derivation here is formal, and we assume that all functions are smooth enough so that the formal manipulations are justified. The precise choice of the function spaces is given in Subsection 2.2.

We first deal with the porous rock matrix flow problem.

We start with the continuity equation (1.2a), where no integration by parts is required.

$$ \sum_{\alpha \in \mathcal{C}} \int_{\omega_{\alpha}} \nabla \cdot u_{\alpha}^{\mathfrak{m}} \Psi_{\alpha}^{\mathfrak{m}} = \int_{\Omega^{\mathfrak{m}}} f^{\mathfrak{m}} \Psi^{\mathfrak{m}}, \quad \forall \Psi^{\mathfrak{m}} \in L^2(\Omega^{\mathfrak{m}}). \quad (A.1) $$

Then we turn to Darcy's law (1.2b). For a given $\alpha \in \mathcal{C}$ and test functions $v_{\alpha}^{\mathfrak{m}}$ such that $\gamma_{n}^{\mathfrak{m}, \Gamma_{\alpha}^{N}} v_{\alpha}^{\mathfrak{m}} = 0$ on $\Gamma_{\alpha}^{N}$, we multiply equation (1.2b) by $v_{\alpha}^{\mathfrak{m}}$, integrate over $\omega_{\alpha}$ and use Green's formula. After splitting the boundary as $\partial\omega_{\alpha} = \Gamma_{\alpha}^{D} \cup \Gamma_{\alpha}^{N} \cup \Gamma_{\alpha}^{\gamma} \cup \Gamma_{\alpha}^{\text{fict}}$ and using the boundary conditions, there holds:

$$ \int_{\omega_{\alpha}} (K_{\alpha}^{\mathfrak{m}})^{-1} u_{\alpha}^{\mathfrak{m}} \cdot v_{\alpha}^{\mathfrak{m}} - \int_{\omega_{\alpha}} p_{\alpha}^{\mathfrak{m}} \nabla \cdot v_{\alpha}^{\mathfrak{m}} + \sum_{\ell \in \mathcal{I}_{\alpha}} \int_{\gamma_{\ell}} p_{\alpha}^{\mathfrak{m}} \gamma_{n}^{\mathfrak{m}, \gamma_{\ell}} v_{\alpha}^{\mathfrak{m}} + \sum_{\gamma \in \Gamma_{\alpha}^{\text{fict}}} \int_{\gamma} p_{\alpha}^{\mathfrak{m}} \gamma_{n}^{\mathfrak{m}, \gamma} v_{\alpha}^{\mathfrak{m}} = - \int_{\Gamma_{\alpha}^{D}} g_{\alpha}^{\mathfrak{m}, D} \gamma_{n}^{\mathfrak{m}, \Gamma_{\alpha}^{D}} v_{\alpha}^{\mathfrak{m}}. \quad (A.2) $$

We then sum (A.2) over $\alpha$, and reorder the double sums that appear in the last two terms on the left-hand side so as to make the jump term (1.1) appear. We then use the continuity condition (1.4) for the fractures:

$$ \sum_{\alpha \in \mathcal{C}} \sum_{\ell \in \mathcal{I}_{\alpha}} \int_{\gamma_{\ell}} p_{\alpha}^{\mathfrak{m}} \gamma_{n}^{\mathfrak{m}, \gamma_{\ell}} v_{\alpha}^{\mathfrak{m}} = \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} p_{\ell}^{\mathfrak{f}} [[v^{\mathfrak{m}}]]|_{\gamma_{\ell}}, $$

while the corresponding term across a fictitious face disappears, since the jump across such a face is zero.

Finally, we obtain, for any sufficiently smooth $v^{\mathfrak{m}} = (v_{\alpha}^{\mathfrak{m}})_{\alpha \in \mathcal{C}}$:

$$ \int_{\Omega^{\mathfrak{m}}} (K^{\mathfrak{m}})^{-1} u^{\mathfrak{m}} \cdot v^{\mathfrak{m}} - \sum_{\alpha \in \mathcal{C}} \int_{\omega_{\alpha}} p_{\alpha}^{\mathfrak{m}} \nabla \cdot v_{\alpha}^{\mathfrak{m}} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} p_{\ell}^{\mathfrak{f}} [[v^{\mathfrak{m}}]]|_{\gamma_{\ell}} = - \int_{\Gamma^{D}} g^{\mathfrak{m}, D} \gamma_{n}^{\mathfrak{m}, \Gamma^{D}} v^{\mathfrak{m}}. \quad (A.3) $$

We proceed in a similar way for the fracture network flow problem.

The continuity equation (1.3a) gives:

$$ \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} \nabla_{\tau} \cdot u_{\ell}^{\mathfrak{f}} \Psi_{\ell}^{\mathfrak{f}} = \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} [[u^{\mathfrak{m}}]]|_{\gamma_{\ell}} \Psi_{\ell}^{\mathfrak{f}} + \int_{\Omega^{\mathfrak{f}}} f^{\mathfrak{f}} \Psi^{\mathfrak{f}}, \quad \forall \Psi^{\mathfrak{f}} \in L^2(\Omega^{\mathfrak{f}}). \quad (A.4) $$
Next we turn to Darcy's law (1.3b). For a given $\ell \in \mathcal{I}$ and test functions $v_{\ell}^{\mathcal{I}}$ such that $\gamma_{n}^{\mathcal{I}} v_{\ell}^{\mathcal{I}} = 0$ on $\Sigma_{\ell}^{N}$, we multiply (1.3b) by $v_{\ell}^{\mathcal{I}}$, use Green's formula, split the boundary as $\partial\gamma_{\ell} = \Sigma_{\ell} = \Sigma_{\ell}^{D} \cup \Sigma_{\ell}^{N} \cup \Sigma_{\ell}^{\text{inter}}$ and use the boundary conditions. There holds:

$$ \int_{\gamma_{\ell}} (\mathcal{K}_{\ell}^{\mathcal{I}})^{-1} u_{\ell}^{\mathcal{I}} \cdot v_{\ell}^{\mathcal{I}} - \int_{\gamma_{\ell}} p_{\ell}^{\mathcal{I}} \nabla_{\tau} \cdot v_{\ell}^{\mathcal{I}} + \sum_{\sigma \subset \Sigma_{\ell}^{\text{inter}}} \int_{\sigma} \gamma_{\ell}^{\mathcal{I},\sigma} p_{\ell}^{\mathcal{I}} \gamma_{n}^{\mathcal{I},\sigma} v_{\ell}^{\mathcal{I}} = - \int_{\Sigma_{\ell}^{D}} g_{\ell}^{\mathcal{I},D} \gamma_{n}^{\mathcal{I},\Sigma_{\ell}^{D}} v_{\ell}^{\mathcal{I}}. \quad (A.5) $$

We then sum (A.5) over the fractures $\ell \in \mathcal{I}$:

$$ \int_{\Omega^{\mathcal{I}}} (\mathcal{K}^{\mathcal{I}})^{-1} u^{\mathcal{I}} \cdot v^{\mathcal{I}} - \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} p_{\ell}^{\mathcal{I}} \nabla_{\tau} \cdot v_{\ell}^{\mathcal{I}} + \sum_{\ell \in \mathcal{I}} \sum_{\sigma \subset \Sigma_{\ell}^{\text{inter}}} \int_{\sigma} \gamma_{\ell}^{\mathcal{I},\sigma} p_{\ell}^{\mathcal{I}} \gamma_{n}^{\mathcal{I},\sigma} v_{\ell}^{\mathcal{I}} = - \sum_{\ell \in \mathcal{I}} \int_{\Sigma_{\ell}^{D}} g_{\ell}^{\mathcal{I},D} \gamma_{n}^{\mathcal{I},\Sigma_{\ell}^{D}} v_{\ell}^{\mathcal{I}}. \quad (A.6) $$

Because of the continuity condition (1.5a), the traces on $\sigma$ of all the functions $p_{\ell}^{\mathcal{I}}$, $\ell \in \mathcal{I}_{\sigma}$, take a common value that we denote by $p^{\sigma}$. We then reorder the two sums over fractures and intersections to obtain:

$$ \int_{\Omega^{\mathcal{I}}} (\mathcal{K}^{\mathcal{I}})^{-1} u^{\mathcal{I}} \cdot v^{\mathcal{I}} - \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} p_{\ell}^{\mathcal{I}} \nabla_{\tau} \cdot v_{\ell}^{\mathcal{I}} + \sum_{\sigma \subset \Sigma} \int_{\sigma} p^{\sigma} \left( \sum_{\ell \in \mathcal{I}_{\sigma}} \gamma_{n}^{\mathcal{I},\sigma} v_{\ell}^{\mathcal{I}} \right) = - \sum_{\ell \in \mathcal{I}} \int_{\Sigma_{\ell}^{D}} g_{\ell}^{\mathcal{I},D} \gamma_{n}^{\mathcal{I},\Sigma_{\ell}^{D}} v_{\ell}^{\mathcal{I}}. $$

Because of condition (1.5b), the last term on the left-hand side is zero, so that we finally obtain, for any sufficiently smooth $v^{\mathcal{I}} = (v_{\ell}^{\mathcal{I}})_{\ell \in \mathcal{I}}$:

$$ \int_{\Omega^{\mathcal{I}}} (\mathcal{K}^{\mathcal{I}})^{-1} u^{\mathcal{I}} \cdot v^{\mathcal{I}} - \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} p_{\ell}^{\mathcal{I}} \nabla_{\tau} \cdot v_{\ell}^{\mathcal{I}} = - \sum_{\ell \in \mathcal{I}} \int_{\Sigma_{\ell}^{D}} g_{\ell}^{\mathcal{I},D} \gamma_{n}^{\mathcal{I},\Sigma_{\ell}^{D}} v_{\ell}^{\mathcal{I}}. \quad (A.7) $$

The weak formulation of the coupled problem is obtained by first summing the continuity equation (A.1) for the porous rock matrix domain with its counterpart (A.4) obtained for the fractures:

$$ \sum_{\alpha \in \mathcal{C}} \int_{\omega_{\alpha}} \nabla \cdot u_{\alpha}^{\mathcal{I}} \Psi_{\alpha}^{\mathcal{I}} + \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} (\nabla_{\tau} \cdot u_{\ell}^{\mathcal{I}} - [[u^{\mathcal{I}}]]_{|\gamma_{\ell}}) \Psi_{\ell}^{\mathcal{I}} = \int_{\Omega^{\mathcal{I}}} f^{\mathcal{I}} \Psi^{\mathcal{I}} + \int_{\Omega^{\mathcal{I}}} f^{\mathcal{I}} \Psi^{\mathcal{I}}, \quad \forall \Psi^{\mathcal{I}} \in L^2(\Omega^{\mathcal{I}}), \forall \Psi^{\mathcal{I}} \in L^2(\Omega^{\mathcal{I}}). \quad (A.8) $$

We proceed similarly for Darcy's law, equations (A.3) and (A.7)

$$ \int_{\Omega^{\mathcal{I}}} (\mathcal{K}^{\mathcal{I}})^{-1} u^{\mathcal{I}} \cdot v^{\mathcal{I}} + \int_{\Omega^{\mathcal{I}}} (\mathcal{K}^{\mathcal{I}})^{-1} u^{\mathcal{I}} \cdot v^{\mathcal{I}} - \sum_{\alpha \in \mathcal{C}} \int_{\omega_{\alpha}} p_{\alpha}^{\mathcal{I}} \nabla \cdot v_{\alpha}^{\mathcal{I}} - \sum_{\ell \in \mathcal{I}} \int_{\gamma_{\ell}} p_{\ell}^{\mathcal{I}} (\nabla_{\tau} \cdot v_{\ell}^{\mathcal{I}} - [[v^{\mathcal{I}}]]_{|\gamma_{\ell}}) \\ = - \int_{\Gamma^{\mathcal{I}}} g_{\ell}^{\mathcal{I},D} \gamma_{n}^{\mathcal{I},\Gamma^{\mathcal{I}}} v^{\mathcal{I}} - \sum_{\ell \in \mathcal{I}} \int_{\Sigma_{\ell}^{D}} g_{\ell}^{\mathcal{I},D} \gamma_{n}^{\mathcal{I},\Sigma_{\ell}^{D}} v_{\ell}^{\mathcal{I}}, \quad \forall v^{\mathcal{I}}, \forall v^{\mathcal{I}}. \quad (A.9) $$

## B Matrix formulation of the mixed-hybrid linear system

In this appendix, we gather the derivation of the reduced linear system (3.25) from the mixed-hybrid formulation (3.23).

We proceed in two steps: we first eliminate the fluxes to obtain a first reduced system, and then we further eliminate the hydraulic head unknowns in the rock matrix (but not those in the fracture network) to obtain the final system. Both steps follow the standard procedure used for equi-dimensional equations, with specific adaptations to take into account the coupling between the rock matrix and the fracture network.
## B.1 The mixed-hybrid formulation: algebraic equations

It is convenient to define the following sets that will let us identify the entities involved in various equation (some sets were already defined and used in Section 3.1.3):

* for a given tetrahedron $T \in \mathcal{T}_h$, we let $\mathcal{F}(T) = \{F \in \mathcal{F}_h | F \text{ is a face of } T\}$;

* for a given face $F \in \mathcal{F}_h$, we let $\mathcal{T}(F) = \{T \in \mathcal{T}_h | F \text{ is a face of } T\}$;

* for a given triangle $K \in \mathcal{K}_h$, we let $\mathcal{T}(K) = \{T \in \mathcal{T}_h | K \text{ is a fracture face of } T\}$;

* for a given triangle $K \in \mathcal{K}_h$, we let $\mathcal{E}(K) = \{E \in \mathcal{E}_h | E \text{ is an edge of } K\}$;

* for a given edge $E \in \mathcal{E}_h$, we let $\mathcal{K}(E) = \{K \in \mathcal{K}_h | E \text{ is an edge of } K\}$.

Furthermore, we let $\mathcal{F}(T)^{\text{frac}} = \mathcal{F}(T) \cap \mathcal{F}_h^{\text{frac}}$ and $\mathcal{F}(T)^{\text{in,N}} = \mathcal{F}(T) \cap \mathcal{F}_h^{\text{in,N}}$, as well as $\mathcal{T}(F, F') = \mathcal{T}(F) \cap \mathcal{T}(F')$ (where $K$ can be substituted by $F$) and $\mathcal{K}(E, E') = \mathcal{K}(E) \cap \mathcal{K}(E')$.

Note that the sets $\mathcal{T}(F)$ usually contain two elements, unless the face is on the boundary of the domain. On the other hand $\mathcal{T}(K)$ always contains two elements, as it is assumed that no fracture lies along the boundary of the cube.

Because of the block structures of all matrices involved in (3.23), it is possible to write the equations at the element level. We note that this enables us to recover a formulation that is identical to the formulation described in [44].

**In the rock matrix.** We write the discrete version of Darcy's law and the local mass conservation for each tetrahedron:

$$A_T^m \mathbf{q}_T^m - p_T^m e^m - \sum_{K \in \mathcal{F}(T)^{\text{frac}}} B_{T,K}^{m,\mathcal{F}} p_K^\mathcal{F} + \sum_{F \in \mathcal{F}(T)^{\text{in,N}}} C_{T,F}^m \lambda_F^m = g_T^m, \quad \forall T \in \mathcal{T}_h. \quad (B.1a)$$

$$(e^m)^T \mathbf{q}_T^m = f_T^m, \quad \forall T \in \mathcal{T}_h. \quad (B.1b)$$

and the flux continuity across interior faces gives:

$$\sum_{T \in \mathcal{T}(F)} (C_{T,F}^m)^T \mathbf{q}_T^m = 0, \quad \forall F \in \mathcal{F}_h^{\text{in,N}}. \quad (B.1c)$$

**In the fracture network.** We write the discrete version of Darcy's law and the local mass conservation for each triangle:

$$A_K^\mathcal{F} \mathbf{q}_K^\mathcal{F} - p_K^\mathcal{F} e^\mathcal{F} + \sum_{E \in \mathcal{E}(K)} C_{K,E}^\mathcal{F} \lambda_E^\mathcal{F} = g_K^\mathcal{F}, \quad \forall K \in \mathcal{K}_h, \quad (B.2a)$$

$$(e^\mathcal{F})^T \mathbf{q}_K^\mathcal{F} + \sum_{T \in \mathcal{T}(K)} (B_{T,K}^{\mathcal{F},\mathcal{F}})^T \mathbf{q}_T^m = f_K^\mathcal{F}, \quad \forall K \in \mathcal{K}_h, \quad (B.2b)$$

and the flux continuity across fracture intersections gives:

$$\sum_{K \in \mathcal{K}(E)} (C_{K,E}^\mathcal{F})^T \mathbf{q}_K^\mathcal{F} = 0, \quad \forall E \in \mathcal{E}_h \setminus \mathcal{E}_h^D. \quad (B.2c)$$

## B.2 Elimination of the fluxes

As usual for hybridized formulations, it is now possible to eliminate the flux unknowns at the element level, using equations (B.1a) and (B.2a). Note also that these two equations will be used to recover the fluxes after the global linear system (3.25) has been solved. In the more usual case of an equi-dimensional model, the matrix $\mathcal{B}$ is also block diagonal, but this is not the case here due to the coupling between the rock matrix medium and the fracture network medium.

We give below the system of equations obtained after the elimination, as well as the value of the matrix elements and right-hand side vectors. Only the non-zero elements are defined.
Mean hydraulic head in the rock matrix.

$$D_T^m p_T^m - \sum_{F \in \mathcal{F}(T)^{\text{in,N}}} R_{T,F}^{m,\mathcal{F}} \lambda_F^m + \sum_{K \in \mathcal{F}(T)^{\text{frac}}} R_{T,K}^{m,\mathcal{F}} p_K^{\mathcal{F}} = \mathcal{F}_T^m, \quad \forall T \in \mathcal{T}_h. \quad (B.3)$$

with

$$D_T^m = (e^m)^T (A_T^m)^{-1} e^m = \sum_{(F,F') \in \mathcal{F}(T)^2} ((A_T^m)^{-1})_{FF'}, \quad \forall T \in \mathcal{T}_h$$
$$R_{TF}^{m,\mathcal{F}} = (e^m)^T (A_T^m)^{-1} C_{TF}^m = \sum_{F' \in \mathcal{F}(T)} ((A_T^m)^{-1})_{FF'} \quad \forall T \in \mathcal{T}_h, \forall F \in \mathcal{F}_h^{\text{in,N}},$$
$$R_{TK}^{m,\mathcal{F}} = -(e^m)^T (A_T^m)^{-1} B_{TK}^{m,\mathcal{F}} = \sum_{F \in \mathcal{F}(T)} ((A_T^m)^{-1})_{K,F}, \quad \forall T \in \mathcal{T}_h, \forall K \in \mathcal{K}_h,$$

and

$$\mathcal{F}_T^m = f_T^m - (e^m)^T (A_T^m)^{-1} g_T^m, \quad \forall T \in \mathcal{T}_h.$$

Lagrange multipliers in the rock matrix.

$$- \sum_{T \in \mathcal{T}(F)} (R_{T,F}^{m,\mathcal{F}})^T p_T^m + \sum_{F' \in \mathcal{F}_h^{\text{in,N}}} M_{F,F'}^{m,\mathcal{F}} \lambda_{F'}^m + \sum_{K \in \mathcal{K}_h} M_{F,K}^{m,\mathcal{F}} p_K^{\mathcal{F}} = \mathcal{Y}_F^m, \quad \forall F \in \mathcal{F}_h. \quad (B.4)$$

with

$$M_{F,F'}^{m,\mathcal{F}} = \sum_{T \in \mathcal{T}(F,F')} (C_{TF}^{m})^T (A_T^m)^{-1} C_{TF'}^{m'} = \sum_{T \in \mathcal{T}(F,F')} ((A_T^m)^{-1})_{FF'} \quad \forall (F,F') \in (\mathcal{F}_h^{\text{in,N}})^2,$$
$$M_{FK}^{m,\mathcal{F}} = - \sum_{T \in \mathcal{T}(F,K)} (C_{TF}^{m})^T (A_T^m)^{-1} B_{TK}^{m,\mathcal{F}} = \sum_{T \in \mathcal{T}(F,K)} (A_T^m)^{-1} {}_{FK} \quad \forall F \in \mathcal{F}_h^{\text{in,N}}, \forall K \in \mathcal{K}_h.$$

and

$$\mathcal{Y}_F^m = \sum_{T \in \mathcal{T}(F)} (C_{TF}^m)^T (A_T^m)^{-1} g_T^m, \quad \forall F \in \mathcal{F}_h^{\text{in,N}}.$$

Mean hydraulic head in the fracture network.

$$- \sum_{T \in \mathcal{T}(K)} (R_{T,K}^{m,\mathcal{F}})^T p_T^m + \sum_{F \in \mathcal{F}_h^{\text{in,N}}} (M_{F,K}^{m,\mathcal{F}})^T \lambda_F^m + D_K^{\mathcal{F}} p_K^{\mathcal{F}} + \sum_{K' \in \mathcal{K}_h} M_{K,K'}^{m,\mathcal{F}} p_{K'}^{\mathcal{F}} \\ - \sum_{E \in \mathcal{E}(K)} R_{K,E}^{\mathcal{F}} \lambda_E^{\mathcal{F}} = \mathcal{F}_K^{\mathcal{F}} + \mathcal{Y}_K^{m,\mathcal{F}}, \quad \forall K \in \mathcal{K}_h. \quad (B.5)$$

with

$$D_K^{\mathcal{F}} = (e^{\mathcal{F}})^T (A_K^{\mathcal{F}})^{-1} e^{\mathcal{F}} = \sum_{(E,E') \in \mathcal{E}(K)^2} ((A_K^{\mathcal{F}})^{-1})_{E,E'} \quad \forall K \in \mathcal{K}_h,$$
$$M_{K,K'}^{\mathcal{F},\mathcal{F}} = \sum_{T \in \mathcal{T}(K,K')} (B_{TK}^{m,\mathcal{F}})^T (A_T^m)^{-1} B_{TK'}^{m,\mathcal{F}} = \sum_{T \in \mathcal{T}(K,K')} ((A_T^m)^{-1})_{K,K'} \quad \forall (K,K') \in (\mathcal{K}_h)^2,$$
$$R_{KE}^{\mathcal{F}} = (e^{\mathcal{F}})^T (A_K^{\mathcal{F}})^{-1} C_{KE}^{\mathcal{F}} = \sum_{E' \in \mathcal{E}(K)} ((A_K^{\mathcal{F}})^{-1})_{E,E'} \quad \forall K \in \mathcal{K}_h, \forall E \in \mathcal{E}_h.$$

and

$$\mathcal{F}_K^{\mathcal{F}} = f_K^{\mathcal{F}} + (e^{\mathcal{F}})^T (A_K^{\mathcal{F}})^{-1} g_K^{\mathcal{F}}, \quad \forall K \in \mathcal{K}_h,$$
$$\mathcal{Y}_K^{m,\mathcal{F}} = - \sum_{T \in \mathcal{T}(K)} (B_{TK}^{m,\mathcal{F}})^T (A_T^m)^{-1} g_T^m, \quad \forall K \in \mathcal{K}_h.$$
Lagrange multipliers in the fracture network.

$$ - \sum_{K \in \mathcal{K}(E)} (R_{K,E}^{\mathcal{F}})^T p_K^{\mathcal{F}} + \sum_{E' \in \mathcal{E}_h \setminus \mathcal{E}_h^D} M_{E,E'}^{\mathcal{F}} \lambda_{E'}^{\mathcal{F}} = \gamma_E^{\mathcal{F}}, \quad \forall E \in \mathcal{E}_h \setminus \mathcal{E}_h^D. \quad (B.6) $$

with

$$ M_{E,E'}^{\mathcal{F}} = \sum_{K \in \mathcal{K}(E,E')} (C_{KE}^{\mathcal{F}})^T (A_K^{\mathcal{F}})^{-1} C_{KE'}^{\mathcal{F}} = \sum_{K \in \mathcal{K}(E,E')} ((A_K^{\mathcal{F}})^{-1})_{E,E'}, \quad \forall (E,E') \in (\mathcal{E}_h)^2. $$

and

$$ \gamma_E^{\mathcal{F}} = \sum_{K \in \mathcal{K}(E)} (C_{KE}^{\mathcal{F}})^T (A_K^{\mathcal{F}})^{-1} g_K^{\mathcal{F}}, \quad \forall E \in \mathcal{E}_h \setminus \mathcal{E}_h^D. $$

The matrices defined above have the following dimensions:

$$ \begin{align*} D^{\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{T}_h) \times \mathcal{N}(\mathcal{T}_h)}, & R^{\mathcal{F},\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{T}_h) \times \mathcal{N}(\mathcal{T}_h^{\text{in},N})}, \\ R^{\mathcal{F},\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{T}_h) \times \mathcal{N}(\mathcal{K}_h)}, & M^{\mathcal{F},\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{F}_h^{\text{in},N}) \times \mathcal{N}(\mathcal{F}_h^{\text{in},N})}, \\ M^{\mathcal{F},\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{K}_h) \times \mathcal{N}(\mathcal{F}_h^{\text{in},N})}, & D^{\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{K}_h) \times \mathcal{N}(\mathcal{K}_h)}, \\ M^{\mathcal{F},\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{K}_h) \times \mathcal{N}(\mathcal{K}_h)}, & R^{\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{K}_h) \times \mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)}, \\ M^{\mathcal{F}} &\in \mathbb{R}^{\mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D) \times \mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)}. \end{align*} $$

At the matrix level, after reordering the unknowns to group the rock matrix and fracture network unknowns together, one obtains a reduced system for the hydraulic head and Lagrange multiplier unknowns:

$$ \begin{pmatrix} D^{\mathcal{F}} & -R^{\mathcal{F},\mathcal{F}} & -R^{\mathcal{F},\mathcal{F}} & 0 \\ -(R^{\mathcal{F},\mathcal{F}})^T & M^{\mathcal{F},\mathcal{F}} & M^{\mathcal{F},\mathcal{F}} & 0 \\ -(R^{\mathcal{F},\mathcal{F}})^T & (M^{\mathcal{F},\mathcal{F}})^T & M^{\mathcal{F},\mathcal{F}} + D^{\mathcal{F}} & -R^{\mathcal{F}} \\ 0 & 0 & -(R^{\mathcal{F}})^T & M^{\mathcal{F}} \end{pmatrix} \begin{pmatrix} \mathbf{P}^{\mathcal{F}} \\ \mathbf{\Lambda}^{\mathcal{F}} \\ \mathbf{P}^{\mathcal{F}} \\ \mathbf{\Lambda}^{\mathcal{F}} \end{pmatrix} = \begin{pmatrix} \mathcal{F}^{\mathcal{F}} \\ \gamma^{\mathcal{F}} \\ \gamma^{\mathcal{F},\mathcal{F}} + \mathcal{F}^{\mathcal{F}} \\ \gamma^{\mathcal{F}} \end{pmatrix}. \quad (B.7) $$

The left-hand side matrix in (B.7) is square and sparse. It can be shown to be symmetric and positive definite (see [43, 25]).

## B.3 Elimination of the hydraulic head in the rock matrix

The matrix $D^{\mathcal{F}}$ is again diagonal, so that the unknowns $\mathbf{P}^{\mathcal{F}}$ can be eliminated at the element level. As explained above, this is not true for the $\mathbf{P}^{\mathcal{F}}$ unknowns. After elimination, the following linear system is obtained:

$$ \begin{pmatrix} A^{\mathcal{F},\mathcal{F}} & A^{\mathcal{F},\mathcal{F}} & 0 \\ (A^{\mathcal{F},\mathcal{F}})^T & A^{\mathcal{F},\mathcal{F}} + D^{\mathcal{F}} & -R^{\mathcal{F}} \\ 0 & -(R^{\mathcal{F}})^T & M^{\mathcal{F}} \end{pmatrix} \begin{pmatrix} \mathbf{\Lambda}^{\mathcal{F}} \\ \mathbf{P}^{\mathcal{F}} \\ \mathbf{\Lambda}^{\mathcal{F}} \end{pmatrix} = \begin{pmatrix} \tilde{\mathcal{F}}^{\mathcal{F}} + \gamma^{\mathcal{F}} \\ \tilde{\mathcal{F}}^{\mathcal{F},\mathcal{F}} + V^{\mathcal{F},\mathcal{F}} + \mathcal{F}^{\mathcal{F}} \\ \gamma^{\mathcal{F}} \end{pmatrix}, \quad (B.8) $$

As it the Schur complement of the left-hand side matrix in (B.7), the matrix in (B.8) is also square, sparse, symmetric and positive definite.
The blocks are given by:

$$A_{F,F'}^{m,m} = M_{F,F'}^{m,m} - \sum_{T \in \mathcal{T}(F,F')} (R_{T,F}^{m,m})^T (D_T^m)^{-1} R_{T,F'}^{m,m}, \quad \forall (F,F') \in (\mathcal{F}_h^{\text{in},N})^2,$$

$$A_{F,K}^{m,f} = M_{F,K}^{m,f} - \sum_{T \in \mathcal{T}(F,K)} (R_{T,F}^{m,m})^T (D_T^m)^{-1} R_{T,K}^{m,f}, \quad \forall F \in \mathcal{F}_h^{\text{in},N}, \forall K \in \mathcal{K}_h,$$

$$A_{K,K'}^{f,f} = M_{K,K'}^{f,f} - \sum_{T \in \mathcal{T}(K,K')} (R_{T,K}^{m,f})^T (D_T^m)^{-1} R_{T,K'}^{m,f}, \quad \forall (K,K') \in (\mathcal{K}_h)^2,$$

$$\tilde{\mathcal{F}}_F^m = \sum_{T \in \mathcal{T}(F)} (R_{T,F}^{m,m})^T (D_T^m)^{-1} \mathcal{F}_T^m, \quad \forall F \in \mathcal{F}_h^{\text{in},N},$$

$$\tilde{\mathcal{F}}_K^{m,f} = \sum_{T \in \mathcal{T}(K)} (R_{T,K}^{m,f})^T (D_T^m)^{-1} \mathcal{F}_K^f, \quad \forall K \in \mathcal{K}_h.$$

## References

[1] C. Alboin. Deux outils mathématiques pour modéliser l'écoulement et le transport de polluants dans un milieu poreux fracturé. Phd thesis (in french), Paris IX Dauphine University (France), 2000. URL: https://bu.dauphine.psl.eu/fileviewer/view.php?doc=2000PA090024.

[2] C. Alboin, J. Jaffré, J. E. Roberts, and C. Serres. Modeling fractures as interfaces for flow and transport in porous media. In *Fluid flow and transport in porous media: mathematical and numerical treatment (South Hadley, MA, 2001)*, volume 295 of *Contemp. Math.*, pages 13–24. Amer. Math. Soc., Providence, RI, 2002. doi:10.1090/conm/295/04999.

[3] L. Amir, M. Kern, Z. Mghazli, and J. E. Roberts. Intersecting fractures in porous media: mathematical and numerical analysis. *Appl. Anal.*, 102(5):1312–1334, 2023. doi:10.1080/00036811.2021.1981878.

[4] L. Amir, M. Kern, J. E. Roberts, and V. Martin. Décomposition de domaine pour un milieu poreux fracturé : un modèle en 3D avec fractures. *ARIMA Rev. Afr. Rech. Inform. Math. Appl.*, Volume 5, Special Issue TAM TAM’05, 2006. doi:10.46298/arima.1851.

[5] P. Angot, F. Boyer, and F. Hubert. Asymptotic and numerical modelling of flows in fractured porous media. *ESAIM Math. Model. Numer. Anal.*, 43(2):239–275, 2009. doi:10.1051/m2an/2008052.

[6] P. F. Antonietti, C. Facciolà, and M. Verani. Polytopic discontinuous Galerkin methods for the numerical modelling of flow in porous media with networks of intersecting fractures. *Comput. Math. Appl.*, 116:116–139, 2022. doi:10.1016/j.camwa.2021.08.015.

[7] M. F. Benedetto, A. Borio, F. Kyburg, J. Mollica, and S. Scialò. An arbitrary order mixed virtual element formulation for coupled multi-dimensional flow problems. *Comput. Methods Appl. Mech. Engrg.*, 391:Paper No. 114204, 27, 2022. doi:10.1016/j.cma.2021.114204.

[8] I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gläser, E. Keilegavlen, A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, S. Burbulla, P. Devloo, O. Duran, M. Favino, J. Hennicker, I.-H. Lee, K. Lipnikov, R. Masson, K. Mosthaf, M. G. C. Nestola, C.-F. Ni, K. Nikitin, P. Schädle, D. Svyatskiy, R. Yanbarisov, and P. Zulian. Verification benchmarks for single-phase flow in three-dimensional fractured porous media. *Adv. in Water Res.*, 147:103759, 2021. doi:10.1016/j.advwatres.2020.103759.

[9] I. Berre, F. Doster, and E. Keilegavlen. Flow in Fractured Porous Media: A Review of Conceptual Models and Discretization Approaches. *Transp. Porous Media*, 130(1):215–236, 2019. doi:10.1007/s11242-018-1171-6.
[10] S. Berrone, S. Pieraccini, and S. Scialò. A PDE-Constrained Optimization Formulation for Discrete Fracture Network Flows. *SIAM J. Sci. Comput.*, 35(2):B487-B510, 2013. doi: 10.1137/120865884.

[11] S. Berrone, S. Pieraccini, and S. Scialò. Flow simulations in porous media with immersed intersecting fractures. *J. Comput. Phys.*, 345:768-791, 2017. doi:10.1016/j.jcp.2017.05.049.

[12] D. Boffi, F. Brezzi, and M. Fortin. *Mixed finite element methods and applications*, volume 44 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013. doi:10.1007/978-3-642-36519-5.

[13] W. M. Boon, J. M. Nordbotten, and I. Yotov. Robust discretization of flow in fractured porous media. *SIAM J. Numer. Anal.*, 56(4):2203-2233, 2018. doi:10.1137/17M1139102.

[14] H. Borouchaki, P. Laug, and P.-L. George. Parametric surface meshing using a combined advancing-front generalized Delaunay approach. *Internat. J. Numer. Methods Engrg.*, 49(1-2):233-259, 2000. doi:10.1002/1097-0207(20000910/20)49:1/2<233::AID-NME931>3.0.CO;2-G.

[15] K. Brenner, M. Groza, C. Guichard, G. Lebeau, and R. Masson. Gradient discretization of hybrid dimensional Darcy flows in fractured porous media. *Numerische Mathematik*, 134(3):569-609, 2016. doi:10.1007/s00211-015-0782-x.

[16] K. Brenner, J. Hennicker, R. Masson, and P. Samier. Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces. *IMA J. Numer. Anal.*, 37(3):1551-1585, 2017. doi:10.1093/imanum/drw044.

[17] K. Brenner, J. Hennicker, R. Masson, and P. Samier. Hybrid-dimensional modelling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions. *J. Comput. Phys.*, 357:100-124, 2018. doi:10.1016/j.jcp.2017.12.003.

[18] S. C. Brenner. A multigrid algorithm for the lowest-order raviart-thomas mixed triangular finite element method. *SIAM J. Numer. Anal.*, 29(3):647-678, 1992. doi:10.1137/0729042.

[19] S. Burbulla, M. Hörl, and C. Rohde. Flow in porous media with fractures of varying aperture. *SIAM J. Sci. Comput.*, 45(4):A1519-A1544, 2023. doi:10.1137/22M1510406.

[20] F. Chave, D. A. Di Pietro, and L. Formaggia. A Hybrid High-Order Method for Darcy Flows in Fractured Porous Media. *SIAM J. Sci. Comput.*, 40(2):A1063-A1094, 2018. doi: 10.1137/17M1119500.

[21] Z. Chen. Equivalence between and multigrid algorithms for nonconforming and mixed methods for second-order elliptic problems. *East-West J. Numer. Math.*, 4(1):1-33, 1996.

[22] P. G. Ciarlet. *The finite element method for elliptic problems*. Studies in Mathematics and its Applications, Vol. 4. North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.

[23] B. Cockburn and J. Gopalakrishnan. A characterization of hybridized mixed methods for second order elliptic problems. *SIAM J. Numer. Anal.*, 42(1):283-301, 2004. doi:10.1137/S0036142902417893.

[24] B. Cockburn and J. Gopalakrishnan. Error analysis of variable degree mixed methods for elliptic problems via hybridization. *Math. Comput.*, 74(252):1653-1677, 2005. doi:10.1090/S0025-5718-05-01741-2.
[25] F. Coulomb and C. Fedon-Magnaud. Mixed and mixed-hybrid elements for the diffusion equation. *Nucl. Sci. Eng.*, 100(3):218-225, 1988. doi:10.13182/NSE88-A29034.

[26] M. Del Pra, A. Fumagalli, and A. Scotti. Well Posedness of Fully Coupled Fracture/Bulk Darcy Flow with XFEM. *SIAM J. Numer. Anal.*, 55(2):785-811, 2017. doi:10.1137/15M1022574.

[27] P. R. B. Devloo, W. Teng, and C.-S. Zhang. Multiscale Hybrid-Mixed Finite Element Method for Flow Simulation in Fractured Porous Media. *Comp. Model. Eng. Sci.*, 119(1):145-163, 2019. doi:10.32604/cmes.2019.04812.

[28] V. Dolean, P. Jolivet, and F. Nataf. *An Introduction to Domain Decomposition Methods*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015. doi:10.1137/1.9781611974065.

[29] O. Durán, P. R. B. Devloo, S. M. Gomes, and F. Valentin. A multiscale hybrid method for Darcy's problems using mixed finite element local solvers. *Comput. Methods Appl. Mech. Engrg.*, 354:213-244, 2019. doi:10.1016/j.cma.2019.05.013.

[30] A. Ern, T. Gudi, I. Smears, and M. Vohralík. Equivalence of local- and global-best approximations, a simple stable local commuting projector, and optimal hp approximation estimates in H(div). *IMA J. Numer. Anal.*, 42(2):1023-1049, 2022. doi:10.1093/imanum/draa103.

[31] A. Ern and J.-L. Guermond. *Finite elements I. Approximation and interpolation*, volume 72 of *Texts Appl. Math*. Cham: Springer, 2020. doi:10.1007/978-3-030-56341-7.

[32] A. Ern and J.-L. Guermond. *Finite elements II. Galerkin approximation, elliptic and mixed PDEs*, volume 73 of *Texts Appl. Math*. Cham: Springer, 2021. doi:10.1007/978-3-030-56923-5.

[33] A. Ern, F. Hédin, G. Pichot, and N. Pignet. Hybrid high-order methods for flow simulations in extremely large discrete fracture networks. *SMAI J. Comput. Math.*, 8:375-398, 2022. doi:10.5802/smai-jcm.92.

[34] B. Flemisch, A. Fumagalli, and A. Scotti. A review of the XFEM-based approximation of flow in fractured porous media. In G. Ventura and E. Benvenuti, editors, *Advances in Discretization Methods: Discontinuities, Virtual Elements, Fictitious Domain Methods*, pages 47-76. Springer International Publishing, Cham, 2016. doi:10.1007/978-3-319-41246-7_3.

[35] N. Frih, V. Martin, J. E. Roberts, and A. Saâda. Modeling fractures as interfaces with nonmatching grids. *Comput. Geosci.*, 2012. doi:10.1007/s10596-012-9302-6.

[36] A. Fumagalli, E. Keilegavlen, and S. Scialò. Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations. *J. Comput. Phys.*, 376:694-712, 2019. doi:10.1016/j.jcp.2018.09.048.

[37] A. Fumagalli and A. Scotti. A Reduced Model for Flow and Transport in Fractured Porous Media with Non-matching Grids. In A. Cangiani, R. L. Davidchack, E. Georgoulis, A. N. Gorban, J. Levesley, and M. Tretyakov, editors, *Numerical Mathematics and Advanced Applications 2011*, pages 499-507, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg. doi:10.1007/978-3-642-33134-3_53.

[38] G. Gatica. *A simple introduction to the mixed finite element method. Theory and applications*. SpringerBriefs in Mathematics. Springer, 2014. doi:10.1007/978-3-319-03695-3.

[39] P. L. George, H. Borouchaki, F. Alauzet, P. Laug, A. Loseille, and L. Maréchal. *Meshing, geometric modeling and numerical simulation. 2*. Numerical Methods in Engineering Series. John Wiley & Sons, Inc., Hoboken, NJ, 2019. Metrics, meshes and mesh adaptation, Geometric modeling and applications set. Vol. 2. doi:10.1002/9781119384380.
[40] P. Grisvard. *Elliptic problems in nonsmooth domains*, volume 69 of *Classics in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011. Reprint of the 1985 original [MR0775683], With a foreword by Susanne C. Brenner. doi:10.1137/1.9781611972030.ch1.

[41] F. Hellman, A. Målqvist, and M. Mosquera. Well-posedness and finite element approximation of mixed dimensional partial differential equations. *BIT Numerical Mathematics*, 64(1):2, 2023. doi:10.1007/s10543-023-01001-w.

[42] M. Hörl and C. Rohde. Rigorous derivation of discrete fracture models for Darcy flow in the limit of vanishing aperture. *Netw. Heterog. Media*, 19(1):114–156, 2024. doi:10.3934/nhm.2024006.

[43] H. Hoteit, J. Erhel, R. Mosé, B. Philippe, and P. Ackerer. Numerical Reliability for Mixed Methods Applied to Flow Problems in Porous Media. *Comput. Geosci.*, 6:161–194, 2002. doi:10.1023/A:1019988901420.

[44] H. Hoteit and A. Firoozabadi. An efficient numerical model for incompressible two-phase flow in fractured media. *Adv. in Water Res.*, 31(6):891–905, 2008. doi:10.1016/j.advwatres.2008.02.004.

[45] J. D. Hyman, M. R. Sweeney, C. W. Gable, D. Svyatsky, K. Lipnikov, and J. D. Moulton. Flow and transport in three-dimensional discrete fracture matrix models using mimetic finite difference on a conforming multi-dimensional mesh. *J. Comput. Phys.*, 466:111396, 2022. doi:10.1016/j.jcp.2022.111396.

[46] P. Jolivet, M. Kern, F. Nataf, G. Pichot, and D. Z. Vasquez. Domain decomposition preconditioners for efficient parallel simulations of single-phase flow in three-dimensional fractured porous media with a very large number of fractures. Preprint, 2025. URL: https://inria.hal.science/hal-05029676.

[47] M. Karimi-Fard, L. J. Durlofsky, and K. Aziz. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. *SPE Journal*, 9(02):227–236, 06 2004. doi:10.2118/88812-PA.

[48] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, and I. Berre. Porepy: an open-source software for simulation of multiphysics processes in fractured porous media. *Comput. Geosci.*, 25:243–265, 2021. doi:10.1007/s10596-020-10002-5.

[49] M. Kern, G. Pichot, and D. Z. Vasquez. Performance of algebraic preconditioners for large-scale simulations of single-phase flow in three-dimensional fractured porous media. Preprint, 2025. URL: https://inria.hal.science/hal-05029652.

[50] K. Kumar, F. List, I. S. Pop, and F. A. Radu. Formal upscaling and numerical validation of unsaturated flow models in fractured porous media. *J. Comput. Phys.*, 407:109138, 21, 2020. doi:10.1016/j.jcp.2019.109138.

[51] P. Laug and G. Pichot. Mesh Generation and Flow Simulation in Large Tridimensional Fracture Networks. In S. Carillo, C. Conti, D. Mansutti, F. Pitolli, and R. M. Spitaleri, editors, *IMACS Series in Computational and Applied Mathematics*, volume 22, pages 71–80, 2019. URL: https://hal.archives-ouvertes.fr/hal-02102811.

[52] H. Leng and H. Chen. An adaptive hybridizable discontinuous Galerkin method for Darcy-Forchheimer flow in fractured porous media. *Math. Models Methods Appl. Sci.*, 34(13):2497–2536, 2024. doi:10.1142/S0218202524500532.
[53] H. Leng and H. Chen. Adaptive interior penalty hybridized discontinuous Galerkin methods for Darcy flow in fractured porous media. *IMA J. Numer. Anal.*, 44(4):2165-2197, 2024. doi:10.1093/imanum/drad054.

[54] F. List, K. Kumar, I. S. Pop, and F. A. Radu. Rigorous upscaling of unsaturated flow in fractured porous media. *SIAM J. Math. Anal.*, 52(1):239-276, 2020. doi:10.1137/18M1203754.

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

[56] J. M. Nordbotten, W. M. Boon, A. Fumagalli, and E. Keilegavlen. Unified approach to discretization of flow in fractured porous media. *Comput. Geosci.*, 23(2):225-237, 2019. doi: 10.1007/s10596-018-9778-9.

[57] B. Nœtinger. A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow. *J. Comput. Phys.*, 283:205-223, 2015. doi:10.1016/j.jcp.2014.11.038.

[58] P. F. Antonietti, L. Formaggia, A. Scotti, M. Verani, and N. Verzott. Mimetic finite difference approximation of flows in fractured porous media. *ESAIM Math. Model. Numer. Anal.*, 50(3):809-832, 2016. doi:10.1051/m2an/2015087.

[59] G. Pichot, J. Erhel, and J. R. de Dreuzy. A mixed hybrid mortar method for solving flow in discrete fracture networks. *Appl. Anal.*, 89(10):1629-1643, 2010. doi:10.1080/00036811.2010.495333.

[60] G. Pichot, J. Erhel, and J.-R. de Dreuzy. A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks. *SIAM J. Sci. Comput.*, 34(1):B86-B105, 2012. doi:10.1137/100804383.

[61] V. Reichenberger, H. Jakobs, P. Bastian, and R. Helmig. A mixed-dimensional finite volume method for two-phase flow in fractured porous media. *Adv. in Water Res.*, 29(7):1020-1036, 2006. doi:10.1016/j.advwatres.2005.09.001.

[62] J. E. Roberts and J.-M. Thomas. Mixed and hybrid methods. In *Handbook of numerical analysis, Vol. II*, Handb. Numer. Anal., II, pages 523-639. North-Holland, Amsterdam, 1991.

[63] N. Spillane, V. Dolean, P. Hauret, F. Nataf, C. Pechstein, and R. Scheichl. Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps. *Numerische Mathematik*, 126(4):741-770, 2014. doi:10.1007/s00211-013-0576-y.

[64] L. Sun, X. Tang, M. Li, A. Abdelaziz, K. Aboayanah, Q. Liu, and G. Grasselli. Flow simulation in 3D fractured porous medium using a generalized pipe-based cell-centered finite volume model with local grid refinement. *Geomech. Energy Environ.*, 36:100505, 2023. doi:10.1016/j.gete.2023.100505.

[65] A. Veeser. Best error localizations for piecewise polynomial approximation of gradients, functions and functionals. In *Numerical mathematics and advanced applications—ENUMATH 2017*, volume 126 of *Lect. Notes Comput. Sci. Eng.*, pages 357-365. Springer, Cham, 2019. doi:10.1007/978-3-319-96415-7_31.

[66] M. Vohralík, J. Maryška, and O. Severýn. Mixed and nonconforming finite element methods on a system of polygons. *Appl. Numer. Math.*, 57(2):176-193, 2007. doi:10.1016/j.apnum.2006.02.005.

[67] L. Zhao and E. Chung. An adaptive discontinuous Galerkin method for the Darcy system in fractured porous media. *Comput. Geosci.*, 26(6):1581-1596, 2022. doi:10.1007/s10596-022-10171-5.