# Performance of algebraic preconditioners for large-scale simulations 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. Performance of algebraic preconditioners for large-scale simulations of single-phase flow in three-dimensional fractured porous media. 2025. [hal-05029652]

HAL Id: hal-05029652

https://inria.hal.science/hal-05029652v1

Preprint submitted on 11 Apr 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
# Performance of algebraic preconditioners for large-scale simulations of single-phase flow in three-dimensional fractured porous media

Michel Kern*†

Géraldine Pichot*†

Daniel Zegarra Vasquez*†

[April 11, 2025 ; first draft, long version]

## Abstract

A three dimensional fractured-porous medium is classically described by a reduced model, called the discrete fracture matrix (DFM), where the rock matrix remains three-dimensional, while the fracture network is of co-dimension one. Solving single-phase flow in such a medium is particularly challenging, especially when a large number of fractures is present. This difficulty primarily arises from the intricate complexity introduced by the fracture network. This work demonstrates that solving flow in such systems is feasible, even for discrete fracture networks as large as 90,000 fractures. The mixed-hybrid finite element method is employed to discretize the flow problem, resulting in a large linear system that must be solved. However, the efficiency of direct linear solvers is severely hindered by the three-dimensional nature of the problem, requiring a substantial amount of RAM. Iterative methods, when preconditioned with algebraic multigrid techniques, offer a better alternative in terms of memory usage. However, their effectiveness declines significantly as the fracture network grows and transmissivity within the network becomes highly heterogeneous, as both lead to poor conditioning of the linear system. Even more critically, we show that achieving an adequately low tolerance for iterative solvers is essential for accurate flow simulations. The choice of tolerance strongly influences the quality of the solution, insufficient precision can prevent local mass conservation in some cases, resulting in non-physical values in fully advective transport simulations.

## Keywords

large-scale simulations, single-phase flow, heterogeneous fractured porous media, discrete fracture-matrix model, conforming meshes, mixed-hybrid finite element method, conjugate gradient, algebraic multigrid preconditioners, transport.

## Introduction

We are interested in solving single-phase flow in the subsurface. Imaging the subsurface is a challenge as measurements of boreholes, outcrops, lineament maps, etc., provide data only at some specific scales. Moreover, the measured data carries a lot of uncertainties [12, 71]. It is observed that fractures are ubiquitous in the subsurface. They are discontinuities in the rock matrix medium in the form of narrow zones. There exists many relevant geometric representations of fractures. On the one hand, we have implicit representation of fractures via single-continuum models [76, 82] and multi-continuum models [80, 2]. On the other hand, we have explicit representation of fractures thanks to Discrete Fracture Networks (DFN) [28, 34, 33] when the rock matrix is impervious, and Discrete Fracture Matrix (DFM), when the rock matrix and the fracture network have some interaction [20, 14, 25]. This work focuses specifically on the DFM representation, which is the closest to the original fractured porous medium in terms of accuracy of representation of fractures and accuracy of representation of the background medium. One particularity of the DFM representation is that fractures are modeled as lower dimensional objects [3, 72, 8]. The rock matrix is a three-dimensional domain and fractures, which are conductive, are considered as two-dimensional flat surfaces inside a three-dimensional space. To be more precise, fractures are modeled as ellipses, disks or polygons. The distribution of the fractures either comes from field observations [47] or follows some probability law, which is the core of genetic algorithms [33].

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

<footnote>†CERMICS, École Nationale des Ponts et Chaussées, 6 avenue Blaise Pascal, 77420 Champs-sur-Marne, France.
E-mails: michel.kern@inria.fr, geraldine.pichot@inria.fr, daniel.zegarra-vasquez@inria.fr.</footnote>
Fractures play a major role in many different applications, such as oil and gas extraction [22, 89], geothermal energy production [68, 81], nuclear waste storage [79, 83], CO2 sequestration [57, 30], water extraction from aquifers [42, 86], and cannot be neglected. Fractures can be conductive [1] or acting as barriers [72, 8].

These applications in fractured porous media involve various physical and chemical phenomena. For example, there could be flow, which can be single-phase [14, 40], two-phase [18] or multi-phase [91]. There could also be transport [56, 11], chemical reactions [49, 46], and/or mechanical processes (e.g. fracture growth [78]). In this work we focus on single-phase flow and purely advective transport in conductive fracture networks.

With a DFM representation, 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 [72, 75]. Additionally, the hydraulic head can either exhibit discontinuities between the two sides of a fracture or remain continuous. Here, we assume that the hydraulic head is continuous across fractures, implying that fractures do not act as barriers and do not obstruct flow. As a result, we work with a fully linear elliptic model for steady-state single-phase flow, incorporating a transfer term between the rock matrix and the fractures.

For meshing the domain, we opt for a conforming mesh approach [59, 63, 73] rather than a non-conforming one [69, 48, 15, 16, 45]. In the case of a conforming mesh, the surface mesh of the fracture network consists of a set of faces that act as constraints for the volume mesh. In this work, we consider simplicial meshes. The fractures often intersect at very small angles, making the meshing of the volume domain relative to the surface fracture networks challenging. Typically, the surface mesh must maintain a good aspect ratio for all triangles to ensure that the volume mesh can be generated [63]. One might consider simplifying the geometry [58]. In this work, we retain the original geometry to preserve the number of fractures. To ensure a minimum quality of the surface mesh while minimizing alterations to the original geometry, we apply automatic corrections, which are implemented in the software MODFRAC [21, 67]. Then a robust implementation, such as GHS3D [51], is used to generate the 3D volume mesh. This approach enables us to explore large 3D DFM models, potentially containing a large number of fractures. While it is possible to generate geometries with many fractures, meshing such complex geometries, where intersections and small angles are prominent, has not been widely explored in the recent literature.

There exists many numerical methods for flow in 3D DFM [14]. For instance, there is the mimetic finite difference method [70, 77], the Vertex Approximate Gradient scheme [23, 24], the Hybrid Finite Volumes [23, 24], the multi-point flux approximation [74, 64, 66], virtual element method [74, 64], mixed finite element method [74, 64, 20], the multi-scale mixed-hybrid finite element method [39, 41], a Hybrid High-Order Method [29], and many others.

In this work, we adopt a mixed-hybrid finite element method (MHFEM) of the lowest order, combined with conforming meshes for the DFM. The unknowns in this approach are the velocity, the hydraulic head computed at the center of each tetrahedron for the rock matrix (or at the center of each triangle for the fracture network), and the trace of the hydraulic head computed at the center of each tetrahedron face (or at the center of each triangle edge for the fracture network) [27, 50, 17]. This approach originated in [55], where it was applied to incompressible two-phase flow. The advantages of this method are as follows. First, due to the use of conforming simplicial meshes, the unknowns at the matrix/fracture interfaces are easily coupled, particularly in the case of conductive fractures where the hydraulic head remains continuous across the fractures. Second, the resulting linear system is square, sparse, symmetric, and positive definite (SPD), making it amenable to solution by both direct and iterative methods. However, there are some drawbacks. The resulting linear system can be very large and ill-conditioned, a consequence of the intricate 2D-3D mesh due to the complex fracture network geometry. This can lead to nearly flat tetrahedral elements and/or triangles in the mesh. Ultimately, these issues negatively affect the quality of the velocity field, which is computed in a post-processing phase from the hydraulic heads.

Here, since the linear system is SPD, it can be solved using black-box solvers, including direct
solvers such as UMFPACK [32] and MUMPS [5, 4], as well as iterative algorithms like the conjugate gradient (CG) method. Of course, CG requires an efficient preconditioner. In this work, we investigate black-box preconditioners such as Incomplete Cholesky Factorization or algebraic multigrid methods (AMG) [9, 92, 36]. Investigation of more advanced preconditioners is the topic of a companion paper [61]. Assessing solver performances for an extremely large number of fractures in 2D DFN has already been studied in [44] with the Hybrid High Order (HHO) method, but never for 3D DFM. To the best of our knowledge, most 3D DFM flow simulations in recent studies handle up to thousands of fractures [13, 59, 87].

The contributions of this paper are

* to validate a workflow based on the MHFEM to simulate time-independent single-phase flow in fractured porous media;

* to extend the simulation capabilities to handle tens of thousands of fractures;

* to assess the performance of different black-box solvers, both direct and iterative (CG with different preconditioners such as Incomplete Cholesky Factorization or AMG) on these large-scale simulations;

* to study the impact of the mesh quality/heterogeneity contrasts on the solver convergence;

* to investigate how the chosen iterative solver tolerance affects the quality of the computed velocity by running purely advective transport simulations.

The outline is the following. Section 1 presents the flow problem and recalls the MHFEM and the linear system to be solved. Section 2 validates of our flow and transport codes on a benchmark from [13]. Sections 3 and 4 are the main contributions of the paper, presenting test cases of increasing complexity and the results of various solvers on these cases. In particular, the impact of the geometry of the fracture network, of the mesh quality, of the heterogeneity contrasts and of the impact of the tolerance of the iterative solver on the quality of the fluid flow solution are investigated in details. The last section provides conclusions.

# 1 Mathematical setting of the flow problem

## 1.1 Discrete fracture matrix (DFM) notation

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^f_\ell$. The rock matrix domain is defined as $\Omega^m = \Omega \setminus \bar{\Omega}^f$.

For the purpose of the mathematical description, we divide each polygon $\Omega^f_\ell$, $\ell = 1, \dots, N^f$, into subpolygons $\gamma_\ell$ along each intersection line, such that $\Omega^f = \cup_{\ell \in \mathcal{I}} \gamma_\ell$, where $\mathcal{I}$ is the index set of the subpolygons. The subpolygons are built so as to satisfy: $\gamma_\ell \cap \gamma_k = \emptyset$ if $\ell \neq k$ and that $\bar{\gamma}_\ell \cap \bar{\gamma}_k = \emptyset$ if $\ell \neq k$, is either an edge or a point or an empty set, as described in [88, 65]. The subpolygons $\gamma_\ell$, $\ell \in \mathcal{I}$, further partition the rock matrix domain $\Omega^m$ into connected components $\omega_\alpha$, such that $\Omega^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: $\omega_\alpha \cap \omega_\beta = \emptyset$ if $\alpha \neq \beta$ and that $\bar{\omega}_\alpha \cap \bar{\omega}_\beta$, $\alpha \neq \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 subpolygons on $\partial \omega_\alpha$. We further assume that both the subpolygons and the subdomains are Lipschitz. This is always achievable in practice: for example, after generating a mesh for both $\Omega^f$ and $\Omega^m$, the subpolygons and the subdomains may be built by agglomeration of mesh elements.

Dirichlet and Neumann boundary conditions are set on the porous domain on $\Gamma^D$ and $\Gamma^N$ respectively, with $\partial\Omega^m = \Gamma_0$ with $\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}}$ gathers the faces introduced by the
subdivision process (they will only play a role in Section 1.2). We call such a face a “fictitious”
face, and note that it belongs to two subdomains.

We assume that Dirichlet and Neumann boundary conditions also hold for the fracture networks
on the domain sides: $\partial\Omega^f \cap \partial\Omega = \Sigma^D \cup \Sigma^N$. For each $\ell \in \mathcal{I}$, $\Sigma^D_\ell = \partial\gamma_\ell \cap \Sigma^D$ and $\Sigma^N_\ell = \partial\gamma_\ell \cap \Sigma^N$.
In this work, homogeneous Neumann boundary conditions are applied at the fracture tips (border
of the subpolygon $\gamma_\ell$, $\ell \in \mathcal{I}$, that are not intersection and not on the border of the domain). The
fracture tips are part of $\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}$. For each $\ell \in \mathcal{I}$, we let $\Sigma_{\ell}^{\text{inter}}$ be the set of intersections of the sub-polygon $\gamma_{\ell}$, and we define $\Sigma_{\ell} = \Sigma_{\ell}^{D} \cup \Sigma_{\ell}^{N} \cup \Sigma_{\ell}^{\text{inter}}$. Let $\Sigma$ be the set of intersections between the sub-polygons. The intersections may be segments shared by two or more sub-polygons. Let $\mathcal{I}_{\sigma} = \{\ell \in \mathcal{I}, \sigma \subset \Sigma\}$ the index set of polygons that contain $\sigma$ as intersection.

1.2 Physical governing equations

We consider a single-phase, incompressible, and time-independent flow problem, as detailed for
instance in [72, 23, 24, 13, 14, 7, 16]. The unknowns are:

• the fluid flow velocity in the rock matrix domain $u^m = (u_\alpha^m)_{\alpha \in \mathcal{C}}$, obtained from the fluid velocity $u_\alpha^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 subpolygon $\gamma_\ell$, in m²/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 $\tilde{K}^f = (\tilde{K}_\ell^f)_{\ell \in \mathcal{I}}$. Their physical unit is m/s. Furthermore, in the fracture network, the aperture is a scalar denoted $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²/s. In the sequel, we only work with $K^f$ (no closing of fractures).

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)).

Let $q_{\alpha}^{m,N}$ (resp. $q_{\ell}^{f,N}$) be the Neumann boundary conditions applying on $\Gamma_{\alpha}^{N}$ (resp. $\Sigma_{\ell} \cap \Sigma^{N}$) for all $\alpha \in \mathcal{C}$ (resp. for all $\ell \in \mathcal{I}$). The Neumann condition in the rock matrix (resp. in the fracture network) is denoted $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²/s)).

Let $g_{\alpha}^{m,D}$ (resp. $g_{\ell}^{f,D}$) be the Dirichlet boundary conditions applying on $\Gamma_{\alpha}^{D}$ (resp. $\Sigma_{\ell} \cap \Sigma^{D}$), for all $\alpha \in \mathcal{C}$ (resp. for all $\ell \in \mathcal{I}$). 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 p^m$ (resp. $\gamma^m 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 u^m$ denote the normal trace of $u^m$ along (a part of) the boundary of $\Omega^m$. Informally, $\gamma_n^m u^m = u^m \cdot \nu|_{\partial\Omega^m}$, where $\nu$ 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 write $\gamma_n^{m,\perp} u_\alpha^m$. If the part of the boundary is an immersed polygon $\gamma_\ell$, $\ell \in \mathcal{I}$, we write $\gamma_n^{m,\perp} u_{\alpha_\ell \pm}$ depending on which side of the polygon $\omega_{\alpha_\ell \pm}$ we are looking at.

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

For a face $\gamma$ shared by two subdomains $\alpha^+$ and $\alpha^-$, we define the jump across $\gamma$ as

$$ [[u^m]]|_{\gamma} = \gamma_n^{m,+} u_{\alpha^+}^m + \gamma^{m,-} u_{\alpha^-}^m. \qquad (1.1) $$

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

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

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

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

$$ \gamma_n^{\mathfrak{m}} u^{\mathfrak{m}} = q_{\alpha}^{\mathfrak{m},N} \quad \text{on } \Gamma_{\alpha}^{N}. \qquad (1.2d) $$

Additionally, across a fictitious face shared by two subdomains both the hydraulic head and the jump of the flux is zero.

In each polygon $\gamma_\ell$, $\ell \in \mathcal{I}$, the flow problem is written as follows:

$$ \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}} = -K_{\ell}^{\mathfrak{f}} \nabla_{\tau} p_{\ell}^{\mathfrak{f}} \quad \text{in } \gamma_{\ell}, \qquad (1.3b) $$

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

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

where $[[u^m]]|_{\gamma_\ell}$ is an exchange term that couples the flow between the porous subdomains $\omega_{\alpha_\ell^+}$ and $\omega_{\alpha_\ell^-}$ and the fracture $\gamma_\ell = \bar{\omega}_{\alpha_\ell^+} \cap \bar{\omega}_{\alpha_\ell^-}$. This term acts as a source term for the fracture flow problem and is defined as follows:

$$ [[u^m]]|_{\gamma_\ell} = \gamma_n^{m,+} u_{\alpha_\ell^+}^m + \gamma^{m,-} u_{\alpha_\ell^-}^m, \quad \forall \ell \in \mathcal{I}. \qquad (1.4) $$

We assume the continuity of the hydraulic heads between the subdomains and the polygons [8]:

$$ \gamma^{\mathfrak{m}} p_{\alpha_\ell^+}^{\mathfrak{m}} = \gamma^{\mathfrak{m}} p_{\alpha_\ell^-}^{\mathfrak{m}} = p_{\ell}^{\mathfrak{f}}, \quad \forall \ell \in \mathcal{I}. \qquad (1.5) $$

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

$$ \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}, \forall \sigma \subset \Sigma, \qquad (1.6a) $$

$$ \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.6b) $$

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

We assume that $K^{\mathfrak{m}}(x)$ (resp. $K^{\mathfrak{f}}(x)$) is a uniformly symmetric and positive definite tensor so 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 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}}, \qquad (1.7a) $$

$$ K_{\min}^{\mathfrak{f}} \|\xi\|_2^2 \le \xi^T K^{\mathfrak{f}}(x) \xi \le K_{\max}^{\mathfrak{f}} \|\xi\|_2^2, \quad \forall \xi \in \mathbb{R}^2, \forall x \in \Omega^{\mathfrak{f}}. \qquad (1.7b) $$

We further assume that the source terms $f_\alpha^m$ (resp. $f_\ell^\mathfrak{f}$) live in $L^2(\Omega_\alpha^m)$ (resp. $L^2(\gamma_\ell^\mathfrak{f})$), that the Dirichlet boundary conditions functions $g_\alpha^{m,D}$ (resp. $g_\ell^{\mathfrak{f},D}$) are in $H^{\frac{1}{2}}(\Gamma_\alpha^D)$ (resp. $H^{\frac{1}{2}}(\Sigma_\ell^D)$) and that the Neumann boundary conditions functions $q_\alpha^{m,N}$ (resp. $q_\alpha^{\mathfrak{f},N}$) are in $H^{-\frac{1}{2}}(\Gamma_\alpha^N)$ (resp. $H^{-\frac{1}{2}}(\Sigma_\ell^N)$).
## 1.3 Mixed-hybrid finite element method

We recall that the rock matrix is three-dimensional while the fracture-network is two-dimensional in a 3D space. Thus, the governing physical equations (1.2) and (1.3) are approximated almost independently from each other, but not totally since there is a coupling term (1.4) from the rock matrix-fracture interaction acting as a source term in (1.3a). It accounts for the volumetric transfer between the rock matrix and the fracture network. We apply the mixed-hybrid finite element discretization to solve this flow problem. Here we recall only the main notation, the unknowns and the linear system. More details on the mathematical and numerical analysis can be found in [65].

### 1.3.1 Discretization

We first describe the fracture mesh. Let $K_h$ be a mesh of the fracture network, composed of $\mathcal{N}(K_h)$ triangles, and let $\mathcal{E}_h$ be the set of edges of $K_h$, with the number of edges denoted by $\mathcal{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^D$ the set of Dirichlet edges and $\mathcal{E}_h^N$ the set of Neumann edges. We have $\mathcal{E}_h = \mathcal{E}_h^{\text{in}} \cup \mathcal{E}_h^{\text{inter}} \cup \mathcal{E}_h^D \cup \mathcal{E}_h^N$. The number of elements in each of these sets is denoted respectively by $\mathcal{N}(\mathcal{E}_h^{\text{in}})$, $\mathcal{N}(\mathcal{E}_h^{\text{inter}})$, $\mathcal{N}(\mathcal{E}_h^D)$ and $\mathcal{N}(\mathcal{E}_h^N)$. Let $\mathcal{E}_h \setminus \mathcal{E}_h^D$ the set of non-Dirichlet edges of cardinal $\mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)$.

In the same way, let $\mathcal{T}_h$ be a mesh of the porous domain composed of $\mathcal{N}(\mathcal{T}_h)$ tetrahedral elements. Let $\mathcal{F}_h$ be the set of faces of $\mathcal{T}_h$, whose number is denoted by $\mathcal{N}(\mathcal{F}_h)$. We define $\mathcal{F}_h^{\text{in}}$ the set of inner faces, not on the boundary and not fractures (rock faces), $\mathcal{F}_h^{\text{frac}}$ the set of tetrahedral faces that are fracture faces (triangles, see below), $\mathcal{F}_h^D$ the set of tetrahedral faces that are Dirichlet faces and $\mathcal{F}_h^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^D \cup \mathcal{F}_h^N$. The number of inner faces is denoted $\mathcal{N}(\mathcal{F}_h^{\text{in}})$, the number of fracture faces is denoted $\mathcal{N}(\mathcal{F}_h^{\text{frac}})$ and the number of Dirichlet and Neumann faces $\mathcal{N}(\mathcal{F}_h^D)$ and $\mathcal{N}(\mathcal{F}_h^N)$ respectively. Let $\mathcal{F}_h \setminus \mathcal{F}_h^D$ the set of non-Dirichlet faces of cardinal $\mathcal{N}(\mathcal{F}_h \setminus \mathcal{F}_h^D)$. Let $\mathcal{F}_h^{\text{in},N}$ the set of non-fracture and non-Dirichlet faces of cardinal $\mathcal{N}(\mathcal{F}_h^{\text{in},N})$.

We assume that $K_h$ and $\mathcal{T}_h$ both satisfy the usual properties of a (simplicial) finite element mesh (cf. for instance [43, Sec. 8.2]) and are built so as to accommodate the boundary conditions.

The degrees of freedom are the fluxes on each face of each tetrahedron (resp. each edge of each triangle) for the matrix and the fractures respectively:

$$Q^{\text{m}} = (q_{T,F}^{\text{m}})_{(T \in \mathcal{T}_h, F \in \{1,2,3,4\})} \in \mathbb{R}^{4\mathcal{N}(\mathcal{T}_h)}, \quad Q^{\text{f}} = (q_{K,E}^{\text{f}})_{(K \in \mathcal{K}_h, E \in \{1,2,3\})} \in \mathbb{R}^{3\mathcal{N}(\mathcal{K}_h)};$$

The spaces for the pressure are piecewise constant functions over each tetrahedra or triangles respectively:

$$P^{\text{m}} = (p_T^{\text{m}})_{T \in \mathcal{T}_h} \in \mathbb{R}^{\mathcal{N}(\mathcal{T}_h)}, \quad P^{\text{f}} = (p_K^{\text{f}})_{K \in \mathcal{K}_h} \in \mathbb{R}^{\mathcal{N}(\mathcal{K}_h)};$$

The spaces for the multipliers are piecewise constant on the inner faces of the tetrahedral mesh, and on the edges of the triangular mesh respectively:

$$\Lambda^{\text{m},\text{m}} = (\lambda_F^{\text{m}})_{F \in \mathcal{N}(\mathcal{F}_h^{\text{in},N})} \in \mathbb{R}^{\mathcal{N}(\mathcal{F}_h^{\text{in},N})}, \quad \Lambda^{\text{f}} = (\lambda_E^{\text{f}})_{E \in \mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)} \in \mathbb{R}^{\mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)}.$$

An illustration with two tetrahedra elements $T$ and $T'$ (rock matrix) and one triangle $K$ (fracture network) is displayed on Figure 1).
<figure><img src="image_2.png" /><figcaption>Figure 1: Mixed-hybrid finite element method unknowns</figcaption></figure>

### 1.3.2 Linear system

After discretization with the mixed-hybrid finite element method and elimination of the flux unknowns, one first obtains a system with $P^m$, $\Lambda^{m,m}$, $P^f$ and $\Lambda^f$ as the unknowns. Algebraically, we can eliminate the mean hydraulic heads in the rock matrix $P^m$. Then we obtain a Schur complement matrix for the traces of hydraulic head $\Lambda^{m,m}$ and $\Lambda^f$ and the hydraulic head $P^f$ unknowns. The resulting linear system is square, sparse, symmetric, positive and definite [54]

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

where $A^{m,m}$, $A^{m,f}$, $A^{f,f}$, $\mathcal{F}^{m,m}$ and $\mathcal{F}^{m,f}$ are matrices and vectors defined in Appendix A.2.

The vector $\mathcal{V}^{m,m}$ (resp. $\mathcal{V}^{m,f}$) corresponds to Dirichlet and Neumann boundary conditions on matrix-matrix triangle interfaces (resp. matrix-fracture triangle interfaces). The vector $\mathcal{V}^f$ corresponds to Dirichlet and Neumann boundary conditions on fracture edges.

Once the linear system (1.8) is solved, the rock matrix mean hydraulic heads $P^m$ are easily by inverting the diagonal matrix $D^m$:

$$ D^m \mathbf{P}^m = R^{m,m} \Lambda^{m,m} + R^{m,f} \mathbf{P}^f + \mathcal{F}^m. \quad (1.9) $$

Finally the vectors $Q^m$ (resp. $Q^f$) containing the fluxes $q_{T,F}^m$, on each tetrahedral element $T$ through its face $F$ (resp. $q_{K,E}^m$ on each triangle $K$ through its edge $E$) are derived locally, independently from one element to another.

In the rock matrix,

$$ q_{T,F}^{m} = \left( \sum_{F' \subset \partial T} ((B_{T}^{m})^{-1})_{F,F'} \right) p_{T}^{m} - \sum_{F' \subset \partial T} ((B_{T}^{m})^{-1})_{F,F'} \lambda_{F'}^{m}. \quad (1.10) $$

where $B_T^m$ is the local finite element matrix of the tetrahedral element $T$, which is a 4×4 symmetric, positive and definite matrix and whose inverse is denoted $(B_T^m)^{-1}$ and we have extended the definition of the traces of hydraulic head at the fracture/matrix interface (case of a conductive fracture) by using the continuity of hydraulic heads:

$$ \lambda_F^m = p_K^f, \quad \text{if } F \in \mathcal{F}_h^{\text{frac}} \text{ (which coincides geometrically with } K). \quad (1.11) $$

In the fracture,

$$ q_{K,E}^{f} = \left( \sum_{E' \subset \partial K} ((B_{K}^{f})^{-1})_{E,E'} \right) p_{K}^{f} - \sum_{E' \subset \partial K} ((B_{K}^{f})^{-1})_{E,E'} \lambda_{E'}^{f}. \quad (1.12) $$
where $B_K^\ell$ the local finite element matrix of the triangle $K$, which is a $3 \times 3$ symmetric, positive and definite matrix and whose inverse is denoted $(B_K^\ell)^{-1}$.

The velocities in the rock matrix (resp. in the fractures) can then be computed locally using
the Raviart-Thomas(-Nédélec) basis functions [65] (see Appendix A.1).

From this point onward, whenever we refer to “the linear system”, we mean the linear system (1.8).

## 2 Software implementation and validation

In this Section, we detail the workflow and software we are using to perform the mesh generation and flow and transport simulations. As described below, we are using the surface mesh generator **MODFRAC**<fnref n="1" /> for the DFN. The **.mesh** file generated by **MODFRAC** serves as input for the volume mesher **GHS3D**<fnref n="2" /> [51]. We have implemented our own software to solve flow and transport in DFM, called **nef-flow-fpm** and **nef-transport-fpm** respectively. This Section also presents the validation of our workflow on a test case from [13].

### 2.1 Mesh generation

The meshing of the fractured porous medium, modeled as a Discrete Fracture Matrix, is performed in two stages.

<figure><img src="image_3.png" /><figcaption>Figure 2: (left) Example of a surface mesh performed at Stage 1; (right) Volume mesh performed at Stage 2.</figcaption></figure>

Stage 1. The cube faces and the fractures are meshed with triangles. Notice the mesh is conforming to the intersections between fractures. See Figure 2 (left). The mesh of the fractured network is performed with the software **MODFRAC** [21, 67]. The outputs generated by **MODFRAC** follow the description provided in [35]. To the best of our knowledge, **MODFRAC** is the only mesh generator capable of handling large scale fracture networks (in this paper, up to 90k fractures), which enables us to consider very large linear systems (up to 38.9M dofs) arising from the discretization.

Stage 2. Based on the surface mesh obtained from stage 1, the volume mesh is generated with the robust and reliable software **GHS3D** [51]. See Figure 2 (right). In this work, **GHS3D** is used as a blackbox.

We assess the quality of the generated mesh through two formulas:

<footnote n="1">¹**MODFRAC** is a proprietary software owned by Inria project-teams GAMMA3 and SERENA and the University of Technology of Troyes.</footnote>

<footnote n="2">²**GHS3D** is a proprietary software owned by Inria project-team GAMMA3.</footnote>
* For a triangle $K \in \mathcal{K}_h$, the chosen formula to compute the triangle quality $\delta(K)$ is the following:

$$ \delta(K) = 4\sqrt{3} \frac{|K|}{\sum_{i=1}^{3} l_i^2}, \qquad (2.1) $$

with $|K|$ the area of the triangle $K$ and $l_i$ the size of the edge $i$ of the triangle $K$, $i = 1, \dots, 3$.

* For a tetrahedron $T \in \mathcal{T}_h$, the chosen formula to compute the quality $\delta(T)$ is the following:

$$ \delta(T) = 72\sqrt{3} \frac{|T|}{(\sum_{i=1}^{6} l_i^2)^{3/2}}, \qquad (2.2) $$

with $|T|$ the volume of the tetrahedron $T$ and $l_i$ the size of the edge $i$ of the element $T$, $i = 1, \dots, 6$.

If a triangle is regular (equilateral) then, $\delta(K) = 1$. If a tetrahedron is regular, then $\delta(T) = 1$. In both cases, a value of $\delta$ close to zero translates into an element of very poor quality.

In MODFRAC, a user quality threshold can be imposed as input so that the quality of the 2D surface mesh is not too low. In this work, for all considered meshes, the quality threshold for MODFRAC is chosen equal to $10^{-4}$, and we do not impose a quality threshold in GHS3D. Also, the user can choose a gradation parameter in MODFRAC. It indicates the maximum ratio of lengths between any two adjacent edges. In this work, unless indicated otherwise, gradation is equal to 2.0.

All the 2D surface meshes are generated with MODFRAC and are always generated on a Laptop Intel Core i7 4 cores CPU 32GiB RAM, hereafter referred to as [PC], with 4 POSIX-threads. All volume meshes are generated using GHS3D in sequential mode on a Mac cluster with a 16-core Intel(R) Xeon(R) W-3245 CPU @ 3.20GHz and 64 GiB of RAM, hereafter referred to as [MAC]

## 2.2 Workflow for flow and transport

All the flow simulations presented in this paper are performed with our MATLAB code nef-flow-fpm, which implements the mixed-hybrid finite element method described in the previous section. It massively uses MATLAB vectorization to speed up the computations. The code nef-flow-fpm is structured as follows. First, it loads and pre-processes the fracture network and the rock matrix corresponding mesh files, computes the local matrices and assembles the global linear system. Second, it calls a direct or an iterative solver (with a preconditioner if requested), or alternatively stores the linear system in PETSc binary format in order to solve it externally with PETSc built-in solvers. Third, it post-processes the solution to compute the hydraulic head in the rock matrix, the fluxes and the velocity fields. Finally, it runs a series of tests to check whether mass conservation, flux continuity, and boundary conditions are satisfied.

To further validate the computed velocity fields, we propose to run purely advective transport simulations. These simulations are run with our Eulerian MATLAB code nef-transport-fpm. It implements a finite volume discretization as spatial discretization and a finite difference scheme with constant time step as temporal discretization method (backward Euler scheme). In the following simulations, we use the fluxes computed with nef-flow-fpm as inputs for nef-transport-fpm.

## 2.3 Matrix scaling

Due to the potentially heterogeneous properties of fractured porous media, the matrices computed have elements that can vary over several orders of magnitude. We use the matrix scaling algorithm proposed in [6] in order to reduce the potentially high condition number of the linear system matrix. This scaling technique is applied systematically to all the subsequent linear system matrices, for instance, the left hand side of equation (1.8).
Consider the general linear system of the form $Ax = b$ with A a m×n real matrix. If we denote by $A_S$ the scaled matrix, $b_S$ the scaled second member and $x_S$ the scaled solution, the scaling method consists in solving the scaled problem $A_Sx_S = b_S$ with $A_S = D_R A D_C$, $b_S = D_R b$ and $x_S = D_C^{-1} x$. The two matrices $D_R$ and $D_C$ are the m×m and n×n diagonal matrices defined

by $D_R = \text{diag} \left( \frac{1}{\sqrt{\max_i r_i}} \right)_{i=1,...,m}$ and $D_C = \text{diag} \left( \frac{1}{\sqrt{\max_j c_j}} \right)_{j=1,...,n}$ with $r_i = A_i^T \in \mathbb{R}^{nx1}$, $i = 1, ..., m$, the row-vectors from A and $c_j = A_j^T \in \mathbb{R}^{mx1}$, $j = 1, ..., n$ the column-vectors from A.

Such a matrix scaling is an efficient technique to reduce the condition number of the matrix A when values in A span several orders of magnitude. In our present study, we observed that using one iteration of Algorithm 1 proposed in [6] has shown to be already very efficient also when the rock and fractures have orders of magnitude differences in their permeability/transmissivity properties (see Section 3.4).

## 2.4 Validation: Benchmark Field case

To validate our MATLAB code nef-flow-fpm, we consider the fourth benchmark test case from [13] called [Field]. It consists of 52 fractures embedded in a domain of dimensions $(-500\text{m}, 350\text{m}) \times (100\text{m}, 1500\text{m}) \times (-100\text{m}, 500\text{m})$ as depicted in Figure 3. The data comes from the repository https://git.iws.uni-stuttgart.de/benchmarks/fracture-flow-3d.git.

For both flow and transport simulations, the hydraulic conductivity tensor of the rock matrix is $K_{\alpha}^{\text{m}} = I_3 \text{ m/s}, \forall \alpha$ (the same everywhere). For all fractures, the tangential hydraulic conductivity tensor $K_{\ell}^{\text{f}}$ is constant equal to 10,000 $I_2 \text{ m/s}$ and the scalar aperture $d_{\ell}^{\text{f}}$ is constant equal to 0.01 m. Thus, the fracture transmissivity tensor is also constant and the same for all fractures: $K_{\ell}^{\text{f}} = d_{\ell}^{\text{f}} K_{\ell}^{\text{f}} = 100 I_2 \text{ m}^2/\text{s}, \forall \ell$.

For the flow problem, zero Dirichlet boundary conditions are applied on $\partial\Omega_{\text{in}}$. Neumann conditions are applied on $\partial\Omega_{\text{in}}$ where $q^{\text{f},N}$ is constant equal to $-1 \text{ m/s}$, which induces an inflow of $-1.25 \text{ m}^3/\text{s}$. Zero Neumann conditions are applied everywhere on $\partial\Omega_N = \partial\Omega \setminus (\partial\Omega_{\text{in}} \cup \partial\Omega_{\text{out}})$.

The inlet $\partial\Omega_{\text{in}}$ boundaries and outlet boundaries $\partial\Omega_{\text{out}}$ are shown on Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Benchmark test case field from [13].</figcaption></figure>

For the purely advective transport problem, the final time is 5000s and the time step is chosen equal to 50 s. The matrix porosity $\psi_{\alpha}^{\text{m}}$ is chosen constant for each subdomain $\alpha$ and equal to 0.2. The fracture porosity $\psi_{\ell}^{\text{f}}$ is also chosen constant for each fracture $\ell$ and equal to 0.2. The initial
condition on concentration is zero in $\Omega$, and the boundary condition is a unit concentration at $\partial\Omega_{in}$.

The geometry is downloaded from the repository. The `fracture_network.csv` file contains the data points with the contour of the 52 fractures. The surface mesh of the fracture network and the volume mesh are generated as described in Section 2.1. The fracture network mesh is made of 16,196 triangles and 24,493 edges. The rock matrix mesh contains 263,154 tetrahedra and 533,043 faces. The target number of cells was approximately 260k [13].

The flow problem is solved thanks to our MATLAB code `nef-flow-fpm`. The linear system is solved in 19s with an LU direct solver on [PC]. The hydraulic heads on $K_h$ and the traces of hydraulic head on $\mathcal{F}_h^D \cup \mathcal{F}_h^N$ are shown on Figure 4. We recall that a post-processing phase in `nef-flow-fpm` computes the mean hydraulic heads in the rock matrix, then the velocity fields, both in the rock matrix and the fracture network from their respective mean hydraulic heads and traces of hydraulic head. Note that here the linear system (1.8) is solved with a direct solver so that the qualities of the two velocity fields do not depend on the choice of an iterative solver tolerance, as discussed later on in Section 4.6.

<figure><img src="image_5.png" /><figcaption>Figure 4: Solution computed with `nef-flow-fpm`: (left) hydraulic heads on fracture network; (right) traces of hydraulic head on rock matrix external boundary.</figcaption></figure>

Figure 5 shows the hydraulic heads along two specific lines as in [13]. The results obtained with `nef-flow-fpm` are shown with the black thick line. The results are compared to other methods in [13] using the Python script `compare_my_method.py` from the repository. Our results fall within the published spread, thus validating the hydraulic heads computed with `nef-flow-fpm` on this test case.
<figure><img src="image_6.png" /><figcaption>Figure 5: Hydraulic head profiles across the domain. The shaded region depicts the area between the 10th and the 90th percentile of the published results from [13]. The match number results from evaluating at 100 evenly distributed evaluation points if our result value is between the respective lower and upper value. Top: Line 1. Bottom: Line 2.</figcaption></figure>

In order to validate the velocity fields obtained with nef-flow-fpm we use them as inputs in a purely advective transport problem, solved with our eulerian code nef-transport-fpm. We plot the time series of the average concentrations in the rock matrix and in different fractures. We compare the fracture network results with the ones in [13].

Let the simulation time interval be discretized into $N_t$ points. Let $K_h^\ell$ be the subset of triangles from $K_h$ that discretize the fracture $\ell$. The concentration at time $n \in [1, N_t]$ on triangle $K \in K_h^\ell$ is denoted by $c_{K,\ell}^{\ell,n}$. The time series $\bar{c}_\ell^\ell = (\bar{c}_\ell^{\ell,n})_{n \in [1,N_t]}$ of the average concentration in fracture $\ell$ over time is given by:

$$ \bar{c}_{\ell}^{\ell,n} = \frac{\sum_{K \in \mathcal{K}_{h}^{\ell}} c_{K,\ell}^{\ell,n} |K| d_{\ell}^{\ell} \psi_{\ell}^{\ell}}{\sum_{K \in \mathcal{K}_{h}^{\ell}} |K| d_{\ell}^{\ell} \psi_{\ell}^{\ell}}, \quad \forall n \in [1, N_t]. \quad (2.3) $$

The evolution of the concentration in the fracture network and in the rock matrix at different times is displayed on Figure 6. Figure 7 displays the time series of the average concentration in fractures 16, 46 and 49 (same as in [13]).

Thanks to figures 6 and 7, we show that our results are in good agreement with the ones obtained from [13]. Thus we validate our MATLAB code nef-flow-fpm on this test case.
<figure><img src="image_7.png" /><figcaption>Figure 6: Concentration in the rock matrix and in the fracture network at different times: 100s, 400s, 800s, 1200s, 2000s and 5000s (final time).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Mean concentration over time in selected fractures.</figcaption></figure>

# 3 Generation of test cases with a large number of fractures

## 3.1 Geometrical description

As a reminder, a possible modeling strategy of the network of fractures is the Discrete Fracture Network (DFN) approach. Fractures are modelled as ellipses or disks whose properties (orientation, size, position, transmissivity) are governed by statistical laws derived from field observations [28, 19, 34, 33]. In this work, DFNs are generated using the DFN.lab software, developed by the LabCom Fractory. Examples of such networks are displayed in Figure 8. By considering the
surrounding rock matrix, we obtain a Discrete Fracture Matrix model. Here, the domain under consideration is a cube.

<figure><img src="image_9.png" /><figcaption>Figure 8: Fracture networks of increasing complexity and cube size. Number of fractures: (a) L20: 1.8k; (b) L50: 23k; (c) L60: 39k; (d) L80: 90k.</figcaption></figure>

Flow in such large fracture networks has already been studied in DFNs (impervious rock matrix), for example in [44], with up to 1 million fractures. Nonetheless, such a study has never been carried out within a Discrete Fracture Matrix (DFM) framework, where the 3D rock matrix is considered. The number of fractures in a fractured porous media (FPM) geometry is potentially larger than in an equivalent DFN geometry. This is because, in the DFN case, isolated fractures are not retained, as they do not contribute to flow. This work highlights the ability to mesh such large FPM and to perform physical simulations within them, including flow and transport. Hereafter, the geometries are named according to the domain size (here a cube).

<table><caption>Table 1: DFM “FPM” geometrical properties.</caption><thead><tr><td>Geometry name</td><td>N<sup>fl</sup></td><td>#intersections</td></tr></thead><tbody><tr><td>[L20-FPM]</td><td>1,809</td><td>2,542</td></tr><tr><td>[L50-FPM]</td><td>23,430</td><td>38,528</td></tr><tr><td>[L60-FPM]</td><td>39,456</td><td>66,181</td></tr><tr><td>[L80-FPM]</td><td>90,447</td><td>155,362</td></tr></tbody></table>

#intersections stands for the number of fracture-fracture intersections.

For [L20-FPM] and [L60-FPM], we also consider the equivalent DFN related geometries (impervious rock matrix). As expected, Table 2 shows that, even if the cube size remains the same, the number of fractures is smaller in the DFN geometry compared to the FPM geometry.

<table><caption>Table 2: DFN geometrical properties.</caption><thead><tr><td>Geometry name</td><td>N<sup>fl</sup></td><td>#intersections</td></tr></thead><tbody><tr><td>[L20-DFN]</td><td>1,397</td><td>2,481</td></tr><tr><td>[L60-DFN]</td><td>33,899</td><td>65,485</td></tr></tbody></table>

## 3.2 Meshes

Once we have the geometrical description of a 2D-3D fractured porous medium (FPM) with a large number of fractures, one of the greatest challenges is the mesh generation. Thanks to MODFRAC and GHS3D, in this work it is perfectly possible to mesh such DFNs.
We add to the geometry names found in Section 3.1 either “geo” or “fine/coarse”. “geo” means that the mesh step is automatically set by MODFRAC according to the geometry. “fine” (resp. “coarse”) means that we manually chose a smaller (resp. larger) mesh step. The quality threshold imposed to MODFRAC, equal to $10^{-4}$ and is respected for all the meshes.

Table 3 gives the mesh properties and the time for the mesh generation of FPM geometries with MODFRAC and GHS3D.

<table><caption>Table 3: FPM mesh properties and generation time</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>MODFRAC time hh:mm:ss</th><th>GHS3D time hh:mm:ss</th></tr></thead><tbody><tr><td>[L20geo-FPM]</td><td>1,917,757</td><td>3,860,570</td><td>217,659</td><td>337,699</td><td>00:00:08</td><td>00:00:29</td></tr><tr><td>[L20fine-FPM]</td><td>16,492,814</td><td>33,091,249</td><td>930,339</td><td>1,420,306</td><td>00:00:12</td><td>00:02:22</td></tr><tr><td>[L50geo-FPM]</td><td>5,196,937</td><td>10,447,882</td><td>1,133,823</td><td>1,749,717</td><td>00:00:19</td><td>00:01:49</td></tr><tr><td>[L50fine-FPM]</td><td>16,517,976</td><td>33,140,576</td><td>2,277,976</td><td>3,523,188</td><td>00:00:50</td><td>00:03:49</td></tr><tr><td>[L60geo-FPM]</td><td>7,876,897</td><td>15,825,076</td><td>1,790,168</td><td>2,751,907</td><td>00:00:55</td><td>00:03:02</td></tr><tr><td>[L80geo-FPM]</td><td>16,544,945</td><td>33,206,750</td><td>3,825,145</td><td>5,860,352</td><td>00:02:02</td><td>00:08:25</td></tr></tbody></table>

Table 3 shows that the mesh generation time increases as the fracture network becomes larger, with more fractures and more intersections. Nevertheless, the computational times remain reasonable. The total time, which corresponds to the sum of the time for MODFRAC and the time for GHS3D, is around 10min for the generation of 16.5M of tetrahedral elements of the largest mesh [L80geo-FPM] (90k fractures, 155k intersections).

Tables 4 and 5 give the mesh quality for the above examples with MODFRAC and GHS3D respectively, where $|K|$ (resp. $|T|$) denotes the area (resp. volume) of the triangle $|K|$ (resp. tetrahedral element $T$).

<table><caption>Table 4: Properties of the 2D surface FPM mesh performed with MODFRAC.</caption><thead><tr><th>Mesh name</th><th>mean δ(K)</th><th>min<sub>K</sub> δ(K)</th><th>min<sub>K</sub> |K|</th></tr></thead><tbody><tr><td>[L20geo-FPM]</td><td>0.94</td><td>1.27e-03</td><td>5.74e-07</td></tr><tr><td>[L20fine-FPM]</td><td>0.98</td><td>4.25e-04</td><td>3.03e-09</td></tr><tr><td>[L50geo-FPM]</td><td>0.86</td><td>1.54e-04</td><td>1.39e-09</td></tr><tr><td>[L50fine-FPM]</td><td>0.93</td><td>1.40e-04</td><td>1.39e-09</td></tr><tr><td>[L60geo-FPM]</td><td>0.85</td><td>1.05e-04</td><td>1.39e-09</td></tr><tr><td>[L80geo-FPM]</td><td>0.84</td><td>1.11e-04</td><td>3.58e-09</td></tr></tbody></table>

<table><caption>Table 5: Properties of the 3D volume FPM mesh performed with GHS3D.</caption><thead><tr><th>Mesh name</th><th>mean δ(T)</th><th>min<sub>T</sub> δ(T)</th><th>min<sub>T</sub> |T|</th><th>median δ(T)</th><th>N(T<sub>h</sub>) s.t. δ(T) &lt; 0.1</th><th>% N(T<sub>h</sub>) s.t. δ(T) &lt; 0.1</th></tr></thead><tbody><tr><td>[L20geo-FPM]</td><td>0.76</td><td>4.43e-05</td><td>4.54e-10</td><td>0.78</td><td>4.4k</td><td>0.23%</td></tr><tr><td>[L20fine-FPM]</td><td>0.78</td><td>2.47e-06</td><td>1.22e-11</td><td>0.79</td><td>11.0k</td><td>0.07%</td></tr><tr><td>[L50geo-FPM]</td><td>0.64</td><td>4.41e-07</td><td>1.19e-12</td><td>0.67</td><td>50.5k</td><td>0.97%</td></tr><tr><td>[L50fine-FPM]</td><td>0.74</td><td>7.87e-07</td><td>2.12e-13</td><td>0.76</td><td>50.0k</td><td>0.30%</td></tr><tr><td>[L60geo-FPM]</td><td>0.62</td><td>1.69e-06</td><td>1.53e-12</td><td>0.65</td><td>88.9k</td><td>1.13%</td></tr><tr><td>[L80geo-FPM]</td><td>0.61</td><td>6.90e-08</td><td>1.36e-13</td><td>0.63</td><td>216k</td><td>1.30%</td></tr></tbody></table>
Despite the fact that surface meshes respect the minimum quality of $10^{-4}$, as seen in Table 4, Table 5 shows that the volume meshes have a lower minimum quality. For instance, for [L80geo-FPM], the quality is below $10^{-7}$. Moreover, as shown in Table 5, as the number of fractures increases, there is a slight decrease in the median quality value and a slight increase in the number of tetrahedral elements with a quality below 0.1.

Nevertheless, the number of low quality elements in the 3D mesh ($\delta(T) < 0.1$) remains relatively small, equal or less than 1.30% of the total number of tetrahedra. As a matter of fact, in Section 4 we show that these FPM meshes perform correctly for flow and transport simulations.

For [L20-FPM] and [L60-FPM] geometries we also consider the related “geometry” meshes:

* [L20geo-DFN] and [L60geo-DFN], 2D only, impervious matrix, as presented before;

* [L20coarse-PM-GQ] and [L60geo-PM-GQ], 3D only (“PM” stands for “Porous Media”, the rock matrix only), with a good quality (GQ) mesh without constraints.

* [L20geo-PM-LQ] and [L60geo-PM-LQ], 3D only again, but this time with a low quality (LQ) mesh. The low quality is achieved by keeping the fractures as geometric constraints, resulting in the same mesh as in the FPM case.

Table 6 gives the PM and DFN mesh properties and the time for the mesh generation.

<table><caption>Table 6: “PM” and “DFN” mesh properties and generation time</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>MODFRAC time hh:mm:ss</th><th>GHS3D time hh:mm:ss</th></tr></thead><tbody><tr><td>[L20geo-DFN]</td><td>0</td><td>0</td><td>178,377</td><td>276,114</td><td>00:00:38</td><td>-</td></tr><tr><td>[L20coarse-PM-GQ]</td><td>618,428</td><td>1,258,229</td><td>0</td><td>0</td><td>00:00:00.35</td><td>00:00:06</td></tr><tr><td>[L20geo-PM-LQ]</td><td>1,917,757</td><td>3,860,570</td><td>0</td><td>0</td><td>00:00:08</td><td>00:00:29</td></tr><tr><td>[L60geo-DFN]</td><td>0</td><td>0</td><td>1,726,854</td><td>2,637,594</td><td>00:00:24</td><td>-</td></tr><tr><td>[L60geo-PM-GQ]</td><td>672,721</td><td>1,368,422</td><td>0</td><td>0</td><td>00:00:00.3</td><td>00:00:06</td></tr><tr><td>[L60geo-PM-LQ]</td><td>7,876,897</td><td>15,825,076</td><td>0</td><td>0</td><td>00:00:55</td><td>00:03:02</td></tr></tbody></table>

<table><caption>Table 7: Properties of the 2D surface “DFN” mesh performed with MODFRAC.</caption><thead><tr><th>Mesh name</th><th>mean δ(K)</th><th>min<sub>K</sub> δ(K)</th><th>min<sub>K</sub> |K|</th></tr></thead><tbody><tr><td>[L20geo-DFN]</td><td>0.95</td><td>5.65e-04</td><td>5.74e-07</td></tr><tr><td>[L60geo-DFN]</td><td>0.85</td><td>1.05e-04</td><td>2.01e-08</td></tr></tbody></table>

<table><caption>Table 8: Properties of the 3D volume “PM” mesh performed with GHS3D.</caption><thead><tr><th>Mesh name</th><th>mean δ(T)</th><th>min<sub>T</sub> δ(T)</th><th>min<sub>T</sub> |T|</th><th>median δ(T)</th><th>N(T<sub>h</sub>) s.t. δ(T) < 0.1</th><th>% N(T<sub>h</sub>) s.t. δ(T) < 0.1</th></tr></thead><tbody><tr><td>[L20coarse-PM-GQ]</td><td>0.81</td><td>3.19e-01</td><td>1.38e-03</td><td>0.82</td><td>0</td><td>0%</td></tr><tr><td>[L20geo-PM-LQ]</td><td>0.76</td><td>4.43e-05</td><td>4.54e-10</td><td>0.78</td><td>4.4k</td><td>0.23%</td></tr><tr><td>[L60geo-PM-GQ]</td><td>0.81</td><td>3.21e-01</td><td>3.93e-02</td><td>0.82</td><td>0</td><td>0%</td></tr><tr><td>[L60geo-PM-LQ]</td><td>0.62</td><td>1.69e-06</td><td>1.53e-12</td><td>0.65</td><td>88.9k</td><td>1.13%</td></tr></tbody></table>

For “PM” meshes there is a MODFRAC time because we generate the 2D surface mesh of the cube faces with the mentioned software. Furthermore, specifically for the “PM-LQ” meshes, we also
have the geometrical constraints inside the cube from the fracture network. Tables 7 and 8 give the mesh quality of the three mentioned meshes. We can see that there is indeed a large difference in quality between the “PM-GQ” and “PM-LQ” meshes.

## 3.3 Flow hydrogeological setting

### 3.3.1 Boundary conditions and source terms

In this work, we use permeameter boundary conditions (BCs) [35], both in the rock matrix and in the fracture network. They work as follows:

We impose a gradient of the hydraulic heads on two opposite sides of the cube. Thus, in the fracture network, $\Sigma^D$ is the union of two sets of edges $\Sigma^{D,in}$ and $\Sigma^{D,out}$ living on two opposite faces of the cube: $\Sigma^D = \Sigma^{D,in} \cup \Sigma^{D,out}$. Also, in the rock matrix, $\Gamma^D$ is the union of two opposite cube faces $\Gamma^{D,in}$ and $\Gamma^{D,out}$.

At the same time, no flux boundary conditions are imposed in the other faces of the cube. Thus, $\Sigma^N$ and $\Gamma^N$ are defined as: $\Sigma^N = (\partial\Omega^f)\cap\Sigma^D$ and $\Gamma^N = (\partial\Omega^m)\cap\Gamma^D$.

Thanks to this configuration, we obtain a rock matrix flow and a fracture network flow in the $x$-direction. See an illustration on figure 9.

<figure><img src="image_10.png" /><figcaption>Figure 9: Permeameter boundary conditions with a flow in the $x$-direction. We give an example of a network with 3 fractures.</figcaption></figure>

Thus, in this work we apply the same Dirichlet and Neumann BCs conditions to the rock matrix and to the fracture network. We choose the following constant BC functions: $p^{m,D,in} = 10$ m, $p^{m,D,out} = 0$, $p^{f,D,in} = 10$ m, $p^{f,D,out} = 0$, $q^{m,N} = 0$ and $q^{f,N} = 0$.

With respect to source terms, even if nef-flow-fpm is able to handle any source term function $f^m$ and $f^f$, in this work we choose them to be equal to 0 in the rock matrix and in the fracture network.

### 3.3.2 Rock matrix hydraulic conductivity and fracture network transmissivity

One of the goals of this work is to assess the influence of permeability/transmissivity contrast on solver efficiency. Thus, we choose two settings for the hydraulic conductivity in the rock matrix and the transmissivity in the fractures in the “FPM” geometries :

* “1v100”, where the hydraulic conductivity in the rock matrix is chosen equal to $K_\alpha^m = I_3$ m/s, $\forall\alpha$, and where the fracture tangential transmissivity is chosen equal to $K_\ell^f = 100I_2$ m$^2$/s, $\forall\ell$ (the same for all fractures);

* “heter”, where the hydraulic conductivity in the rock matrix is chosen equal to $K_\alpha^m = 10^{-8}I_3$ m/s, $\forall\alpha$, and where the fracture tangential transmissivity $K_\ell^f$ takes one value
per fracture $ℓ$ in the range [$10^{-6}$, $20$] m$^2$/s. The transmissivity data is read from a file given by the LabCom Fractory.

Exclusively for the FPM case, we denote the maximum ratio of heterogeneity by $K_{ratio}^{m,f}$. For the “1v100” setting, this ratio is equal to 100. For the “heter” setting, even if the exact value depends on the chosen geometry, it is always close to $10^7$.

In this work, we choose the fractures to be more hydraulically conductive than the surrounding rock matrix. The goal is to make the fractures preferential channels for flow, avoiding the case where fractures are considered as barriers for the flow in the rock matrix. This is consistent with the assumption that the hydraulic head is continuous between the rock matrix and the fractures (Equation (1.5)).

The histograms of the transmissivities with the “heter” setting applied to [L20-FPM] and [L60-FPM] geometries are displayed on figure 10.

<figure><img src="image_11.png" /><figcaption>Figure 10: Histogram of transmissivities of the “heter” setting applied to different geometries.</figcaption></figure>

For the “PM” case, the hydraulic conductivity in the rock matrix is set to $K = I_3$ m/s, and it is indicated in the simulation name with “1”.

For the “DFN” case, the fracture tangential transmissivity is chosen equal to $K_{ℓ}^{f} = 100 I_2$ m$^2$/s, $\forall ℓ$ (the same for all fractures), and it is indicated in the simulation name with “100”.

## 3.4 Large linear systems and condition number estimation

After choosing a geometry, a mesh and a permeability/transmissivity setting and boundary conditions (here permeameter), nef-flow-fpm is able to assemble the following linear systems, each one corresponding to a “test case”. They can be found in table 9.

We propose to estimate the condition number (non-scaled and scaled) of the matrix of the linear system (1.8) for all the presented test cases, whenever the computational ressources allow it. The scaling technique is the one described in subsection 2.3. The estimated condition numbers (non-scaled and scaled) are obtained with MATLAB’s routine condest. Table 9 gathers the main characteristics of the various matrices.

For the test cases related to [L20-FPM] geometry, as expected from the values of $K_{ratio}^{m,f}$ ($10^{-2}$ versus $10^7$), the condition number in the case [L20geo-FPM-heter] is 5 orders of magnitude larger than the one in the case [L20geo-FPM-1v100] [54]. The condition number also depends on the mesh quality as suggested in subsection 3.2: the lower the mesh quality, the higher the condition number.
Table 9 also shows that the scaling method, described in subsection 2.3, allows to significantly reduce the condition number, especially for the “heter” test cases: up to 5 orders of magnitude for [L20geo-FPM-heter] and 4 orders of magnitude for [L60geo-FPM-1v100].

These estimations of the condition number on the test cases related to [L20-FPM] and [L60-FPM] geometries underline how ill-conditioned the linear systems can be. The larger and more heterogeneous the test cases are, the more ill-conditioned the associated linear systems are. Moreover, the scaling method we are using helps in reducing the condition number.

<table><caption>Table 9: Properties of the linear systems and their estimated condition numbers.</caption><thead><tr><th>Test case</th><th># dofs</th><th>nnz</th><th>condest<br>w/o scaling</th><th>condest<br>with scaling</th></tr></thead><tbody><tr><td>[L20fine-FPM-1v100]</td><td>34,437,027</td><td>243,067,445</td><td>**</td><td>**</td></tr><tr><td>[L20geo-FPM-heter]</td><td>4,179,938</td><td>29,694,870</td><td>6.86e+16</td><td>6.25e+11</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>4,179,938</td><td>29,694,870</td><td>1.71e+11</td><td>2.13e+07</td></tr><tr><td>[L20phy-PM-GQ-1]</td><td>1,243,973</td><td>26,872,005</td><td>1.55e+05</td><td>7.70e+04</td></tr><tr><td>[L20geo-PM-LQ-1]</td><td>3,860,470</td><td>8,579,364</td><td>1.01e+09</td><td>6.22e+07</td></tr><tr><td>[L20geo-DFN-100]</td><td>274,708</td><td>1,339,346</td><td>4.16e+08</td><td>5.47e+06</td></tr><tr><td>[L50geo-FPM-1v100]</td><td>12,158,838</td><td>87,895,320</td><td>4.57e+13</td><td>1.61e+08</td></tr><tr><td>[L50fine-FPM-1v100]</td><td>36,586,026</td><td>261,670,914</td><td>**</td><td>**</td></tr><tr><td>[L60geo-FPM-heter]</td><td>18,525,914</td><td>134,224,135</td><td>**</td><td>**</td></tr><tr><td>[L60geo-FPM-1v100]</td><td>18,525,914</td><td>134,224,135</td><td>5.35e+12</td><td>1.73e+08</td></tr><tr><td>[L60geo-PM-GQ-1]</td><td>1,353,104</td><td>9,333,699</td><td>1.73e+05</td><td>8.54e+04</td></tr><tr><td>[L60geo-PM-LQ-1]</td><td>15,825,049</td><td>110,347,069</td><td>5.27e+11</td><td>2.37e+09</td></tr><tr><td>[L60geo-DFN-100]</td><td>2,628,345</td><td>12,952,420</td><td>4.56e+09</td><td>5.37e+07</td></tr><tr><td>[L80geo-FPM-1v100]</td><td>38,986,205</td><td>282,941,795</td><td>**</td><td>**</td></tr></tbody></table>

The symbol ** in Table 9 indicates that the computation of the estimation of the condition number requires a RAM greater than 3TB (the maximum available in [MEM]) and therefore cannot be performed with the available resources.

Thanks to specialized software like MODFRAC and GHS3D, as well as nef-flow-fpm, this work demonstrates the ability to assemble 2D-3D single-phase flow test cases in fractured porous media with a large number of fractures. For instance, the test case [L80geo-FPM-1v100], created from the [L80-FPM] geometry with 90k fractures, is the largest linear system with ~39M dofs. Such linear systems are very ill-conditioned, in particular under the “heter” permeability/transmissivity configuration. For example, the estimation of the condition number of [L20geo-FPM-heter] test case reaches $7 \cdot 10^{16}$ without scaling and $6 \cdot 10^{11}$ with scaling.

# 4 Performance of linear solvers with respect to large-scale simulations

In this section, we present a variety of numerical experiments on flow in fractured porous media. Given the large number of test cases and the unique properties of each case, we systematically summarize their key characteristics in each subsection. For full details on these properties, readers should refer back to Section 3.

## 4.1 Hardware

Depending on the problem size and solvers, the simulations presented from now are run on different computational resources:
* on the local Inria cluster, CLEPS, on two different partitions:

- the first one with large RAM memory: 4x Intel Xenon, E7-4860 v2; 12 cores, 2.6-3.2 GHz; RAM: 2.95 TB, 1,333 MHz, hereafter referred to as [MEM]

- the second with 2x Cascade Lake Intel Xeon 5218, 16 cores, 2.4GHz, 192 GB, 2,667 MHz, hereafter referred to as [HOMO].

* on supercomputer Joliot-Curie located at CEA's Very Large Computing Center (TGCC), AMD Irene ROME, hereafter referred to as [ROME], 2292 dual-processor AMD Epyc Rome computer nodes at 2.6 GHz with 64 cores per processor, for a total of 293,376 computing cores and a power of 11.75 Pflop/s, 256 GB DDR4 memory/node.

We do not use hyper threading in our simulations.

## 4.2 Direct solvers

As the linear system (1.8) is square, sparse, symmetric, positive and definite, as a first test, we propose to assess direct solvers that can benefit from the structure of the linear system. The most natural direct solver in our case would be the Cholesky factorization (called CHOL in the following). We also use the LU decomposition (called LU in the following) despite the linear system symmetry, since some LU implementations, like the one in MATLAB, have shown more robustness than the CHOL implementation with respect to floating-point arithmetic.

We recall that nef-flow-fpm allows to save the linear system in PETSc [10, 9] binary format. Thus, in order to test CHOL and LU, the linear system resolution is done externally in parallel with MUMPS [5, 4] via a PETSc interface.

In order to assess these two direct solvers, we choose the following test cases: [L20geo-FPM-1v100] and [L60geo-FPM-1v100]. We can find the details of their geometry, hydraulic conductivity in the rock matrix, transmissivity in the fracture network, boundary conditions, # dofs and condition number estimation in Section 3. We recall their main properties in table 10.

<table><caption>Table 10: Main properties of [L20geo-FPM-1v100] and [L60geo-FPM-1v100] test cases.</caption><thead><tr><th>Test case</th><th>N<sup>fl</sup></th><th>N(T<sub>h</sub>)</th><th>K<sup>rm,fl</sup><sub>ratio</sub></th><th># dofs</th><th>condest<br>w/o scaling</th></tr></thead><tbody><tr><td>[L20geo-FPM-1v100]</td><td>1.8k</td><td>1.9M</td><td>10<sup>2</sup></td><td>4.2M</td><td>2 · 10<sup>11</sup></td></tr><tr><td>[L60geo-FPM-1v100]</td><td>40k</td><td>7.9M</td><td>10<sup>2</sup></td><td>19M</td><td>5 · 10<sup>12</sup></td></tr></tbody></table>

The time and computational resources (RAM) of the simulations run in parallel with 48 MPI processes can be found in Table 11.

<table><caption>Table 11: Time and computational resources (RAM) of [L20geo-FPM-1v100] and [L60geo-FPM-1v100] test cases.</caption><thead><tr><th>Test case</th><th>Solver</th><th>Hardware</th><th>Solver time<br>hh:mm:ss</th><th>RAM<br>(GB)</th></tr></thead><tbody><tr><td>[L20geo-FPM-1v100]</td><td>CHOL</td><td>[ROME]</td><td>00:00:57</td><td>27.3</td></tr><tr><td>[L60geo-FPM-1v100]</td><td>CHOL</td><td>[ROME]</td><td>00:06:26</td><td>107</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>LU</td><td>[ROME]</td><td>00:01:21</td><td>41.3</td></tr><tr><td>[L60geo-FPM-1v100]</td><td>LU</td><td>[ROME]</td><td>00:08:50</td><td>159</td></tr></tbody></table>

We read the RAM usage from the MUMPS log file at the INFOG(22) line, which represents the memory effectively used during factorization, summed over all 48 processors. The indicated
solver time is the total “main stage” execution time read from PETSc log file, accounting almost exclusively for MUMPS clock time.

As expected, CHOL has a lower solver time and uses significantly fewer computational resources than LU, as CHOL takes advantage of the symmetric nature of the linear system.

For [L20geo-FPM-1v100], the smallest test case with 1.8k fractures, both CHOL and LU have reasonable solver times. However, [L60geo-FPM-1v100], which is a larger test case with 40k fractures, requires significantly more computational resources (RAM), and the solver time is no longer as favorable, especially with LU.

Direct solvers would be impractical for even larger test cases, such as the one corresponding to the [L80-FPM] geometry with 90k fractures, or [L20geo-FPM-heter], which features a very strong $K_{ratio}^{m,f}$ contrast, leading to a high condition number of the linear system matrix. We now turn to iterative solvers, which require fewer computational resources and can potentially offer more favorable CPU times.

## 4.3 Conjugate Gradient algebraic preconditioners

We consider the Preconditioned Conjugate Gradient (PCG) as the iterative solver to be tested. It is implemented in MATLAB and in AMGCL's C++ library [36, 37, 38].

From each implementation we choose one of the following preconditioners, namely: MATLAB's Incomplete Cholesky factorization (hereafter named PCG-iCHOL) and AMGCL's Algebraic Multi-Grid (hereafter named PCG-AMGCL). The latter is an external library encoded in a MEX-file, thanks to MRST, integrated directly to nef-flow-fpm.

In order to assess these preconditioners, we choose the test case [L20geo-FPM-1v100]. Its properties can be found in the previous subsection 4.2, Table 10.

The PCG tolerance is set to $10^{-14}$. At first glance, this value might seem very demanding, but this choice is thoroughly discussed later in Section 4.6. Besides this parameter, we set the maximum number of iterations to 10,000. Furthermore, each preconditioner has its own set of parameters in order to run properly:

* **PCG-iCHOL:** in order to avoid PCG breakdowns, we use the `diacomp` option that constructs a shifted incomplete Cholesky factor of matrix $A+m \times \text{diag}(\text{diag}(A))$, where A is the linear system matrix from equation (1.8) and where $m = \max(\text{sum}(\text{abs}(A), 2) / \text{diag}(A)) - 2$ (in MATLAB notation); the simulations are sequential;

* **PCG-AMGCL:** `nef-flow-fpm` parses the parameters to AMGCL C++'s backend presented in appendix A.3; in particular, we run the simulations sequentially; in this work, these parameters are the same for all simulations.

Table 12 presents the performances of the mentioned PCG preconditioners: number of iterations and time. Figure 11 displays the relative residual versus the number of iterations, truncated at 4,000 iterations, for PCG-iCHOL and PCG-AMGCL.

<table><caption>Table 12: Time and computational resources (RAM) of PCG with different preconditioners for test case [L20geo-FPM-1v100].</caption><thead><tr><th>Precond.</th><th>Hardware</th><th># it.</th><th>Solver time<br>hh:mm:ss</th></tr></thead><tbody><tr><td>iCHOL</td><td>[HOMO]</td><td>9,972 &lt;!&gt;</td><td>02:55:36 &lt;!&gt;</td></tr><tr><td>AMGCL</td><td>[HOMO]</td><td>2,869</td><td>00:28:01</td></tr></tbody></table>

The `<!&gt;` symbol indicates that PCG has not converged in 9,972 iterations to the requested tolerance $10^{-14}$ because two iterations were the same. The simulation stopped at iteration 9,972 with the following residual: $7.26 \cdot 10^{-13}$.
For PCG-iCHOL and PCG-AMGCL, the solver time is read from nef-flow-fpm log file at "solver time" line.

We do not show the RAM for these three preconditioners since there is no reliable way to compare them. As a reference, we give the RAM from the AMG preconditioner's memory footprint computed by AMGCL's C++ backend, which is 1.62 GB.

<figure><img src="image_12.png" /><figcaption>Figure 11: Convergence of PCG for test case [L20geo-FPM-1v100] with different preconditioners, namely PCG-iCHOL and PCG-AMGCL: relative residual (log-scale) versus the number of iterations.</figcaption></figure>

When it comes to RAM, iterative methods are much better than direct methods.

However, MATLAB's PCG-iCHOL has a very slow convergence rate. It does not even converge to the desired tolerance of $10^{-14}$, although the final residual is small equal to $7.28 \cdot 10^{-13}$. It was not surprising that iCHOL did not work well for such a difficult linear system. On the other hand, PCG-AMGCL converges down to $10^{-14}$

In the following, we only consider AMGCL as the preconditioner for the iterative method PCG.

## 4.4 Geometry complexity impact on solver convergence

We now analyze the effect of the fracture network geometric complexity on the convergence of the iterative solver PCG preconditioned with AMGCL.

We consider three test cases: [L20fine-FPM-1v100], [L50fine-FPM-1v100] and [L80geo-FPM-1v100], whose descriptions can be found in Section 3. As a reminder, we recall their main characteristics in Table 13.

These test cases are characterized by an increasing number of fractures and cube size. Each mesh is built so as to contain roughly 17M tetrahedral elements. As a result, the linear system size for the three cases is roughly the same, between 34M and 39M. This setup allows us to compare the effect of the fracture network geometry as the linear systems have roughly the same sizes.

<table><caption>Table 13: Main properties of [L20fine-FPM-1v100], [L50fine-FPM-1v100] and [L80geo-FPM-1v100] test cases.</caption><thead><tr><th>Test case</th><th>N<sup>f</sup></th><th>N(T<sub>h</sub>)</th><th>K<sup>m,f</sup><sub>ratio</sub></th><th># dofs</th></tr></thead><tbody><tr><td>[L20fine-FPM-1v100]</td><td>1.8k</td><td>17M</td><td>10<sup>2</sup></td><td>34M</td></tr><tr><td>[L50fine-FPM-1v100]</td><td>23k</td><td>17M</td><td>10<sup>2</sup></td><td>37M</td></tr><tr><td>[L80geo-FPM-1v100]</td><td>90k</td><td>17M</td><td>10<sup>2</sup></td><td>39M</td></tr></tbody></table>
The simulations are run with PCG-AMGCL with the same parameters as in Section 4.3, except that the maximum number of iterations is set to 50,000.

Figure 12 displays the relative residual versus the number of iterations for the three test cases. Table 14 shows the number of iterations at given residual values. Table 15 shows the computational times and the consumption of computational resources (RAM).

<figure><img src="image_13.png" /><figcaption>Figure 12: Convergence of PCG-AMGCL for test cases [L20fine-FPM-1v100], [L50fine-FPM-1v100] and [L80geo-FPM-1v100]: relative residual (log-scale) versus the number of iterations.</figcaption></figure>

<table><caption>Table 14: Number of iterations at given residuals with PCG-AMGCL for test cases [L20fine-FPM-1v100], [L50fine-FPM-1v100] and [L80geo-FPM-1v100].</caption><thead><tr><th>Test case</th><th>10<sup>-2</sup></th><th>10<sup>-4</sup></th><th>10<sup>-6</sup></th><th>10<sup>-8</sup></th><th>10<sup>-10</sup></th><th>10<sup>-12</sup></th><th>10<sup>-14</sup></th></tr></thead><tbody><tr><td>[L20fine-FPM-1v100]</td><td>11</td><td>106</td><td>301</td><td>909</td><td>1,324</td><td>2085</td><td>2,252</td></tr><tr><td>[L50fine-FPM-1v100]</td><td>11</td><td>169</td><td>594</td><td>2,030</td><td>4,039</td><td>4,502</td><td>5,325</td></tr><tr><td>[L80geo-FPM-1v100]</td><td>12</td><td>250</td><td>1,125</td><td>3,768</td><td>20,324</td><td>34,505</td><td>39,464</td></tr></tbody></table>

<table><caption>Table 15: Time and computational resources (RAM) of PCG-AMGCL, with tolerance $10^{-14}$, for test cases [L20fine-FPM-1v100], [L50fine-FPM-1v100] and [L80geo-FPM-1v100].</caption><thead><tr><th>Test case</th><th>Hardware</th><th>#it.</th><th>Solver time<br>dd-hh:mm:ss</th><th>RAM<br>(GB)</th></tr></thead><tbody><tr><td>[L20fine-FPM-1v100]</td><td>[HOMO]</td><td>2,252</td><td>00-03:09:38</td><td>13.38</td></tr><tr><td>[L50fine-FPM-1v100]</td><td>[HOMO]</td><td>5,325</td><td>00-08:09:28</td><td>14.32</td></tr><tr><td>[L80geo-FPM-1v100]</td><td>[HOMO]</td><td>39,464</td><td>02-18:36:54</td><td>15.73</td></tr></tbody></table>

The indicated solver time is read from nef-flow-fpm log file at "solver time" line. We read the RAM from the AMG preconditioner's memory footprint computed by AMGCL's C++ backend.

As seen on Table 15, the three chosen test cases converge down to $10^{-14}$.

However, in Figure 12, we see that the convergence of PCG-AMGCL is slowed down as the complexity of the fracture network increases, for a given linear system size. This is visible especially with the largest test case [L80geo-FPM-1v100], with 90k fractures and 39M dofs, where both the number of iterations and the solver time are very important. The RAM remains rather constant for all test cases and much smaller with respect to direct solvers, as shown in Table 15.

Furthermore, depending on the desired precision on the solution, the number of iterations can be prohibitively high. Typically, for [L80geo-FPM-1v100], the largest test case, according to
Table 14, PCG-AMGCL needs more than 20k iterations to reach a tolerance of $10^{-10}$ and almost 40k iterations to converge down to $10^{-14}$.

We conclude that linear iterative solvers, with the help of algebraic preconditioners, can solve the linear system (1.8) corresponding to single-phase flow in three-dimensional fractured porous media with a large number of fractures, up to 90k. Even though it is possible for very strict tolerances as low as $10^{-14}$, the iterative solver struggles to converge for the largest test case. Note such large-scale simulations present at least two challenges for black-box multigrid preconditioners. First, the presence of a small number of low-quality cells in the mesh may degrade convergence, as recently investigated by [31]. Second, black-box implementations of multigrid solvers may be affected by the chosen hybridized discretization. Future work could explore more advanced implementations with proper intergrid transfer operators designed for the hybridized version of mixed methods, as proposed in [26, 52, 90].

## 4.5 Mesh quality and ratio of heterogeneity impact on convergence

We investigate the impact of the mesh quality and the effect of the ratio of heterogeneity between the rock matrix hydraulic conductivity and the fracture network transmissivity, $K_{\text{ratio}}^{\text{m,f}}$, defined in subsection 3.3.2, on the number of iterations of the iterative solver PCG preconditioned with AMGCL.

We recall some notation. The test cases denominated with “PM” correspond to exclusively 3D rock matrix porous flow simulations, without any fractures. “FPM” test cases correspond to porous-fractured simulations with 2D and 3D (fracture network and rock matrix). “DFN” test cases correspond to exclusively 2D fracture network flow simulations. The test cases denominated with “PM-GQ” are “PM” simulations with a good quality mesh. The test cases denominated with “PM-LQ” are “PM” simulations with a low quality mesh, which is simply the volume mesh from the corresponding “FPM” test case. Only for “FPM” simulations, $K_{\text{ratio}}^{\text{m,f}}$ is the contrast between the hydraulic conductivity in the rock matrix and the transmissivity in the fracture network.

We consider two cube size L20 and L60 with their respective test cases : [L20geo-DFN-100], [L20coarse-PM-GQ-1], [L20geo-PM-LQ-1], [L20geo-FPM-1v100], [L20geo-FPM-heter] and [L60geo-DFN-100], [L60geo-PM-GQ-1], [L60geo-PM-LQ-1], [L60geo-FPM-1v100], [L60geo-FPM-heter]. Their description is given in Section 3. Their main characteristics are recalled in Table 16.

<table><caption>Table 16: Main properties of [L20-FPM] and [L60-FPM] related test cases.</caption><thead><tr><th>Test case</th><th>N<sup>f</sup></th><th>N(T<sub>h</sub>)</th><th>min<sub>T</sub><br>δ(T)</th><th>N(K<sub>h</sub>)</th><th>min<sub>K</sub><br>δ(K)</th><th>K<sub>ratio</sub><sup>m,f</sup></th><th>condest<br>w/o scaling</th><th># dofs</th></tr></thead><tbody><tr><td>[L20geo-DFN-100]</td><td>1.4k</td><td>0</td><td>-</td><td>0.18M</td><td>6 · 10<sup>-4</sup></td><td>-</td><td>4 · 10<sup>8</sup></td><td>0.27M</td></tr><tr><td>[L20coarse-PM-GQ-1]</td><td>0</td><td>0.62M</td><td>3 · 10<sup>-1</sup></td><td>0</td><td>-</td><td>-</td><td>2 · 10<sup>5</sup></td><td>1.2M</td></tr><tr><td>[L20geo-PM-LQ-1]</td><td>0</td><td>1.9M</td><td>4 · 10<sup>-5</sup></td><td>0</td><td>-</td><td>-</td><td>1 · 10<sup>9</sup></td><td>3.9M</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>1.8k</td><td>1.9M</td><td>4 · 10<sup>-5</sup></td><td>0.22M</td><td>1 · 10<sup>-3</sup></td><td>10<sup>2</sup></td><td>2 · 10<sup>11</sup></td><td>4.2M</td></tr><tr><td>[L20geo-FPM-heter]</td><td>1.8k</td><td>1.9M</td><td>4 · 10<sup>-5</sup></td><td>0.22M</td><td>1 · 10<sup>-3</sup></td><td>10<sup>7</sup></td><td>7 · 10<sup>16</sup></td><td>4.2M</td></tr><tr><td>[L60geo-DFN-100]</td><td>34k</td><td>0</td><td>-</td><td>1.7M</td><td>1 · 10<sup>-4</sup></td><td>-</td><td>5 · 10<sup>9</sup></td><td>2.6M</td></tr><tr><td>[L60geo-PM-GQ-1]</td><td>0</td><td>0.67M</td><td>3 · 10<sup>-1</sup></td><td>0</td><td>-</td><td>-</td><td>2 · 10<sup>5</sup></td><td>1.4M</td></tr><tr><td>[L60geo-PM-LQ-1]</td><td>0</td><td>7.9M</td><td>2 · 10<sup>-6</sup></td><td>0</td><td>-</td><td>-</td><td>5 · 10<sup>11</sup></td><td>16M</td></tr><tr><td>[L60geo-FPM-1v100]</td><td>40k</td><td>7.9M</td><td>2 · 10<sup>-6</sup></td><td>1.8M</td><td>1 · 10<sup>-4</sup></td><td>10<sup>2</sup></td><td>5 · 10<sup>12</sup></td><td>19M</td></tr><tr><td>[L60geo-FPM-heter]</td><td>40k</td><td>7.9M</td><td>2 · 10<sup>-6</sup></td><td>1.8M</td><td>1 · 10<sup>-4</sup></td><td>10<sup>7</sup></td><td>**</td><td>19M</td></tr></tbody></table>

The symbol ** indicates that the computation of the estimation of the condition number requires a RAM greater than 3TB (the maximum available in [MEM]) and therefore cannot be performed with the available resources.
The simulations are run with PCG-AMGCL with the same parameters as in Section 4.3. In particular, PCG tolerance is set to $10^{-14}$ and the maximum number of iterations is equal to 250,000.

Table 17 shows the computational time, the number of iterations and the consumption of computational resources for the test cases presented in Table 16. Figure 13 displays the relative residual versus the number of iterations. For [L60-FPM], the number of iterations is truncated at 150k.

<table><caption>Table 17: Time and computational resources (RAM) of PCG-AMGCL for test cases related to [L20-FPM] and [L60-FPM] geometries. PCG tolerance is set to $10^{-14}$.</caption><thead><tr><th>Test case</th><th># it.</th><th>Hardware</th><th>Solver time<br>dd-hh:mm:ss</th><th>RAM<br>(GB)</th></tr></thead><tbody><tr><td>[L20geo-DFN-100]</td><td>692</td><td>[HOMO]</td><td>00-00:00:17</td><td>0.08</td></tr><tr><td>[L20coarse-PM-GQ-1]</td><td>66</td><td>[HOMO]</td><td>00-00:00:18</td><td>0.497</td></tr><tr><td>[L20geo-PM-LQ-1]</td><td>2,929</td><td>[HOMO]</td><td>00-00:26:18</td><td>1.49</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>2,869</td><td>[HOMO]</td><td>00-00:28:01</td><td>1.62</td></tr><tr><td>[L20geo-FPM-heter]</td><td>122,798</td><td>[HOMO]</td><td>00-20:01:26</td><td>1.65</td></tr><tr><td>[L60geo-DFN-100]</td><td>3,082</td><td>[HOMO]</td><td>00-00:16:35</td><td>0.76</td></tr><tr><td>[L60geo-PM-GQ-1]</td><td>69</td><td>[HOMO]</td><td>00-00:00:18</td><td>0.51</td></tr><tr><td>[L60geo-PM-LQ-1]</td><td>5,121</td><td>[HOMO]</td><td>00-03:16:46</td><td>6.33</td></tr><tr><td>[L60geo-FPM-1v100]</td><td>5,249</td><td>[HOMO]</td><td>00-04:02:27</td><td>7.45</td></tr><tr><td>[L60geo-FPM-heter]</td><td>215,072 ××</td><td>[HOMO]</td><td>07-00:00:00 ××</td><td>××</td></tr></tbody></table>

The symbol ×× indicates the simulation stopped at the requested time limit for a job in [HOMO], which is 7 days. At that time, the lowest residual was equal to $4.93 \cdot 10^{-12}$ and was reached after 212,237 iterations. Since the simulation stopped, RAM cannot be indicated.

The indicated solver time is read from nef-flow-fpm log file as "solver time" line. We read the RAM from the AMG preconditioner's memory footprint computed by AMGCL's C++ backend.

<figure><img src="image_14.png" /><figcaption>Figure 13: Convergence of PCG-AMGCL - relative residual (log-scale) versus the number of iterations (log-scale).</figcaption></figure>

For the “PM-GQ” simulations, we consider we have [L20phy-PM-GQ-1] and [L60geo-PM-GQ-1]. Their iteration counts are the lowest among all simulations, with 66 and 69 iterations, respectively, as shown in Table 17.

For the “PM-LQ” simulations, we consider [L20geo-PM-LQ-1] and [L60geo-PM-LQ-1]. Their iteration counts are significantly higher than those of the “PM-GQ” test cases, reaching approximately 3k and 5k iterations, respectively, according to Table 17. These values are almost
comparable to the “FPM-1v100” test cases, where [L20geo-FPM-1v100] and [L60geo-FPM-1v100] also require around 3k and 5k iterations, despite the presence of the 2D-3D coupling.

This confirms that the number of iterations of multigrid solvers is very sensitive to the mesh quality [31]. As shown in Table 17, all “FPM” meshes have some low quality elements. Even if the surface mesher respects a quality threshold, the volume mesh quality deteriorates as the number of fractures grows. As the number of fractures grows, there are many more fracture-fracture intersections, which increases the complexity of the geometry. The effects of the degradation in “FPM” mesh quality are visible in these flow simulations.

For low $K_{\text{ratio}}^{m,f}$ “FPM” simulations, we have [L20geo-FPM-1v100] and [L60geo-FPM-1v100]. According to Table 17, their number of iterations remain almost the same as their equivalent in “PM-LQ”, even with the 2D-3D coupling.

For high $K_{\text{ratio}}^{m,f}$ “FPM” simulations, we have [L20geo-FPM-heter] and [L60geo-FPM-heter]. The number of iterations is much higher than their “FPM-1v100” and “PM-LQ” counterparts due to the heterogeneous properties between the rock matrix and the fractures. Typically, according to Table 17, for the test case [L20geo-FPM-heter], PCG needs 42 times more iterations than [L20geo-FPM-1v100] in order to converge down to $10^{-14}$. The struggle of PCG-AMGCL to converge down to $10^{-14}$ is visible in figure 13a, where we observe that the residual highly oscillates starting from residual values lower than $10^{-6}$. The increase of $K_{\text{ratio}}^{m,f}$ makes the linear system (1.8) much harder to solve [54], as the number of iterations and the solver time of PCG explode.

Even worse, in test case [L60geo-FPM-heter], PCG does not even converge down to $10^{-14}$ after one week of solver time and more than 215k iterations, both of which are unfeasible in terms of time and an excessively high number of iterations for a single simulation.

To conclude, from Section 4.4, we showed that PCG-AMGCL is able to solve the linear system (1.8) for three dimensional fractured porous media with a large number of fractures.

However, the complexity of such geometries yield inevitably a poor quality 2D-3D mesh. Furthermore, in more realistic scenarios, the heterogeneity of the permeability/transmissivity might be very high. Because of these two factors, the condition number of the linear system is large (see Table 16) and the iterative linear solver struggles to converge, as seen with the oscillations in the test case [L20geo-FPM-heter], and even does not converge at all, as seen with the test case [L60geo-FPM-heter].

And even when PCG-AMGCL converges, it does so with a prohibitive number of iterations and a very long solver time. This motivates the need for an ad-hoc preconditioner for linear solvers in future work to improve both solver time and iteration count.

## 4.6 How should one choose the tolerance of the iterative solver?

In Sections 4.3, 4.4 and 4.5, we set a tolerance of $10^{-14}$ for PCG. At first glance, this might look as a very demanding tolerance. This subsection aims at justifying the choice of such a low tolerance value through the analysis of the computed velocity fields and transport simulations.

### 4.6.1 Effect on the computed velocity field

We propose to study the effect of tolerance of the iterative solver on the quality of the computed velocity field on the following test cases related to L20 geometry: [L20phy-PM-GQ-1], [L20geo-FPM-LQ-1], [L20geo-FPM-1v100] and [L20geo-FPM-heter]. Their description is given in Section 3, but we recall their main characteristics in Table 18.
<table><caption>Table 18: Main properties of [L20phy-PM-GQ-1], [L20geo-PM-LQ-1], [L20geo-FPM-1v100] and [L20geo-FPM-heter] test cases.</caption><thead><tr><th>Test case</th><th>N<sup>f</sup></th><th>N(T<sub>h</sub>)</th><th>min<sub>T</sub><br>δ(T)</th><th>N(K<sub>h</sub>)</th><th>min<sub>K</sub><br>δ(K)</th><th>K<sup>m,f</sup><sub>ratio</sub></th><th>condest<br>w/o scaling</th><th># dofs</th></tr></thead><tbody><tr><td>[L20phy-PM-GQ-1]</td><td>0</td><td>0.62M</td><td>3·10<sup>-1</sup></td><td>0</td><td>-</td><td>-</td><td>2·10<sup>5</sup></td><td>1.2M</td></tr><tr><td>[L20geo-PM-LQ-1]</td><td>0</td><td>1.9M</td><td>4·10<sup>-5</sup></td><td>0</td><td>-</td><td>-</td><td>1·10<sup>9</sup></td><td>3.9M</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>1.8k</td><td>1.9M</td><td>4·10<sup>-5</sup></td><td>0.22M</td><td>1·10<sup>-3</sup></td><td>10<sup>2</sup></td><td>2·10<sup>11</sup></td><td>4.2M</td></tr><tr><td>[L20geo-FPM-heter]</td><td>1.8k</td><td>1.9M</td><td>4·10<sup>-5</sup></td><td>0.22M</td><td>1·10<sup>-3</sup></td><td>10<sup>7</sup></td><td>7·10<sup>16</sup></td><td>4.2M</td></tr></tbody></table>

First, we simulate each test case with a direct solver, MATLAB's LU for instance, via nef-flow-fpm. The aim of this subsection is not to assess the performance of such direct solvers, whose impracticality has already been discussed in Section (4.2). However, when solving with LU, the solution of the linear system (1.8) is considered exact up to floating point errors. As a reminder, machine epsilon (or precision) is equal to $2.22 \cdot 10^{-16}$.

Thus, we obtain a reference solution for the three following unknowns (See Appendix A.1): the trace of hydraulic head in the matrix denoted $\lambda_{h,LU}^m$, the mean hydraulic head in the fracture network denoted $p_{h,LU}^f$ and the trace of hydraulic head in the fracture network denoted $\lambda_{h,LU}^f$. Thanks to the post-processing phase implemented in nef-flow-fpm, we obtain a reference solution for the other unknowns: the mean hydraulic head in the rock matrix denoted $p_{h,LU}^m$, the velocity in the rock matrix denoted $u_{h,LU}^m$ and the velocity in the fracture network denoted $u_{h,LU}^f$.

After obtaining a reference solution with the direct solver LU, each test case is then simulated via nef-flow-fpm with the iterative solver PCG-AMGCL, whose performances have been studied thoroughly in Subsections 4.4 and 4.5. For a given PCG tolerance denoted “tol”, similarly to LU, we obtain the following PCG-AMGCL solutions: $\lambda_{h,tol}^m$, $p_{h,tol}^f$, $\lambda_{h,tol}^f$, $p_{h,tol}^m$, $u_{h,tol}^m$ and $u_{h,tol}^f$.

For “PM” simulations, we only have $\lambda_{h,tol}^m$, $p_{h,tol}^m$ and $u_{h,tol}^m$ for LU and $\lambda_{h,tol}^f$, $p_{h,tol}^f$ and $u_{h,tol}^f$ for PCG-AMGCL.

With these quantities available, we compute the following relative $L^2$ errors on:

* the trace of hydraulic head in the rock matrix $\text{err}(\lambda_{h,tol}^m)$ (resp. $\text{err}(\lambda_{h,tol}^f)$ in the fracture network) ($\text{err}(\lambda_{h,tol}^m)$ for “PM” simulations is defined similarly):

$$ \text{err}(\lambda_{h,tol}^m) = \frac{\|\lambda_{h,LU}^m - \lambda_{h,tol}^m\|_{L^2(\mathcal{F}_h)}}{\|\lambda_{h,LU}^m\|_{L^2(\mathcal{F}_h)}}; \quad \text{err}(\lambda_{h,tol}^f) = \frac{\|\lambda_{h,LU}^f - \lambda_{h,tol}^f\|_{L^2(\mathcal{E}_h)}}{\|\lambda_{h,LU}^f\|_{L^2(\mathcal{E}_h)}}; \quad (4.1) $$

* the mean hydraulic head in the rock matrix $\text{err}(p_{h,tol}^m)$ (resp. $\text{err}(p_{h,tol}^f)$ in the fracture network):

$$ \text{err}(p_{h,tol}^m) = \frac{\|p_{h,LU}^m - p_{h,tol}^m\|_{L^2(\mathcal{T}_h)}}{\|p_{h,LU}^m\|_{L^2(\mathcal{T}_h)}}; \quad \text{err}(p_{h,tol}^f) = \frac{\|p_{h,LU}^f - p_{h,tol}^f\|_{L^2(\mathcal{K}_h)}}{\|p_{h,LU}^f\|_{L^2(\mathcal{K}_h)}}; \quad (4.2) $$

* the velocity in the rock matrix $\text{err}(u_{h,tol}^m)$ (resp. $\text{err}(u_{h,tol}^f)$ in the fracture network):

$$ \text{err}(u_{h,tol}^m) = \frac{\|u_{h,LU}^m - u_{h,tol}^m\|_{L^2(\mathcal{T}_h)}}{\|u_{h,LU}^m\|_{L^2(\mathcal{T}_h)}}; \quad \text{err}(u_{h,tol}^f) = \frac{\|u_{h,LU}^f - u_{h,tol}^f\|_{L^2(\mathcal{K}_h)}}{\|u_{h,LU}^f\|_{L^2(\mathcal{K}_h)}}. \quad (4.3) $$

To interpret these $L^2$ relative errors $\text{err}$, we compute $\log_{10}(\frac{1}{\text{err}})$, which represents the number of common digits between the reference solution LU and the one computed by PCG-AMGCL when $\text{err} < 1$, and is equal to zero when $\text{err} \ge 1$.

In this study, PCG-AMGCL is run with the same parameters as in Subsection 4.3. But this time, for each test case, we have as many simulations as tolerances we want to study, namely: $\{10^{-2}, 10^{-4}, 10^{-6}, 10^{-8}, 10^{-10}, 10^{-12}, 10^{-14}\}$.
The following figures display the relative errors with respect to PCG-AMGCL tolerance:

• [L20phy-PM-GQ-1] (Figure 14a) and [L20geo-PM-LQ-1] (Figure 14b): err($\lambda_{h,tol}^m$), err($p_{h,tol}^m$) and err($u_{h,tol}^m$);

• [L20geo-FPM-1v100] (Figure 14c) and [L20geo-FPM-heter] (Figure 14d): err($\lambda_{h,tol}^{m,m}$), err($\lambda_{h,tol}^f$), err($p_{h,tol}^m$), err($p_{h,tol}^f$), err($u_{h,tol}^m$) and err($u_{h,tol}^f$).

<figure><img src="image_15.png" /><figcaption>Figure 14: Relative $L^2$ errors with respect to PCG-AMGCL tolerance on the traces of hydraulic head, the mean hydraulic heads and the velocities.</figcaption></figure>

We observe that, as expected, the quality of the solution deteriorates with the increase of the tolerance of the iterative solver.

Moreover, on Figures 14b, 14c and 14d, corresponding respectively to [L20geo-PM-LQ-1], [L20geo-FPM-1v100] and [L20geo-FPM-heter], even with a tolerance of $10^{-14}$, the relative $L^2$ error on the traces of hydraulic head, on the mean hydraulic heads and on the velocities between LU and PCG-AMGCL is larger than $10^{-14}$. On Figure 14a, which corresponds to [L20phy-PM-GQ-1], where the mesh has good quality according to Table 18 and where the hydraulic conductivity in the rock matrix is equal to 1, there is no such effect.

This phenomenon is reflected in the large condition number of the linear system, as quantified in Section 3.4. Moreover, when the condition number is high, the computed solution becomes
more sensitive to rounding errors due to floating point arithmetic. For example, let us take a solver tolerance of $10^{-6}$.

In [L20phy-PM-GQ-1] with a good quality mesh (Figure 14a), the relative $L^2$ error on the velocity err($u_{h,tol}^m$) is around $10^{-5}$. In [L20geo-PM-LQ-1] with a low quality mesh (Figure 14b), err($u_{h,tol}^m$) takes a value around $10^{-4}$. This is once again due to the increase of the condition number and to the decrease in mesh quality $\delta(T)$ as seen in Table 18.

In [L20geo-FPM-1v100] (Figure 14c) we add the coupling between the flow in the rock matrix and the flow in fractured networks. The mesh has an intrinsic low quality and the condition number further increases with respect to [L20geo-PM-LQ-1], as shown in Table 18. The relative $L^2$ error err($u_{h,tol}^m$) reaches around $10^{-3}$.

In [L20geo-FPM-heter] (Figure 14d), we increase $K_{ratio}^{m,m}$, which is the contrast between the hydraulic conductivity in the rock matrix and the transmissivity in the fracture network, with respect to [L20geo-FPM-1v100]. This further increases the condition number as shown in Table 18 and the relative error err($u_{h,tol}^m$) reaches 1. This means that the computed solution with PCG-AMGCL with a tolerance of $10^{-6}$ is no longer valid for [L20geo-FPM-heter] because its norm has no common digits with the norm of LU reference solution. In general, for “FPM” simulations, especially with very high $K_{ratio}^{m,f}$ contrasts, the only way to keep as many common digits as possible with LU is to make the iterative solver converge with very low tolerances, thus motivating the choice of $10^{-14}$.

To summarize this study, the qualities of the post-processed velocities mainly depend on the condition number of the linear system. As shown in Subsection 3.4, this condition number increases as the mesh quality and contrasts in heterogeneities increase. This study helps in understanding the effect of the tolerance chosen for iterative solver, used to compute the traces of hydraulic head, on the post-processed velocities. For “FPM” simulations, the tolerance has to be chosen carefully in order to obtain a solution with satisfactory precision, especially when the velocity fields are used as input data for subsequent simulations. We confirm that this is indeed the case for the transport simulations discussed in Section 4.6.2.

### 4.6.2 Effect on advective transport

We now propose to check the effects of the iterative solver tolerance, used to compute the flow field, on transport simulations.

We choose the following test cases of increasing complexity: [Field], [L20geo-FPM-1v100] and [L50geo-FPM-1v100]. Their properties are recalled in Table 19. In particular, we recall that the [Field] test case was studied in Section 2.4.

<table><caption>Table 19: Main properties of [Field], [L20geo-FPM-1v100] and [L50geo-FPM-1v100] test cases.</caption><thead><tr><th>Test case</th><th>N<sup>f</sup></th><th>N(T<sub>h</sub>)</th><th>min<sub>T</sub><br>δ(T)</th><th>N(K<sub>h</sub>)</th><th>min<sub>K</sub><br>δ(K)</th><th>K<sub>ratio</sub><sup>m,f</sup></th><th>condest<br>w/o scaling</th><th># dofs</th></tr></thead><tbody><tr><td>[Field]</td><td>52</td><td>0.26M</td><td>5 · 10<sup>-6</sup></td><td>0.02M</td><td>4 · 10<sup>-3</sup></td><td>10<sup>2</sup></td><td>2 · 10<sup>9</sup></td><td>0.56M</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>1.8k</td><td>1.9M</td><td>4 · 10<sup>-5</sup></td><td>0.22M</td><td>1 · 10<sup>-3</sup></td><td>10<sup>2</sup></td><td>2 · 10<sup>11</sup></td><td>4.2M</td></tr><tr><td>[L50geo-FPM-1v100]</td><td>23k</td><td>5.2M</td><td>4 · 10<sup>-7</sup></td><td>1.1M</td><td>2 · 10<sup>-4</sup></td><td>10<sup>2</sup></td><td>5 · 10<sup>13</sup></td><td>12M</td></tr></tbody></table>

First, we simulate flow in the mentioned test cases thanks to nef-flow-fpm. We solve the linear system (1.8) thanks to MATLAB's LU and with PCG-AMGCL, called with the same parameters as in Section 4.3.

Figure 15 displays the convergence of PCG-AMGCL down to $10^{-14}$ for the three test cases.
<figure><img src="image_16.png" /><figcaption>Figure 15: Convergence of PCG-AMGCL with $10^{-14}$ tolerance value for test cases [Field], [L20geo-FPM-1v100] and [L50geo-FPM-1v100]: relative residual (log-scale) versus the number of iterations.</figcaption></figure>

For PCG-AMGCL, similarly to the previous Section 4.6.1, for each test case, we have as many simulations as tolerances we want to study, namely: $[10^{-2}, 10^{-4}, 10^{-6}, 10^{-8}, 10^{-10}, 10^{-12}, 10^{-14}]$.

Second, we study the post-processing results of nef-flow-fpm, in particular the flux continuity for each given tolerance.

For the rock matrix volume mesh, we only consider faces which are not in the fracture network nor in the boundary. For the surface mesh, we only consider edges which are not intersections nor in the boundary. In this scenario, fluxes on a given face (resp. edge) should be equal in absolute value but have opposite signs[65].

Nevertheless, depending on the iterative solver tolerance, there exists rock matrix faces and fracture network edges with fluxes on both sides having the same sign, as shown in Table 20. It is worth noting that for the presented test cases, this phenomenon does not happen with a direct solver.

<table><caption>Table 20: Total number and percentage of faces+edges not respecting flux continuity (local mass conservation issue) for different PCG-AMGCL tolerance values and for LU in [Field], [L20geo-FPM-1v100] and [L50geo-FPM-1v100] test cases.</caption><thead><tr><th>Test case</th><th>10<sup>-2</sup></th><th>10<sup>-4</sup></th><th>10<sup>-6</sup></th><th>10<sup>-8</sup></th><th>10<sup>-10</sup></th><th>10<sup>-12</sup></th><th>10<sup>-14</sup></th><th>LU</th></tr></thead><tbody><tr><td>[Field]</td><td>173<br>0.03%</td><td>20<br>0.0038%</td><td>1<br>0.0002%</td><td>0<br>0.0%</td><td>0<br>0.0%</td><td>0<br>0.0%</td><td>0<br>0.0%</td><td>0<br>0.0%</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>233.2k<br>6.01%</td><td>14.8k<br>0.38%</td><td>593<br>0.015%</td><td>11<br>0.0003%</td><td>0<br>0%</td><td>0<br>0%</td><td>0<br>0%</td><td>0<br>0%</td></tr><tr><td>[L50geo-FPM-1v100]</td><td>1.16M<br>11.04%</td><td>53.8k<br>0.51%</td><td>1.6k<br>0.015%</td><td>33<br>0.0003%</td><td>4<br>0.00004%</td><td>0<br>0%</td><td>0<br>0%</td><td>0<br>0%</td></tr></tbody></table>

With a high solver tolerance, we obtain that the flux continuity is no longer guaranteed for some faces and edges of the mesh. Lowering PCG tolerance alleviates the problem and the value from which the flow continuity is preserved everywhere is case-dependent. For test cases with a modest number of fractures, such as [Field], tolerances lower or equal than $10^{-8}$ are enough. However, tolerances lower or equal than $10^{-12}$ are necessary for large test cases, such as [L50geo-FPM-1v100], in order to not present “ill” faces/edges.

This is one of the motivations to choose a very low solver tolerance for “FPM” simulations
with a large number of fractures. We can even push this argument further thanks to transport simulations.

Third, we simulate Eulerian transport thanks to our code nef-transport-fpm. We recall that, as stated in Section 2.2, transport needs, among others, rock matrix and fracture network velocity fields as inputs.

We use the velocity fields computed with nef-flow-fpm. However, since there are test cases where some PCG-AMGCL tolerance values result in faces/edges that do not respect flux continuity (see Table 20), we do not perform the transport simulations on those cases. Indeed, carrying them out with these continuity mismatches would lead to non-physical simulations.

With respect to initial conditions, BCs and porosity, for test case [Field] please refer to Sub-section 2.4. In particular, a depiction can be found in Figure 3.

For test cases [L20geo-FPM-1v100] and [L50geo-FPM-1v100], a unit concentration is imposed both on the fracture network Dirichlet inflow edges and on the rock matrix Dirichlet inflow faces. Please refer to figure 9 in Subsection 3.3.1 for details. The initial condition is zero in $\Omega$. The matrix porosity $\psi_{\alpha}^{m}$ is chosen constant for each subdomain $\alpha$ and equal to 0.8. The fracture porosity $\psi_{\ell}^{f}$ is also chosen constant for each fracture $\ell$ and equal to 1. The fracture aperture $d_{\ell}^{f}$ is chosen the same for each fracture $\ell$ and equal to 0.2m.

For [L20geo-FPM-1v100], the final time is 1000 s and the time step is chosen equal to 50 s. For [L50geo-FPM-1v100], the final time is 150 s and the time step is chosen equal to 5 s. For test case [Field], as in subsection 2.4, the final time is 5000 s.

We compare the breakthrough curves given at the outflow of the rock matrix domain. Let the simulation time interval be discretized into $N_t$ points. Let $\Gamma_{\text{out}}^{D}$ be the portion of the cube or cuboid side with the output flow. The concentration at time $n \in [1, N_t]$ on face $F \in \mathcal{F}_h$ is denoted $c_F^{m,n}$. The quantity of interest is the mean output rock matrix concentration at step $n \in [1, N_t]$ at $\Gamma_{\text{out}}^{D}$ and is denoted $\bar{c}_{\text{out}}^{m,n}$ defined as:

$$ \bar{c}_{\text{out}}^{m,n} = \frac{\sum_{F \subset \Gamma_{\text{out}}^{D}} c_{F}^{m,n} |F|}{\sum_{F \subset \Gamma_{\text{out}}^{D}} |F|}, \quad \forall n \in [1, N_t]. \quad (4.4) $$

In the following, the mean output rock matrix concentrations $\bar{c}_{\text{out,LU}}^{m,n}$, computed from an input velocity field obtained with the direct method LU, are considered as a reference values. The mean output rock matrix concentrations $\bar{c}_{\text{out,tol}}^{m,n}$, computed from an input velocity field obtained with different PCG-AMGCL tolerance values “tol” are compared to the reference value $\bar{c}_{\text{out,LU}}^{m,n}$. We do so by computing the following relative error:

$$ \text{err}(c_{\text{out,tol}}^{m}) = \max_{n \in [1, N_t]} \left| \frac{\bar{c}_{\text{out,tol}}^{m,n} - \bar{c}_{\text{out,LU}}^{m,n}}{\bar{c}_{\text{out,LU}}^{m,n}} \right|. \quad (4.5) $$

The results are displayed in Table 21.

<table><caption>Table 21: Maximum over time of the relative error on the mean output porous concentration: comparison between the concentration obtained with a flow solved with LU and the one obtained for different PCG-AMGCL tolerance.</caption><thead><tr><th>Test case</th><th>10<sup>-2</sup></th><th>10<sup>-4</sup></th><th>10<sup>-6</sup></th><th>10<sup>-8</sup></th><th>10<sup>-10</sup></th><th>10<sup>-12</sup></th><th>10<sup>-14</sup></th></tr></thead><tbody><tr><td>[Field]</td><td>--</td><td>--</td><td>--</td><td>3.80e-09</td><td>1.09e-10</td><td>4.33e-12</td><td>4.72e-12</td></tr><tr><td>[L20geo-FPM-1v100]</td><td>--</td><td>--</td><td>--</td><td>--</td><td>3.93e-09</td><td>6.53e-10</td><td>4.14e-10</td></tr><tr><td>[L50geo-FPM-1v100]</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>1.49e-09</td><td>8.55e-10</td></tr></tbody></table>

The “—” symbol means that the transport simulation has not been carried out because of flux continuity mismatches.
We observe that the relative error err($c_{out,tol}^m$) is low for all tested cases, which is satisfactory. The key ingredient is to choose the iterative solver tolerance so as to guarantee local mass conservation everywhere. Furthermore, the larger the test case size, the lower the admissible tolerance becomes. The difficulty is that the iterative solver tolerance is test case dependent.

## Conclusion

In this work, we meshed 2D-3D fractured porous media, represented by the discrete fracture matrix model, with a large number of fractures, reaching an order of magnitude in the tens of thousands. We then solved time-independent single-phase flow using a mixed-hybrid finite element method. The particularity of this method is that we first solve a linear system with the mean hydraulic heads in the fracture network and the traces of hydraulic heads, both in the rock matrix and in the fracture network, as unknowns. Then, in a post-processing phase, element by element, we compute the mean hydraulic heads in the rock matrix and the velocity fields in both domains.

The major challenge is to solve efficiently the linear system, which may contain tens of millions of unknowns and which may be very ill-conditioned. Direct solvers are no longer viable as the RAM requirements become prohibitive. As the linear system matrix is sparse, square, symmetric, positive and definite, we investigated the conjugate gradient as the iterative solver, preconditioned with different techniques: incomplete Cholesky factorization and algebraic multigrid. We considered test cases of increasing geometrical complexity. The chosen test cases are known to be difficult as they may contain a large number of fractures, up to 90k, and may have very strong hydraulic conductivity/transmissivity contrasts, up to 7 orders of magnitude.

PCG-AMGCL appeared to be the most promising on the investigated black-box preconditioners. We demonstrated that the quality of the mesh and the contrasts in hydraulic conductivity are critical as they control the condition number of the linear system. We have proposed to use matrix scaling to help reduce the condition number. Nevertheless, the conditioning influences the number of iterations and the solver time. In the worst cases, we reach hundreds of thousands iterations and one week of solver time, without convergence. This showed the limits of the algebraic preconditioners. Furthermore, we assessed the sensibility of the quality of the solution with respect to the iterative solver tolerance value. We have shown that the tolerance of the iterative solver must be chosen wisely, especially if the velocity fields are used in different subsequent physical simulations such as transport. As a high solver tolerance prevents from guaranteeing flux continuity on the faces of the mesh, it might lead to non physical concentration values.

Future work could involve developing a correction procedure to restore flux continuity everywhere before running the transport simulations. Another upcoming task is the study of recently developed domain decomposition preconditioners, such as HPDDM [60, 62], or specialized multi-grid solvers for hybridized mixed methods [52], with the aim of reducing computational resource requirements.

## Acknowledgements

The authors are grateful to the CLEPS infrastructure from Inria Paris. This work was granted access to the GENCI-sponsored HPC resources of TGCC@CEA under allocation AD010607519R2. The authors also thank Pierre Jolivet for his support with PETSc, and Fractory, a joint laboratory between ITASCA Consultants SAS, CNRS, and the University of Rennes, for supplying the geometries with a large number of fractures and the heterogeneous transmissivity field.
# A Appendices

## A.1 MHFEM Discretization

Let $T \in \mathcal{T}_h$ and $F$ be a face of $T$. A basis function of $RTN^0(T)$, denoted $w_{T,F}^m$, can be written as:

$$w_{T,F}^{m}(x) = a_{T,F}^{m} x + (b_{T,F}^{m} \quad c_{T,F}^{m} \quad d_{T,F}^{m})^{T}, \quad \forall x \in T. \quad (A.1)$$

This basis function has the following properties:

$$\int_{F'} w_{T,F}^{m} \cdot n_{F'} = \delta_{F,F'}, \quad \forall F' \subset \partial T, \quad (A.2a)$$

$$\nabla \cdot w_{T,F}^{m} = \frac{1}{|T|}. \quad (A.2b)$$

Let $T \in \mathcal{T}_h$. We can decompose the velocity $u_h^m \in V_h^{\text{hyb},m}$ thanks to a linear combination of $RTN^0(T)$ basis functions:

$$u_h^m = \sum_{T \in \mathcal{T}_h} \underbrace{\sum_{F \subset \partial T} q_{T,F}^m w_{T,F}^m}_{:=u_T^m}, \quad (A.3)$$

where the $q_{T,F}^m$ coefficients physically represent the normal flux across $F \subset \partial T$.

In a similar way, we can decompose the hydraulic head $p_h^m \in P_h^m$ and the trace of hydraulic head $\lambda_h^m \in L_h^m$ thanks to a linear combination of $P_0(T)$ basis functions:

$$p_h^m = \sum_{T \in \mathcal{T}_h} p_T^m \mathbb{1}_T, \quad (A.4)$$

$$\lambda_h^m = \sum_{F \in \mathcal{F}_h} \lambda_F^m \mathbb{1}_F. \quad (A.5)$$

We define the following 4 × 4 local matrix:

$$(B_T^m)_{F,F'} = \int_T (w_{T,F}^m)^T (K^m)^{-1}_T w_{T,F'}^m, \quad \forall F' \subset \partial T. \quad (A.6)$$

Let $K \in \mathcal{K}_h$ and $E$ be an edge of $K$. A basis function of $RT^0(K)$, denoted $w_{K,E}^f$, can be written as:

$$w_{K,E}^{f}(x) = a_{K,E}^{f} x + (b_{K,E}^{f} \quad c_{K,E}^{f})^{T}, \quad \forall x \in K. \quad (A.7)$$

This basis function has the following properties:

$$\int_{E'} w_{K,E}^{m} \cdot n_{E'} = \delta_{E,E'}, \quad \forall E' \subset \partial K, \quad (A.8a)$$

$$\nabla \cdot w_{K,E}^{m} = \frac{1}{|K|}. \quad (A.8b)$$

Let $K \in \mathcal{K}_h$. We can decompose the velocity $u_h^f \in V_h^{\text{hyb},f}$ thanks to a linear combination of $RT^0(K)$ basis functions:

$$u_h^f = \sum_{K \in \mathcal{K}_h} \underbrace{\sum_{E \subset \partial K} q_{K,E}^f w_{K,E}^f}_{u_K^f}, \quad (A.9)$$

where the $q_{K,E}^f$ coefficients physically represent the normal flux across the edge $E \subset \partial K$.
In a similar way, we can decompose the hydraulic head $p_h^\mathfrak{f} \in P_h^\mathfrak{f}$ and the trace of hydraulic head $\lambda_h^\mathfrak{f} \in L_h^\mathfrak{f}$ thanks to a linear combination of $P_0(K)$ basis functions:

$$p_h^\mathfrak{f} = \sum_{K \in \mathcal{K}_h} p_K^\mathfrak{f} \mathbb{1}_K, \qquad (A.10)$$

$$\lambda_h^\mathfrak{f} = \sum_{E \in \mathcal{E}_h \setminus \mathcal{E}_h^D} \lambda_E^\mathfrak{f} \mathbb{1}_E. \qquad (A.11)$$

Let $K \in \mathcal{K}_h$ and $E \subset \partial K$. We define the following 3 × 3 local matrix:

$$(B_K^\mathfrak{f})_{E,E'} = \int_K (w_{K,E}^\mathfrak{f})^T (\mathcal{K}^\mathfrak{f})_{|K}^{-1} w_{K,E'}^\mathfrak{f}, \quad \forall E' \subset \partial K, \qquad (A.12)$$

which is symmetric, positive and definite.

## A.2 MHFEM Linear system

The matrices and vectors are defined as:

$$A^{\mathfrak{m},\mathfrak{m}} = M^{\mathfrak{m},\mathfrak{m}} - (R^{\mathfrak{m},\mathfrak{m}})^T (D^{\mathfrak{m}})^{-1} R^{\mathfrak{m},\mathfrak{m}}, \qquad (A.13a)$$

$$A^{\mathfrak{m},\mathfrak{f}} = M^{\mathfrak{m},\mathfrak{f}} - (R^{\mathfrak{m},\mathfrak{m}})^T (D^{\mathfrak{m}})^{-1} R^{\mathfrak{m},\mathfrak{f}}, \qquad (A.13b)$$

$$A^{\mathfrak{f},\mathfrak{f}} = M^{\mathfrak{f},\mathfrak{f}} - (R^{\mathfrak{m},\mathfrak{f}})^T (D^{\mathfrak{m}})^{-1} R^{\mathfrak{m},\mathfrak{f}}, \qquad (A.13c)$$

$$\mathcal{F}^{\mathfrak{m},\mathfrak{m}} = (R^{\mathfrak{m},\mathfrak{m}})^T (D^{\mathfrak{m}})^{-1} \mathcal{F}^{\mathfrak{m}}, \qquad (A.13d)$$

$$\mathcal{F}^{\mathfrak{m},\mathfrak{f}} = (R^{\mathfrak{m},\mathfrak{f}})^T (D^{\mathfrak{m}})^{-1} \mathcal{F}^{\mathfrak{m}}. \qquad (A.13e)$$

where the linear system's global blocks $D^{\mathfrak{m}}$, $D^{\mathfrak{f}}$, $R^{\mathfrak{m},\mathfrak{m}}$, $R^{\mathfrak{m},\mathfrak{f}}$, $M^{\mathfrak{m},\mathfrak{m}}$, $M^{\mathfrak{m},\mathfrak{f}}$, $M^{\mathfrak{f},\mathfrak{f}}$ and $M^{\mathfrak{f}}$ are:

$$(D^{\mathfrak{m}})_T = \sum_{F \subset \partial T} \sum_{F' \subset \partial T} ((B_T^{\mathfrak{m}})^{-1})_{F,F'} \quad \forall T \in \mathcal{T}_h, \qquad (A.14a)$$

$$(D^{\mathfrak{f}})_K = \sum_{E \subset \partial K} \sum_{E' \subset \partial K} ((B_K^{\mathfrak{f}})^{-1})_{E,E'} \quad \forall E \in \mathcal{E}_h, \qquad (A.14b)$$

$$(R^{\mathfrak{m},\mathfrak{m}})_{T,F} = \sum_{F' \subset \partial T} ((B_T^{\mathfrak{m}})^{-1})_{F,F'} \quad \forall T \in \mathcal{T}_h, \forall F \in \mathcal{F}_h^{\text{in}}, \qquad (A.14c)$$

$$(R^{\mathfrak{m},\mathfrak{f}})_{T,K} = \sum_{F \subset \partial T} ((B_T^{\mathfrak{m}})^{-1})_{K,F} \quad \forall T \in \mathcal{T}_h, \forall K \in \mathcal{K}_h, \qquad (A.14d)$$

$$(R^{\mathfrak{f}})_{K,E} = \sum_{E' \subset \partial K} ((B_K^{\mathfrak{f}})^{-1})_{E,E'} \quad \forall K \in \mathcal{K}_h, \forall E \in \mathcal{E}_h, \qquad (A.14e)$$

$$(M^{\mathfrak{m},\mathfrak{m}})_{F,F'} = \sum_{\{T \in \mathcal{T}_h | F, F' \supset \partial T\}} ((B_T^{\mathfrak{m}})^{-1})_{F,F'} \quad \forall (F, F') \in (\mathcal{F}_h^{\text{in}})^2, \qquad (A.14f)$$

$$(M^{\mathfrak{m},\mathfrak{f}})_{F,K} = \sum_{\{T \in \mathcal{T}_h | F, K \supset \partial T\}} ((B_T^{\mathfrak{m}})^{-1})_{F,K} \quad \forall F \in \mathcal{F}_h^{\text{in}}, \forall K \in \mathcal{K}_h, \qquad (A.14g)$$

$$(M^{\mathfrak{f},\mathfrak{f}})_{K,K'} = \sum_{\{T \in \mathcal{T}_h | K, K' \supset \partial T\}} ((B_T^{\mathfrak{m}})^{-1})_{K,K'} \quad \forall (K, K') \in (\mathcal{K}_h)^2, \qquad (A.14h)$$

$$(M^{\mathfrak{f}})_{E,E'} = \sum_{\{K \in \mathcal{K}_h | E, E' \supset \partial K\}} ((B_K^{\mathfrak{f}})^{-1})_{E,E'} \quad \forall (E, E') \in (\mathcal{E}_h)^2. \qquad (A.14i)$$

For the construction of such global blocks, one has to assume that there is some order, numbering and data structures linking $\mathcal{T}_h$, $\mathcal{F}_h$, $\mathcal{F}_h^{\text{in}}$, $\mathcal{K}_h$, $\mathcal{F}_h^{\text{frac}}$ and $\mathcal{E}_h$.
### A.3 AMGCL parameters

In every linear system resolution, the conjugate gradient (CG) takes two parameters: a tolerance and a maximum number of iterations. If we precondition CG with AMGCL [36, 37, 38] algebraic multigrid preconditioner, we need to give the following parameters:

* `nthreads = 1`: number of OpenMP threads (1 is sequential);

* `solver = 'cg'`: our MATLAB code `nef-flow-fpm` uses CG from AMGCL's C++ backend via a MEX file when AMGCL preconditioner is called instead of MATLAB's PCG;

* `preconditioner = 'amg'`: the chosen preconditioner [85];

* `ncycle = 1`: type of cycle; 1 is for V-cycle; 2 is for W-cycle; etc.;

* `pre_cycles = 1`: number of cycles to make as part of preconditioning;

* `maxLevels = -1`: maximum number of levels; -1 means that AMGCL's backend decides automatically;

* `coarsening = 'ruge-stuben'`: a coarsening strategy creates coarse systems in the AMG hierarchy by taking the linear system matrix at the current level, and returning a prolongation operator and its corresponding restriction operator; we choose the classic Ruge-Stuben coarsening with direct interpolation [84].

- `rs_esp_strong = 1/4`: parameter defining strong couplings;

- `rs_trunc = 'true'`: we choose prolongation operator truncation; interpolation operators, and, hence coarse operators may increase substantially towards coarser levels; without truncation, this may become too costly; truncation ignores all interpolatory connections which are smaller (in absolute value) than the largest one by a factor of `rs_esp_trunc` (see below); the remaining weights are rescaled so that the total sum remains unchanged;

- `rs_esp_trunc = 1/5`: the threshold for truncation;

* `coarse_enough = -1`: specifies when a hierarchy level is coarse enough to be solved directly; if -1 then there is none;

* `relaxation = 'spai0'`: sparse approximate inverse (SPAI) smoother [53], which computes an approximate inverse $M$ of the matrix of the linear system $A$ by minimizing the Frobenius norm of $Id - MA$, while also imposing a diagonal ("0") sparsity pattern for $M$;

* `npre = 1`: number of pre-relaxations;

* `npost = 1`: number of post-relaxations.

These parameters were chosen heuristically with small-scale single-phase flow simulations. Other parameters not mentioned are chosen automatically by AMGCL's C++ backend.

Of course, the existing options in AMGCL are much more diverse and the reader can refer to AMGCL documentation for more details.

### References

[1] P. M. Adler, J.-F. Thovert, and V. V. Mourzenko. *Fractured Porous Media*. Oxford University Press, 2012. doi:10.1093/acprof:oso/9780199666515.001.0001.

[2] A. Al-Rudaini, S. Geiger, E.-Mackay, C. Maier, and J. Pola. Comparison of Chemical-Component Transport in Naturally Fractured Reservoirs Using Dual-Porosity and Multiple-Interacting-Continua Models. *Society of Petroleum Engineers Journal*, 25(04):1964-1980, 2020. doi:10.2118/195448-PA.
[3] 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.

[4] P. R. Amestoy, A. Buttari, J.-Y. L'Excellent, and T. Mary. Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures. *ACM Transactions on Mathematical Software*, 45:2:1-2:26, 2019. doi:10.1145/3242094.

[5] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L'Excellent. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling. *SIAM Journal on Matrix Analysis and Applications*, 23(1):15-41, 2001. doi:10.1137/S0895479899358194.

[6] P. R. Amestoy, I. S. Duff, D. Ruiz, and B. Uçar. A parallel matrix scaling algorithm. In J. M. L. M. Palma, P. R. Amestoy, M. Daydé, M. Mattoso, and J. C. Lopes, editors, *High Performance Computing for Computational Science - VECPAR 2008*, pages 301-313, Berlin, Heidelberg, 2008. Springer Berlin Heidelberg. doi:10.1007/978-3-540-92859-1_27.

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

[8] P. Angot, F. Boyer, and F. Hubert. Asymptotic and numerical modelling of flows in fractured porous media. *ESAIM: Mathematical Modelling and Numerical Analysis*, 43(2):239-275, 2009. doi:10.1051/m2an/2008052.

[9] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, H. Suh, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.22, Argonne National Laboratory, 2024. doi:10.2172/2205494.

[10] S. Balay, W. D. Gropp, L. Curfman McInnes, and B. F. Smith. Efficient Management of Parallelism in Object Oriented Numerical Software Libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, *Modern Software Tools in Scientific Computing*, pages 163-202. Birkhäuser Press, 1997.

[11] S. I. Banshoya, I. Berre, and E. Keilegavlen. Simulation of Reactive Transport in Fractured Porous Media. *Transport in Porous Media*, 149(2):643-667, 2023. doi:10.1007/s11242-023-01946-0.

[12] R. Basirat, K. Goshtasbi, and M. Ahmadi. Scaling Geological Fracture Network from a Micro to a Macro Scale. *Fracture and Structural Integrity*, 14(51):71-80, 2019. doi:10.3221/IGF-ESIS.51.06.

[13] 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. *Advances in Water Resources*, 147:103759, 2021. doi:10.1016/j.advwatres.2020.103759.

[14] I. Berre, F. Doster, and E. Keilegavlen. Flow in Fractured Porous Media: A Review of Conceptual Models and Discretization Approaches. *Transport in Porous Media*, 130(1):215-236, 2019. doi:10.1007/s11242-018-1171-6.

[15] S. Berrone, S. Pieraccini, and S. Scialò. A PDE-Constrained Optimization Formulation for Discrete Fracture Network Flows. *SIAM Journal on Scientific Computing*, 35(2):B487-B510, 2013. doi:10.1137/120865884.

[16] S. Berrone, S. Pieraccini, and S. Scialò. Flow simulations in porous media with immersed intersecting fractures. *Journal of Computational Physics*, 345:768-791, 2017. doi:10.1016/j.jcp.2017.05.049.

[17] 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.

[18] I. I. Bogdanov, V. V. Mourzenko, J.-F. Thovert, and P. M. Adler. Two-phase flow through fractured porous media. *Physical Review E*, 68:026703, 2003. doi:10.1103/PhysRevE.68.026703.

[19] E. Bonnet, O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie, and B. Berkowitz. Scaling of fracture systems in geological media. *Reviews of Geophysics*, 39(3):347-383, 2001. doi:10.1029/1999RG000074.

[20] W. M. Boon, J. M. Nordbotten, and I. Yotov. Robust discretization of flow in fractured porous media. *SIAM Journal on Numerical Analysis*, 56(4):2203-2233, 2018. doi:10.1137/17M1139102.
[21] H. Borouchaki, P. Laug, and P.-L. George. Parametric surface meshing using a combined advancing-front generalized delaunay approach. *International Journal for Numerical Methods in Engineering*, 49(1-2):233-259, 2000. doi:10.1002/1097-0207(20000910/20)49:1/2<233::AID-NME931>3.0.CO;2-G.

[22] B. Bourbiaux. Fractured Reservoir Simulation: a Challenging and Rewarding Issue. *Oil and Gas Science and Technology - Rev. IFP*, 65(2):227-238, 2010. doi:10.2516/ogst/2009063.

[23] 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.

[24] 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 Journal of Numerical Analysis*, 37(3):1551-1585, 2017. doi:10.1093/imanum/drw044.

[25] 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. *Journal of Computational Physics*, 357:100-124, 2018. doi:10.1016/j.jcp.2017.12.003.

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

[27] F. Brezzi and M. Fortin. *Mixed and hybrid finite element methods*, volume 15 of Springer Series in Computational Mathematics. Springer-Verlag, New York, 1991. doi:10.1007/978-1-4612-3172-1.

[28] M. C. Cacas, E. Ledoux, G. de Marsily, B. Tillie, A. Barbreau, E. Durand, B. Feuga, and P. Peaudecerf. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. the flow model. *Water Resources Research*, 26(3):479-489, 1990. doi:10.1029/WR026i003p00479.

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

[30] R. Chen, W. Xu, Y. Chen, T. Nagel, C. Chen, Y. Hu, J. Li, and D. Zhuang. Influence of Heterogeneity on Dissolved CO2 Migration in a Fractured Reservoir. *Environmental Earth Sciences*, 83:463, 2024. doi: 10.1007/s12665-024-11753-5.

[31] Y. Chen and G. N. Wells. Multigrid on unstructured meshes with regions of low quality cells, 2024. arXiv: 2402.12947, doi:10.48550/arXiv.2402.12947.

[32] T. A. Davis. Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method. *ACM Transactions on Mathematical Software*, 30(2):196-199, 2004. doi:10.1145/992200.992206.

[33] P. Davy, R. Le Goc, and C. Darcel. A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling. *Journal of Geophysical Research: Solid Earth*, 118(4):1393-1407, 2013. doi: 10.1002/jgrb.50120.

[34] P. Davy, R. Le Goc, C. Darcel, O. Bour, J.-R. de Dreuzy, and R. Munier. A likely universal model of fracture scaling and its consequence for crustal hydromechanics. *Journal of Geophysical Research: Solid Earth*, 115(B10), 2010. doi:10.1029/2009JB007043.

[35] Jean-Raynald de Dreuzy, Géraldine Pichot, Baptiste Poirriez, and Jocelyne Erhel. Synthetic benchmark for modeling flow in 3d fractured media. *Computers and Geosciences*, 50:59-71, 2013. Benchmark problems, datasets and methodologies for the computational geosciences. doi:10.1016/j.cageo.2012.07.025.

[36] D. Demidov. AMGCL: An Efficient, Flexible, and Extensible Algebraic Multigrid Implementation. *Lobachevskii Journal of Mathematics*, 40(5):535-546, 2019. doi:10.1134/S1995080219050056.

[37] D. Demidov. AMGCL – A C++ library for efficient solution of large sparse linear systems. *Software Impacts*, 6:100037, 2020. doi:10.1016/j.simpa.2020.100037.

[38] D. Demidov, L. Mu, and B. Wang. Accelerating linear solvers for Stokes problems with C++ metaprogramming. *Journal of Computational Science*, 49:101285, 2021. doi:10.1016/j.jocs.2020.101285.

[39] P. R. B. Devloo, W. Teng, and C.-S. Zhang. Multiscale Hybrid-Mixed Finite Element Method for Flow Simulation in Fractured Porous Media. *Computer Modeling in Engineering and Sciences*, 119(1):145-163, 2019. doi:10.32604/cmes.2019.04812.
[40] A. Dundun and Y. Feng. Modeling fluid flow in fractured porous media: a comparative analysis between Darcy-Darcy model and Stokes-Brinkman model. *Journal of Petroleum Exploration and Production Technology*, 14(4):909-926, 2024. doi:10.1007/s13202-023-01743-x.

[41] 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. *Computer Methods in Applied Mechanics and Engineering*, 354:213-244, 2019. doi:10.1016/j.cma.2019.05.013.

[42] A. El-Rayes, A. Omran, M. Geriesh, and V. Hochschild. Estimation of Hydraulic Conductivity in Fractured Crystalline Aquifers Using Remote Sensing and Field Data Analyses: An Example from Wadi Nasab Area, South Sinai, Egypt. *Journal of Earth System Science*, 129:129, 2020. doi:10.1007/s12040-020-01470-5.

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

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

[45] 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.

[46] L. Formaggia, A. Fumagalli, and A. Scotti. A multi-layer reactive transport model for fractured porous media. *Mathematics in Engineering*, 4(1):1-32, 2022. doi:10.3934/mine.2022008.

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

[48] 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.V. 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.

[49] A. Fumagalli and A. Scotti. Reactive Flow in Fractured Porous Media. In R. Klöfkorn, E. Keilegavlen, F. A. Radu, and J. Fuhrmann, editors, *Finite Volumes for Complex Applications IX - Methods, Theoretical Aspects, Examples*, pages 55-73, Cham, 2020. Springer International Publishing. doi:10.1007/978-3-030-43651-3_4.

[50] G.N. 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.

[51] P.-L. George, H. Borouchaki, F. Alauzet, A. Loseille, P. Laug, and L. Maréchal. *Maillage, modélisation géométrique et simulation numérique 2*. ISTE Group, 2018. URL: https://www.istegroup.com/fr/produit/maillage-modelisation-geometrique-et-simulation-numerique-2/.

[52] J. Gopalakrishnan and S. Tan. A convergent multigrid cycle for the hybridized mixed method. *Numerical Linear Algebra with Applications*, 16(9):689-714, 2009. doi:10.1002/nla.636.

[53] M. J. Grote and T. Huckle. Parallel Preconditioning with Sparse Approximate Inverses. *SIAM Journal on Scientific Computing*, 18(3):838-853, 1997. doi:10.1137/S1064827594276552.

[54] H. Hoteit, J. Erhel, R. Mosé, B. Philippe, and Ph. Ackerer. Numerical Reliability for Mixed Methods Applied to Flow Problems in Porous Media. *Computational Geosciences*, 6:161-194, 2002. doi:10.1023/A:1019988901420.

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

[56] Y. Huang, Z. Yu, Z. Zhou, J. Wang, and Q. Guo. Modeling Flow and Solute Transport in Fractured Porous Media at Jinping I-Hydropower Station, China. *Journal of Hydrologic Engineering*, 19(9):05014007, 2014. doi:10.1061/(ASCE)HE.1943-5584.0000951.

[57] D. Huo and B. Gong. Discrete Modeling and Simulation on Potential Leakage through Fractures in CO2 Sequestration. In *SPE Annual Technical Conference and Exhibition*, pages SPE-135507-MS. SPE, 2010. doi: 10.2118/135507-MS.
[58] J. D. Hyman, C. W. Gable, S. L. Painter, and N. Makedonska. Conforming Delaunay Triangulation of Stochastically Generated Three Dimensional Discrete Fracture Networks: A Feature Rejection Algorithm for Meshing Strategy. *SIAM Journal on Scientific Computing*, 36(4):A1871–A1894, 2014. doi:10.1137/130942541.

[59] 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. *Journal of Computational Physics*, 466:111396, 2022. doi:10.1016/j.jcp.2022.111396.

[60] P. Jolivet, F. Hecht, F. Nataf, and C. Prud'homme. Scalable domain decomposition preconditioners for heterogeneous elliptic problems. In *Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis*, SC '13, New York, NY, USA, 2013. Association for Computing Machinery. doi:10.1145/2503210.2503212.

[61] P. Jolivet, M. Kern, F. Nataf, G. Pichot, and D. Zegarra 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.

[62] P. Jolivet, J.E. Roman, and S. Zampini. KSPHPDDM and PCHPDDM: Extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. *Computers and Mathematics with Applications*, 84:277–295, 2021. doi:10.1016/j.camwa.2021.01.003.

[63] M. Karimi-Fard and L. J. Durlofsky. A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features. *Advances in Water Resources*, 96:354–372, 2016. doi:10.1016/j.advwatres.2016.07.019.

[64] E. Keilegavlen, R. Berge, and A. Fumagalli et al. 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.

[65] M. Kern, G. Pichot, and D. Zegarra Vasquez. Mathematical and numerical analysis of the mixed formulation of single phase flow in three-dimensional fractured porous media. Preprint, 2025. URL: https://inria.hal.science/hal-05029638.

[66] T. Koch, D. Gläser, K. Weishaupt, S. Ackermann, M. Beck, B. Becker, S. Burbulla, H. Class, E. Coltman, S. Emmert, T. Fetzer, C. Grüninger, K. Heck, J. Hommel, T. Kurz, M. Lipp, F. Mohammadi, S. Scherrer, M. Schneider, G. Seitz, and B. Flemisch. DuMux 3 - an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. *Computers and Mathematics with Applications*, 81:423–443, 2021. doi:10.1016/j.camwa.2020.02.012.

[67] 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.

[68] B. Lepillier, A. Daniilidis, N. Gholizadeh, P.-O. Bruna, J. Kummerow, and D. Bruhn. A fracture flow permeability and stress dependency simulation applied to multi-reservoirs, multi-production scenarios analysis. *Geothermal Energy*, 7, 2019. doi:10.1186/s40517-019-0141-8.

[69] L. Li and S. H. Lee. Efficient Field-Scale Simulation of Black Oil in a Naturally Fractured Reservoir Through Discrete Fracture Networks and Homogenized Media. *SPE Reservoir Evaluation and Engineering*, 11(04):750–758, 2008. doi:10.2118/103901-PA.

[70] K. Lipnikov, G. Manzini, and M. Shashkov. Mimetic finite difference method. *Journal of Computational Physics*, 257:1163–1227, 2014. doi:10.1016/j.jcp.2013.07.031.

[71] Y.-C. Lu, Y.-M. Tien, C. H. Juang, H. Farichah, C.-J. Hsu, and V.-B. Bui. Analytical Model for Uncertainty Characterization of Fracture Intensity Measurement in Rock Masses. *Bulletin of Engineering Geology and the Environment*, 81(1):45, 2021. doi:10.1007/s10064-021-02526-6.

[72] V. Martin, J. Jaffré, and J. E. Roberts. Modeling fractures and barriers as interfaces for flow in porous media. *SIAM Journal on Scientific Computing*, 26(5):1667–1691, 2005. doi:10.1137/S1064827503429363.

[73] R. Merland, G. Caumon, B. Lévy, and P. Collon-Drouaillet. Voronoi grids conforming to 3D structural features. *Computational Geosciences*, 18:373–383, 2014. doi:10.1007/s10596-014-9408-0.

[74] J. M. Nordbotten, W. M. Boon, A. Fumagalli, and E. Keilegavlen. Unified approach to discretization of flow in fractured porous media. *Computational Geosciences*, 23(2):225–237, 2019. doi:10.1007/s10596-018-9778-9.
[75] B. Nøtinger. A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow. *Journal of Computational Physics*, 283:205-223, 2015. doi:10.1016/j.jcp.2014.11.038.

[76] M. Oda. Permeability tensor for discontinuous rock masses. *Géotechnique*, 35(4):483-495, 1985. doi:10.1680/geot.1985.35.4.483.

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

[78] A. Paluszny and R. W. Zimmerman. Numerical fracture growth modeling using smooth surface geometric deformation. *Engineering Fracture Mechanics*, 108:19-36, 2013. doi:10.1016/j.engfracmech.2013.04.012.

[79] Y.-J. Park, F. J. Cornaton, S. D. Normani, J. F. Sykes, and E. A. Sudicky. Use of Groundwater Lifetime Expectancy for the Performance Assessment of a Deep Geologic Radioactive Waste Repository: 2. Application to a Canadian Shield Environment. *Water Resources Research*, 44(4), 2008. doi:10.1029/2007WR006212.

[80] K. Pruess and T. N. Narasimhan. A Practical Method for Modeling Fluid and Heat Flow in Fractured Porous Media. *Society of Petroleum Engineers Journal*, 25(01):14-26, 1985. doi:10.2118/10509-PA.

[81] N. Rangel-Jurado, A. J. Hawkins, and P. M. Fulton. Influence of extreme fracture flow channels on the thermal performance of open-loop geothermal systems at commercial scale. *Geothermal Energy*, 11:19, 2023. doi:10.1186/s40517-023-00261-7.

[82] P. Royer, J.-L. Auriault J. Lewandowska, and C. Serres. Continuum Modelling of Contaminant Transport in Fractured Porous Media. *Transport in Porous Media*, 49(3):333-359, 2002. doi:10.1023/A:1016272700063.

[83] M. C. Saceanu, A. Paluszny, R. W. Zimmerman, and D. Mas Ivars. Fracture growth leading to mechanical spalling around deposition boreholes of an underground nuclear waste repository. *International Journal of Rock Mechanics and Mining Sciences*, 152:105038, 2022. doi:10.1016/j.ijrmms.2022.105038.

[84] K. Stüben. Algebraic multigrid (AMG). An introduction with applications. *GMD Forschungszentrum Informationstechnik*, 1999. doi:10.24406/publica-fhg-290270.

[85] U. Trottenberg, C. W. Oosterlee, and A. Schuller. Multigrid. Academic Press, 2000.

[86] B. Van Riet, S. Six, K. Walraevens, A. Vandenbohede, and T. Hermans. Assessing the Impact of Fractured Zones Imaged by ERT on Groundwater Model Prediction: A Case Study in a Chalk Aquifer in Voort (Belgium). *Frontiers in Water*, 3, 2022. doi:10.3389/frwa.2021.783983.

[87] H. S. Viswanathan, J. Ajo-Franklin, J. T. Birkholzer, J. W. Carey, Y. Guglielmi, J. D. Hyman, S. Karra, L. J. Pyrak-Nolte, H. Rajaram, G. Srinivasan, and D. M. Tartakovsky. From Fluid Flow to Coupled Processes in Fractured Rock: Recent Advances and New Frontiers. *Reviews of Geophysics*, 60(1):e2021RG000744, 2022. doi:10.1029/2021RG000744.

[88] M. Vohralík, J. Maryška, and O. Severýn. Mixed and nonconforming finite element methods on a system of polygons. *Applied Numerical Mathematics. An IMACS Journal*, 57(2):176-193, 2007. doi:10.1016/j.apnum.2006.02.005.

[89] Y. Wang, Y. Zhang, X. Yang, J. Dong, and Y. Ma. Optimizing unconventional gas extraction: The role of fracture roughness. *Physics of Fluids*, 36(3):036611, 03 2024. doiXiv:https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/5.0195281/19727189/036611\_1\_5.0195281.pdf, doi:10.1063/5.0195281.

[90] T. Wildey, S. Muralikrishnan, and T. Bui-Thanh. Unified geometric multigrid algorithm for hybridized high-order finite element methods. *SIAM Journal on Scientific Computing*, 41(5):S172-S195, 2019. doi:10.1137/18M1193505.

[91] Y.-S. Wu. Chapter 9 - multiphase flow in fractured porous media. In *Multiphase Fluid Flow in Porous and Fractured Reservoirs*, pages 207-250. Gulf Professional Publishing, 2016. doi:10.1016/B978-0-12-803848-2.00009-X.

[92] J. Xu and L. Zikatanov. Algebraic multigrid methods. *Acta Numerica*, 26:591-721, 2017. doi:10.1017/S0962492917000083.