# An efficient matrix-equivalent discrete fracture network model for variably-saturated flow in fractured porous media

Anis A. Younes

► To cite this version:

Anis A. Younes. An efficient matrix-equivalent discrete fracture network model for variably-saturated flow in fractured porous media. Journal of Computational Physics, 2025, 539, pp.114260. 10.1016/j.jcp.2025.114260. hal-05364970

HAL Id: hal-05364970

https://hal.science/hal-05364970v1

Submitted on 14 Nov 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
# An efficient Matrix-Equivalent Discrete Fracture Network model for variably-saturated flow in fractured porous media

Anis Younes¹

¹ITES, University of Strasbourg, CNRS, ENGEES, 67000, Strasbourg, France.

Submitted to Journal of Computational Physics

Contact author: Anis younes

E-mail: younes@unistra.fr
34

**Abstract**

Flow in fractured reservoirs is often controlled by a complex high-conductive fracture network. The dual-porosity model (DPM) and the discrete fracture-matrix model (DFM) are among the most widely used models for flow simulation in fractured media. The DPM gained high efficiency by replacing the fracture network by a continuum medium, but cannot be used for a disconnected network of fractures. The DFM model allows to treat explicitly the fractures but can suffer from low efficiency due the poor quality of meshes in the case of complex fracture networks. In this study, we propose an efficient matrix-equivalent discrete fracture network (MEDFN) model for unsaturated flow in fractured porous media. The main idea of the new model is to replace the matrix continuum by an equivalent low-permeable rectangular fracture network. Thus, the flow is solved on a global fracture network formed by the initial fracture network supplemented by the matrix-equivalent fracture network. The main advantages of the MEDFN model are: (i) it accurately accounts for the effect of individual fractures by explicitly considering the high-conductive fracture network, (ii) the role of the low-permeable matrix is taking into account through the added equivalent fracture network and (iii) the new model is highly efficient since the flow is solved on the set of 1D linear elements forming the global network.

The MEDFN model is developed here for the nonlinear variably-saturated flow in fractured porous media using the finite volume (FV) spatial discretization method and high-order time integration methods, via the method of lines (MOL). Numerical experiments are performed to validate the proposed MEDFN model and to assess its efficiency and accuracy, as compared to COMSOL and to the mixed finite element discrete fracture-matrix model (MFE-DFM), for the simulation of flow in fractured reservoirs under saturated and unsaturated conditions.

57

**Key words:**
Fractured porous media, unsaturated flow, discrete fracture network, discrete fracture-matrix model, finite volume, mixed finite element.

# 1. Introduction

Fractured reservoirs play an important role as they represent more than 40% of all reservoirs in the world ([1]). Fluid flow simulation in fractured reservoirs is required in several applications such as water resource management, contaminant transport, geological CO2 sequestration, geothermal energy extraction, nuclear waste management and oil and gas production on natural fractured reservoirs or on unconventional reservoirs with hydraulic fracturing techniques.

Efficient and accurate simulations of the flow in fractured reservoirs is challenging because of (i) the high contrast of permeability between the porous matrix and fractures and (ii) the complexity of the fracture network. Fractures can be empty or containing a filling material [2]. The characteristics of the fractures, such as location, size, orientation, and aperture impact significantly the flow process. Various models have been developed for the simulation of flow in fractured porous media ([3], [4]), using either the explicit or the implicit representation of the fractures [2]. One of the most widely used model based on the implicit representation of the fractures is the dual-porosity model (DPM) (e.g. [5], [6]). In this model, the fracture network is replaced by a continuum of a high permeability and the two (matrix and fracture) continua are linked by a linear exchange term. Due to its high efficiency, gained by simplifying the fracture network to a continuum media, the DPM has become a major approach for flow simulation in naturally fractured reservoirs [3]. The DPM is, however, not adapted when the flow is dominated by some discrete large fractures or when the porous medium contains disconnected fractures.

To account for the effect of individual fractures, several models have been developed using an explicit representation of the fractures. Among them, the discrete fracture network (DFN)
model considers the matrix continuum as completely impermeable with flow occurring only through the interconnected fractures ([7], [8]). This model has limitation in accuracy since it neglects the role of the matrix which serves as source/sink terms for the fracture system. To overcome this limitation, the discrete fracture-matrix model (DFM) has been developed by considering fluid flow through both, the porous matrix and the explicitly represented fractures. Assuming the aperture of the fractures small, compared to the matrix scale, DFM discretizes the fractures with elements of co-dimension one with respect to the dimension of the matrix ([9], [10]). Thus, for bi-dimensional problems, the flow processes in 2D matrix elements and in 1D fracture elements, aligned with the edges of matrix elements (conforming meshes). In the literature, the DFM model has been widely used with several numerical methods, such as Galerkin finite elements ([11], [12], [13]), vertex-centered finite volumes ([14], [15], [16]), cell-centered finite volume ([17], [18], [19]), and Mixed finite elements ([20], [21], [22], [23], [24], [25]). Hoteit and Firoozabadi [26] showed that the hybrid formulation of the mixed finite element (MFE) method, which has the traces of pressures at edges as unknowns ([27], [28]), is an excellent candidate for DFM. The hybrid MFE-DFM model assumes the matrix pressure at an edge and the pressure at the fracture, which coincides with that edge, to be the same. Thus, with MFE-DFM, no constraints are required for the size of the mesh elements near the fractures [26].

The principal limitation of the traditional DFM models is the necessity of the mesh to be conform to the geometry of the fractures which are often randomly distributed because of their geotectonic origins. As a consequence, DFM can suffer from poor efficiency because of difficulty in generating high-quality meshes when the fracture network is complex (involving a high number of fractures and/or small distances or angles between the fractures). To overcome this shortcoming, alternative methods working on non-conforming meshes have been developed, among them the Embedded Discrete Fracture Model (EDFM) ([29], [30],
[31]) where independent grids are constructed for the matrix and fracture domains. The fractures are then coupled with the matrix by modifying the transmissibility factors to reflect the geometry and properties of the fractures [32]. The matrix-fracture transfer is assumed to be proportional to the average normal distance between the fracture cell and the matrix cell.

The EDFM has been extended to the projection-based EDFM (pEDFM) ([33], [34], [35]) to deal with low-conductive fractures. In pEDFM, additional matrix-fracture and fracture-fracture connections are added by projecting the embedded fractures to the interfaces of matrix cells. Other non-conforming methods have been developed in the literature for flow in fractured reservoirs like the non-conforming finite element method based on Lagrange multipliers ([36], [37]), the extended finite element model ([38], [39]) which enriches the finite element space by adding basis functions designed to capture the discontinuity, and the reinterpreted discrete fracture model where the permeability tensor is represented by a one-dimensional line Dirac-δ functions contained in the permeability tensor ([4], [40], [41]).

The non-conforming mesh methods allow to avoid mesh problems of conforming mesh methods but have some disadvantages, such as a large amount of computation on the transmissibility factors and difficulty of implementation and therefore low efficiency for complex fracture networks involving dense and narrow intersecting fractures.

In this work, we propose a new model which takes into account both flow in the explicit fracture network and in the low-permeability matrix. The new model can efficiently solve the nonlinear unsaturated flow problem on complex fracture networks involving a high number of fractures with small distances or angles between the fractures. The key idea of the new model, named matrix-equivalent discrete fracture network (MEDFN) model, is to define a fictitious rectangular fracture network that imitates the matrix continuum. Hence, instead of replacing the fracture network by a continuum media (as with the DPM), we replace the matrix continuum by an equivalent rectangular low-permeable fracture network. The flow problem is
then solved on the global network formed by the initial fracture network supplemented by the additional matrix-equivalent fracture network. The main advantages of the MEDFN model are (i) the effect of individual fractures can be accounted for as the fractures are explicitly taken into account, (ii) the mesh problems encountered with DFM are avoided since the problem is solved only on the global network, (iii) the role of the matrix (which acts as a sink/source term for fractures) is taken into account (contrarily to the standard DFN model) through the equivalent low-permeable fracture network and (iv) the model is highly efficient since the nonlinear flow problem is solved only on a set of 1D elements forming the global fracture network.

The paper is structured as following. In section 2, we present the mathematical equations ruling the unsaturated flow in porous media. In section 3, we describe the MEDFN model and detail the employed spatial and temporal discretization methods. In section 4, numerical experiments are performed to investigate the validity and efficiency of the new MEDFN model by comparison against COMSOL and the MFE-DFM model for the simulation of flow in fractured porous media under saturated and unsaturated conditions. Some conclusions are given in the last section of the article.

## 2. Governing Equations

In this work, we consider unsaturated flow in 2D fractured reservoirs where the fluid flows through both the high-permeable fracture network and the low-permeable matrix which is assumed isotropic. The fractures are assumed infilled by a porous medium whose properties are different from those of the matrix. The fluid mass conservation equation in variably-saturated porous media can be written as following:

$$ \frac{\partial \theta}{\partial t} + \nabla \cdot \mathbf{q} = r \frac{\partial H}{\partial t} + \nabla \cdot \mathbf{q} = q_s \quad (1) $$
156 And the Darcy-Buckingham law writes:

$$157 \qquad q = -k_r K \nabla H \qquad (2)$$

158 where $r = c(h) + S_s \frac{\theta}{\theta_s}$ in which $c(h) = \partial\theta/\partial h$ is the specific moisture capacity [L⁻¹], $S_s$ is the specific mass storativity related to head changes [L⁻¹], $\theta$ is the current water content [L³L⁻³], $\theta_s$ is the saturated water content [L³L⁻³], $H = h + y$ is the water head [L], $h$ is the pressure head [L], $y$ is the upward vertical coordinate [L], $t$ is the time [T], $q$ is the Darcy velocity [LT⁻¹], $k_r$ is the relative conductivity [-], $K = \frac{\rho g}{\mu} k$ is the hydraulic conductivity [LT⁻¹], $\rho$ is the fluid density [ML⁻³], $g$ is the gravity acceleration [LT⁻²], $\mu$ is the fluid dynamic viscosity [ML⁻¹T⁻¹], $q_s$ is the source or sink term per unit area [T⁻¹] and $k$ is the permeability [L²].

166 The water content is related to the pressure head following the standard van Genuchten [42] model:

$$168 \qquad S_e = \frac{\theta(h) - \theta_r}{\theta_s - \theta_r} = \begin{cases} \frac{1}{(1 + |\alpha h|^n)^m} & h < 0 \\ 1 & h \ge 0 \end{cases} \qquad (3)$$

169 where $\alpha$ [L⁻¹] and $n$ [-] are the van Genuchten parameters, $m = 1 - 1/n$, $S_e$ is the effective saturation [-], and $\theta_r$ is the residual water content [L³L⁻³].

171 The conductivity is related to the saturation following the Mualem [43] model,

$$172 \qquad k_r = S_e^{1/2} \left[ 1 - \left( 1 - S_e^{1/m} \right)^m \right]^2 \qquad (4)$$
173

The nonlinear system (1)-(4) has to be solved numerically for a fractured reservoir formed by a low-permeable matrix crossed by a complex high-conductive fracture network.

174

175

### 3. The MEDFN model

176

Noticing that the flow in fractured reservoirs is mainly dominated by the fracture network, the idea of the MEDFN model is to discretize explicitly the high-conductive fracture network (instead of replacing it by a continuum as with the DPM) and to replace the low-permeable matrix (which serves as a sink/source for the fractures) by an additional equivalent low-permeable rectangular fracture network. The flow problem can then be solved on a global fracture network, formed by the combination of the initial fracture network with the matrix-equivalent low-permeable fracture network. This yields a highly efficient model, especially in the case of a complex fracture network. For instance, in the simple case of two crossing fractures with a small angle between them (see Figure 1), the DFM model requires a conforming mesh with a high number of elements (1220 in the Figure 1a) involving a high ratio of the largest to the smallest element (area ratio around 200). This can significantly impact the performance of the numerical model which is known to be tied to the shape and number of elements. The configuration shown in Figure 1a can be much more problematic in the case of a complex fracture network formed by a large number of fractures.

177

178

179

180

181

182

183

184

185

186

187

188

189

190

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

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

191
**Figure 1:** Description of the DFM (a) and MEDFN (b) models in the case of two crossing fractures with a small angle between them.

With the MEDFN model, the matrix continuum is replaced by an equivalent rectangular network (the blue fracture network in the Figure 1b) which is added to the original fracture network (in red in the Figure 1b). The flow is then solved on the global network, formed by both the red and blue networks of the Figure 1b. The MEDFN model can be highly efficient since (i) it requires a limited number of unknowns compared to the DFM model and (ii) it solves the unsaturated flow problem on the set of 1D linear elements forming the global fracture network.

In this work, the MEDFN model is combined with the method of lines (MOL) for the solution of the nonlinear unsaturated flow problem. With the MOL, we discretize first the spatial derivatives and then we integrate in time the semi-discrete problem with high order methods, as a system of Ordinary Differential Equations (ODEs). The MOL was shown to be effective for the solution of the Richards' equation ([44], [45], [46]).

### 3.1 Spatial discretization

The finite volume (FV) method is well-known for its robustness and cost-effectiveness in discretization of conservation laws. FV is chosen here for the spatial discretization of the unsaturated flow with the MEDFN model on the global fracture network.

Attributing the sink/source term to the intersection nodes, the integration of the mass conservation Eq. (1) over the fracture element $k$ yields:

$$r^k \ell^k b^k \underset{t}{\stackrel{\sim}{\overline{H}}}^k + q_1^k + q_2^k = 0 \quad (5)$$
where $\ell$ and $b^k$ are respectively the length and the aperture of the fracture $k$,
$r^k = c^k + S_s \frac{\theta^k}{\theta_s^k}$, with $\theta_s^k$, $c^k$ and $\theta^k$, respectively, the saturated water content, the specific
moisture capacity and the current water content at the fracture $k$, $\overline{H}^k$ is the mean head at $k$,
$q_1^k$ and $q_2^k$ are the outer water fluxes at the nodes 1 and 2 of the fracture $k$ (Figure 2a).

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Flow notations for (a) the fracture $k$ and (b) the intersection of fractures at the node $i$.</figcaption></figure>

The flux $q_i^k$ is obtained using a FV discretization of the Darcy Eq. (2):

$$q_i^k = k_r^k K^k \times \frac{2b^k}{\ell^k} \begin{pmatrix} \overline{H}^k & H_i^k \end{pmatrix} \quad (6)$$

where $k_r^k$ is the relative conductivity of the fracture $k$, $K^k$ is the hydraulic conductivity of fracture element $k$ and $H_i^k$ is the head at the node $i$ of $k$.

The mass conservation at the intersection node $i$ shared by $n$ fractures (see Figure 2b corresponding to $n = 5$), writes:

$$\sum_{l=1}^{n} q_{i}^{l} = Q_{i}^{*} \qquad (7)$$
228 where $Q_i^*$ corresponds to the boundary flux or the sink/source term at the intersection node $i$.

229 Substituting Eq. (6) into Eq. (7) gives

$$ \sum_{l=1}^{n} \left( k_r^l K^l \times \frac{2b^l}{\ell^l} (\overline{H^l} \quad H_i^l) \right) = Q_i^* \qquad (8) $$

231 Denoting by $H_i$ the head at the intersection node $i$ and using the continuity of the head at the node $i$ ($H_i = H_i^l, \forall l = 1, \dots, n$), the Eq. (8) gives

$$ H_i = \frac{\sum_{l=1}^{n} k_r^l K^l \frac{2b^l}{\ell^l} \overline{H^l} Q_i^*}{\sum_{l=1}^{n} k_r^l K^l \frac{2b^l}{\ell^l}} \qquad (9) $$

234 Assuming the node 1 of the fracture $k$ shared by $n_1$ fractures and the node 2 by $n_2$ fractures, the combination of Eq. (5), Eq. (6) and Eq. (9) yields

$$ \begin{aligned} & r^k \ell^k b^k \underset{t}{\stackrel{\sim}{\overline{H^k}}} + k_r^k K^k \frac{2b^k}{\ell^k} \left( \sum_{l=1}^{n_1} k_r^l K^l \frac{2b^l}{\ell^l} \overline{H^l} \right) k_r^k K^k \frac{b^k}{\ell^k} H^k \left( \sum_{l=1}^{n_2} k_r^l K^l \frac{2b^l}{\ell^l} \overline{H^l} \right) \\ & = k_r^k K^k \frac{2b^k}{\ell^k} \left[ \sum_{l=1}^{n_1} k_r^l K^l \frac{2b^l}{\ell^l} + \frac{Q_1^*}{\sum_{l=1}^{n_2} k_r^l K^l \frac{2b^l}{\ell^l}} \right] \end{aligned} \qquad (10) $$

237 This equation is written for all the fracture branches of the global network to form the final system to solve. The unknowns are the mean head $\overline{H}$ at the fractures. The number of unknowns is therefore the total number of fracture branches in the global network. The system of Eq. (10) is highly nonlinear since $r^k = r^k(\overline{H}^k)$ and $K^l = K^l(\overline{H}^l)$ are highly nonlinear functions obtained from Eq. (3) and Eq. (4).
## 3.2 The fracture aperture $b^k$ of the matrix-equivalent network

In the absence of high-conductive fractures, the flow in the unfractured domain is modeled with the MEDFN model using the equivalent rectangular fracture network. For steady-state flow, the value of $b^k$ has no effect on the head distribution, but does affect the total flow rate through the matrix. For the steady-state regime ($r^k = 0$), the Eqs. (5) and (6) can be combined to express the mean head at the fracture $k$ as the average of the heads at the two extremities ($\overline{H}^k = \frac{1}{2}(H_1^k + H_2^k)$). Substituting this relation into Eq. (6) and using Eq. (7) yields a global system with the head at the intersection nodes of the rectangular network ($H_i$) as unknowns. The resulting system can easily be shown to be algebraically equivalent to the well-known five-points finite difference method by choosing the aperture $b^k$ of the fracture $k$ to be the sum of the width of its two adjacent rectangles. Therefore, in this work, the aperture $b^k$ of the matrix equivalent network is chosen as follows. For an interior fracture $k$ of the equivalent rectangular fracture network, $b^k$ is the average of the distance to its two parallel adjacent fractures. For a boundary fracture $k$, $b^k$ is the mid-distance to its first interior adjacent fracture. This choice of $b^k$ insures the equivalence of the MEDFN model with the 2D finite difference method in the case of steady state flow in unfractured domains. It is worth noticing that in the case of transient regime, the system of Eq. (12) cannot be simplified. Eq. (12) should be solved for the head at all interior and boundary fractures (except those corresponding to Dirichlet boundary condition) as unknowns. Hence, for transient flow in unfractured domains, there is no equivalence between the MEDFN model and the finite difference method.

## 3.3 Temporal discretization
High-order integration methods allows to obtain an accurate solution of the pressure-head form of Richards' equation with good mass balance properties [47]. The system of Eq. (10) is integrated in time with the sophisticated DASPK time solver. With this solver, the order of the time integration as well as the time step size are adapted during the simulation to reduce the computational burden while maintaining a fixed accuracy.

The ODE system of Eq. (10) is written in the following form:

$$G(t, y, y') = 0 \qquad (11)$$

where the vector $y = [\overline{H}^k]_{k=1,...,nb\_fractures}$ corresponds to the unknown head at all the fractures of the global network.

This system is solved with the DASPK solver, based on the Fixed Leading Coefficient Backward Difference Formulas (FLCBDF) which has good stability properties. The nonlinear system is solved with the Newton method. To improve efficiency, the Jacobian matrix $J$ is calculated by perturbation using the column grouping technique and the same $J$ is employed for many time steps (it is updated only when the convergence fails). The linearized system is solved with the preconditioned Krylov iterative method which is well adapted for large linear systems.

During the simulation, DASPK adapts the order of the time integration (from 1 to 5) as well as the time step size in order to reduce the computational effort while maintaining the relative and absolute time errors under small values ($10^{-7}$ in this work).

## 4. Numerical experiments

In this section, numerical experiments are performed to validate the new MEDFN model and to assess its performances. The results of the MEDFN model are compared against those of COMSOL as well as those of the MFE-2D [46] and MFE-DFM [24] models. The MFE-2D
model is a 2D-2D model based on the MFE method for the solution of Richards' equation ([48], [28]) where 2D triangular elements are employed for both the spatial discretization of the matrix and of the fractures. The MFE-2D model employs a highly unstructured mesh with local mesh refinement inside and around the fractures in order to correctly discretize the flow process in the fracture branches with small apertures. The MFE-DFM is a 1D-2D model using conforming meshes where the 1D fracture elements are aligned with the interfaces (edges) between the 2D matrix elements. The MFE method is used for the spatial discretization of the unsaturated flow in both the 1D-fracture elements and the 2D matrix elements. With MFE-DFM, the matrix-fracture exchanges are implicitly taken into account without using any transfer function, by assuming the same pressure for the matrix element edge and that of the 1D fracture element which coincides with that edge [24]. Both MFE-DFM and MFE-2D models use high-order time integration methods and efficient automatic time-stepping schemes via the MOL ([46], [24]). Different test problems are investigated hereafter. All simulations are carried out on an Intel Xeno E-2246G, CPU 3.60 GHz PC with 32.00 GB of RAM.

## 4.1 Steady-state and transient flow in an unfractured domain

This first problem deals with steady state and transient flow in the low permeable matrix and aims to validate the developed approach in the absence of high-conductive fractures. The isotropic matrix has a square shape (10m×10m) and a permeability $k_m = 10^{-10} m^2$.

In a first stage, we investigate steady-state flow with the MEDFN model. The boundary conditions for the steady-state flow problem are as following. The head is fixed to $H_L = 1m$ in the left side ($x=0$) for ($1m \le y \le 3m$) and to $H_{up} = 0m$ on the top ($y=10m$) for ($7m \le x \le 9m$). All the rest of the boundaries of the domain are impermeable (no-flow boundary condition). The problem is simulated with the MEDFN model using a coarse
311 (25×25) and a fine (50×50) rectangular fracture networks and the results are compared to a
312 2D reference solution obtained with COMSOL using a fine mesh formed by 6290 triangles.
313 Figure 3 shows the coarse rectangular network and the pressure distribution in the domain
314 calculated with the MEDFN model.

<figure><img src="image_5.png" /><figcaption>316 **Figure 3:** Rectangular (25×25) network and steady-state head distribution obtained with the
317 MEDFN model.</figcaption></figure>

318 The results of the MEDFN model using the coarse and the fine networks are compared to the
319 COMSOL reference solution for the x component of the Darcy velocity $q_x$ along the diagonal
320 $y = -x + 10$ and for the head profile along $x = 5m$. The results of Figure 4 show that the
321 MEDFN model yields identical head and velocity profiles with the coarse and the fine
322 networks. Furthermore, the MEDFN head and velocity profiles are almost similar to the
323 reference COMSOL solution. Moreover, the total flow rate crossing the domain calculated
324 with MEDFN ($4.15 \ m^2/s$) is very close to that obtained with the 2D COMSOL reference
325 solution ($4.16 \ m^2/s$). These results validate the MEDFN model for steady-state flow in
326 unfractured domains.

327
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Head profile along $x = 5m$ and Darcy velocity $q_x$ profile along $y = -x + 10$ with COMSOL and the MEDFN model using the coarse (25×25) and the fine (50×50) networks.</figcaption></figure>

In a second stage, we investigate the transient flow in the unfractured domain with the MEDFN model. The boundary conditions for the transient problem are as following. The head is fixed to $H_{up,R} = 100m$ on the upper ($y=10m$) and right ($x=10$) boundaries. All the rest of the boundaries of the domain are impermeable (no-flow boundary condition). The storage coefficient is $S = 10^{-5}$ [Pa$^{-1}$] and the initial head in the domain is $H_{init} = 100m$. A pumping well is located at the lower left corner of the domain with a pumping rate $Q_p = 10^{-2}$ [m$^3$/s]. The problem is simulated with the MEDFN model using the coarse (25×25) network and the results are compared to a 2D reference COMSOL solution obtained using a fine mesh formed by 6342 triangles.

Due to pumping, the results of simulation show a transient drawdown at the well location. The head distribution in the domain at $t = 6000s$ is plotted in the Figure 5.
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Head distribution at t = 6000s obtained with the MEDFN model using the (25×25) network.</figcaption></figure>

The head profile along the first diagonal (y = x) at t = 6000s as well as the evolution of the head at the pumping well during time, obtained with both the MEDFN and COMSOL models are shown in Figure 6. A very good agreement can be observed between the head profiles at t = 6000s obtained with COMSOL and with the MEDFN model using the coarse mesh (Figure 6a). Further, the head drawdown at the well, calculated with the two models is almost the same (Figure 6b). These results validate the MEDFN model for transient flow in unfractured domains.

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Head distribution along (y = x) at t = 6000s (a) and head evolution at the pumping well with COMSOL and the MEDFN model (b).</figcaption></figure>
## 4.2 Flow in a saturated 2D domain with a horizontal fracture

In this first problem, we consider a horizontal square shape domain (9m×9m) formed by a saturated homogeneous matrix of permeability $k_m = 10^{-10} m^2$ crossed in its middle by a highly conductive horizontal fracture of 8m length. The fracture has an aperture of 0.04m and a permeability $k_f = 10^{-5} m^2$. The main flow occurs from the left side to the right side. The boundary conditions are as follows. The head is fixed to $H_L = 1m$ in the left side ($x=0$) for ($4m \le y \le 5m$) and it is fixed to $H_R = 0m$ on all the right side ($x=9m$). All the rest of the boundaries of the domain are impermeable (no-flow boundary condition).

The domain is discretized with an unstructured triangular mesh formed by 10 339 triangular elements (see Figure 7a) for the MFE-2D model with an element area ratio of the largest to the smallest elements of 42. The MFE-DFM model uses an almost uniform mesh formed by 6376 2D triangular elements with an element area ratio of 2.5 and 45 1D fracture elements (Figure 7b). For the MEDFN model, the matrix is replaced by a rectangular fracture network formed by 46 horizontal fractures crossing 46 vertical fractures of 9m length. The global network for the MEDFN model is formed by 4221 fracture branches (Figure 7c).

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

371
<figure><img src="image_10.png" /><figcaption>**Figure 7:** Flow in a 2D domain with a horizontal fracture: Mesh for MFE-2D (a), MFE-DFM (b) and MEDFN (c) models and head distribution for MFE-2D (d), MFE-DFM (e) and MEDFN (f) models.</figcaption></figure>

The flow in the fractured domain is highly affected by the presence of the fracture. Indeed, in the absence of the horizontal fracture, the total flow rate from the left side to the right side is around $6.3 \times 10^{-4} m^2/s$. When the horizontal fracture is present, the total flow rate increases to $1.34 \times 10^{-3} m^2/s$ and around 80% of the fluid flows through the fracture. The steady-state head distribution obtained with the three models are plotted in the Figure 7. This figure shows that the MEDFN model (Figure 7f) yields almost the same head distribution than the MFE-2D (Figure 7d) and the MFE-DFM (Figure 7e) models. These results validate the MEDFN model in the case of flow in saturated fractured porous media.

### 4.3 Flow in a homogeneous 2D domain with a ‘+’-shaped barrier/fracture network

In this test case, we consider a 2D problem where a $(9m \times 9m)$ square domain is crossed in its middle by a ‘+’-shaped barrier/fracture network. As for the previous test case, a no-flow boundary condition is applied at the top and bottom. The left side is impermeable, except between $y = 4m$ and $y = 5m$, where the head is fixed to $H_L = 1m$. The right side has a fixed
head $H_R = 0$ m. The rock matrix has a permeability $k_m = 10^{-10} m^2$. The horizontal fracture and the vertical barrier have 8m length and 0.04m aperture. The permeability of the fracture is $k_f = 10^{-5} m^2$ and that of the barrier is $k_b = 10^{-15} m^2$.

<figure><img src="image_11.png" /><figcaption>**Figure 8:** Flow in a 2D homogeneous domain with a ‘+’-shaped barrier/fracture network: Mesh for MFE-2D (a), MFE-DFM (b) and MEDFN (c) models and head distribution for MFE-2D (d), MFE-DFM (e) and MEDFN (f) models.</figcaption></figure>

The MFE-2D model is used with an unstructured triangular mesh formed by 13753 elements with local refinement inside and near the fracture and the barrier (Figure 8a). The MFE-DFM is used with a mesh formed by 6594 triangles (Figure 8b). As previously, the MEDFN model replaces the matrix by a rectangular fracture network formed by 46 horizontal fractures and 46 vertical fractures of 9m length. The combination of this rectangular network with the ‘+’-
shaped barrier/fracture network yields a global network for the MEDFN model formed by 403
404 4304 linear elements (Figure 8c).

Due to the barrier, a very slight decrease of the total flow rate, crossing the system from the 405
left side to the right side, is observed compared to the previous test case. Indeed, the total flow 406
rate is now 1.31×10⁻³ m² / s with around 91.4% of the fluid flowing through the fracture. 407

Therefore, the proportion of the fluid passing through the fracture becomes more significant 408
than previously because of the effect of the barrier. The results of simulation with the three 409
models are shown in the Figure 8. Similar head distributions are obtained with the MFE-2D 410
(Figure 8d), the MFE-DFM (Figure 8e) and the MEDFN (Figure 8f) models. These results 411
validate the new model for the simulation of homogeneous saturated domains involving 412
fractures and barriers. 413

### 4.4 Flow in a heterogeneous 2D domain with a '+'-shaped barrier/fracture network

In this test case, the matrix of permeability $k_m^1 = 10^{-10} m^2$, contains a (3m×3m) square shaped 415
zone with a very low permeability $k_m^2 = 10^{-12} m^2$ located downward to the barrier. The 416
horizontal fracture crosses both the barrier and the very low permeability zone (Figure 9a). 417
The three models are used with the same meshes than for the previous homogeneous case. For 418
the MFE-2D and MFE-DFM models, the 2D triangular elements located inside the 419
heterogeneity zone are now attributed the permeability $k_m^2$. For the MEDFN model, all the 420
1D-fracture branches of the rectangular network representing the matrix that are located 421
inside the very low permeability zone are now attributed the permeability $k_m^2$. 422

The results of simulations show that the heterogeneity affects the head distribution (Figure 9) 423
and induces a reduction of the total flow rate crossing the system, which is now 424
1.16×10⁻³ m² / s. The part of the fluid passing through the fracture also reduces to 84%. The 425
head distribution obtained with the MFE-2D (Figure 9a), the MFE-DFM (Figure 9b) and 426
427

MEDFN (Figure 9c) models are in very good agreement which validates the MEDFN model
for the simulation of saturated flow in heterogeneous domains involving barriers and
fractures.

428

429

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

430

<caption>**Figure 9:** Flow in a 2D heterogeneous domain with ‘+’-shaped barrier/fracture network:
Head distribution with the MFE-2D (a), the MFE-DFM (b) and the MEDFN (c) models.</caption>

431

432

433

4.5 Flow in a saturated 2D domain with a regular fracture network

434

435

436

437

438

This test case has been studied in [6] and has been suggested as a benchmark for single-phase
flow in fractured porous media in [10]. The domain is a 1m square with an inlet flux
$q_N = -1m^2/s$ at the left side and a prescribed head $H_R = 1m$ at the right side. The top and
bottom sides are impermeable. The domain is crossed by the regular fracture network shown
in the Figure 10.

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

439
440

**Figure 10:** The 2D saturated flow problem with a regular fracture network and the COMSOL reference solution.

441

442

443

444

445

446

447

448

449

450

451

452

453

The permeability of the matrix is $10^{-7} m^2$. The regular high-conductive fracture network, shown in the Figure 10, has a uniform aperture $b=10^{-3} m$ and a permeability $10^{-3} m^2$. The reference solution for this problem is obtained with COMSOL using a very fine 2D/2D triangular mesh based on 800825 elements where both the matrix and the fractures are discretized with 2D triangular elements. With this fine mesh, each fracture is discretized in its normal direction using 5 triangles and the mesh becomes coarser away from the fractures. The obtained COMSOL reference solution is depicted in the Figure 10.

Four different uniform orthogonal meshes containing 8×8, 14×14, 26×26 and 50×50 low-permeable fractures of 1m length (see Figure 11) are used to represent the matrix with the MEDFN model. The permeability of the equivalent fracture network is $10^{-7} m^2$. The number of low and high permeable fractures for each mesh are given in Table 1. The number of unknowns with the MEDFN model is the total number of fractures (see Table 1).
454

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

455

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

456

<caption>Figure 11: Four different orthogonal fracture networks used with the MEDFN model.</caption>

457

458 The results of the MEDFN model using the 26×26 orthogonal grid are compared to the
459 COMSOL reference solution for the head distribution along the vertical fracture at x = 0.5m
460 in the Figure 12. The results of this figure show a perfect fitting of the head obtained by the
461 MEDFN model and by the COMSOL reference solution along the vertical fracture x = 0.5m.
<figure><img src="image_16.png" /><figcaption>**Figure 12:** Head distribution along the fracture located at $x = 0.5m$ with COMSOL and MEDFN models.</figcaption></figure>

In addition, the head error inside the domain is calculated for the four meshes of the Figure 11 using $err_H = \sqrt{\frac{1}{Nb_N} \sum_{i=1}^{Nb_N} \frac{(H_i^{ref} - H_i)^2}{\Delta H_{ref}^2}}$ where $H_i^{ref}$ and $H_i$ are the head at the intersection node $i$ obtained respectively from the COMSOL reference solution and the MEDFN model and $\Delta H_{ref} = \max_{\Omega} H_{ref} - \min_{\Omega} H_{ref}$.

The head error for the four meshes are given in Table 1. The results in this table show that the head error decreases as the mesh is refined, with and order of convergence of around 1.8 for the investigated test case.

<table><thead><tr><th>mesh</th><th>Low-perm. fractures</th><th>High-perm. fractures</th><th>Total fractures</th><th>Length of fractures</th><th>err<sub>H</sub></th><th>Order of convergence</th></tr></thead><tbody><tr><td>8×8</td><td>140</td><td>40</td><td>180</td><td>0.1429</td><td>3.5 10<sup>-3</sup></td><td></td></tr><tr><td>14×14</td><td>376</td><td>13</td><td>778</td><td>0.0769</td><td>1.14 10<sup>-3</sup></td><td>1.8</td></tr><tr><td>26×26</td><td>1322</td><td>23</td><td>2690</td><td>0.04</td><td>3.4 10<sup>-4</sup></td><td>1.8</td></tr><tr><td>50×50</td><td>4180</td><td>41</td><td>8442</td><td>0.0204</td><td>1.1 10<sup>-4</sup></td><td>1.7</td></tr></tbody></table>
472

**Table 1:** Head errors with the MEDFN model using different orthogonal grids.

473

4.6 Flow in a saturated 2D domain with a random distribution of fractures

474

475

476

477

478

479

480

481

In this problem, the horizontal square (9m×9m) domain is formed by a matrix of low permeability $k_m = 10^{-13} m^2$ and a complex network of high-conductive fractures of permeability $k_f = 10^{-10} m^2$. The fracture network is generated using a random distribution of fracture size and Fisher orientation distribution (Figure 13). The conductive fractures have an aperture of 0.04m. The flow boundary conditions are as follows: all boundaries are impermeable, except at the left side ($x=0m$) for ($4.2m \le y \le 4.8m$) where the head is fixed to $H_L = 1m$ and at the right side ($x=9m$) for ($4.2m \le y \le 4.8m$) where the head is fixed to $H_R = 0m$ (Figure 13).

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

482

483

**Figure 13:** Description of the saturated problem with a random distribution of fractures.

484

485

486

487

488

Because of the complexity of the fracture network, the construction of a 2D-2D mesh for the MFE-2D model is not possible and the results of the MEDFN model are only compared to those of the MFE-DFM model. Besides, due to the small distances and small angles between some fractures, the conforming mesh for the MFE-DFM model (Figure 14a) is highly unstructured and involves a large number of elements (11164 triangles). The quality of the
mesh is also deteriorating with a very high area ratio of the largest to the smallest element (area ratio 845300). For the MEDFN model, the mesh is formed by a global fracture network obtained by adding the complex high-conductive fracture network to the low-permeable rectangular fracture network representing the matrix continuum. As previously, the matrix equivalent fracture network is formed by 46 horizontal fractures crossing 46 vertical fractures of 9m length. The obtained global fracture network is formed by 6501 fracture branches (Figure 14b). The results of the simulation with the MFE-DFM and MEDFN models are plotted in the Figure 14c and the Figure 14d, respectively. A very good agreement can be observed between the head distribution obtained by the two models (Figure 14). These results demonstrate the validity of the new MEDFN model for the simulation of flow in fractured saturated porous media involving a complex fracture network.

<figure><img src="image_18.png" /></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 14:** The saturated problem with a random distribution of fractures. Mesh for the MFE-DFM (a) and MEDFN (b) models and head distribution with the MFE-DFM (c) and MEDFN (d) models.</figcaption></figure>

## 4.7 Water infiltration in a fractured initially dry soil

The simulation of water infiltration into dry soils is challenging since several convergence issues can be encountered due to the presence of sharp wetting fronts [49]. To investigate the robustness and efficiency of the developed model for such situations, we simulate the problem of infiltration into an initially dry soil containing a fracture network (Figure 15). The fractured porous medium is a (9m×9m) vertical square where water is injected at the part ($4m \le x \le 6m$) of the surface ($y=9m$) with a constant flow rate of $5 \times 10^{-6}$ m/s. The head is fixed to $H = 0m$ at the left part ($0 \le x \le 2m$) of the bottom ($y=0$). All the rest of the boundaries are impermeable (Figure 15). Initial conditions correspond to a very dry soil with a hydrostatic pressure distribution inside the domain corresponding to a head $H = 0m$. Two configurations are investigated. In the first configuration, the domain contains a relative sparse distribution of 30 fractures, generated using random size between 2m and 8m and
Fisher orientation distribution with an angle of 45° and a dispersion coefficient of 1.4 (Figure 15a). In the second configuration, the domain contains a dense fracture distribution with a fracture network formed by 80 fractures (Figure 15b), generated using the same size and orientation distributions then for the first configuration. All fractures have an aperture of 4cm. The unsaturated hydraulic parameters of the matrix and fractures are depicted in Table 2.

<table><caption>**Table 2:** Hydraulic parameters for the problem of water infiltration into a fractured dry soil.</caption><thead><tr><th></th><th>K<sub>s</sub> (cm/s)</th><th>S<sub>s</sub> (1/cm)</th><th>θ<sub>s</sub></th><th>θ<sub>r</sub></th><th>α (1/cm)</th><th>n</th></tr></thead><tbody><tr><td>Matrix</td><td>10<sup>-7</sup></td><td>10<sup>-8</sup></td><td>0.4</td><td>0.1</td><td>0.017</td><td>1.7</td></tr><tr><td>Fracture</td><td>10<sup>-4</sup></td><td>10<sup>-8</sup></td><td>0.4</td><td>0.1</td><td>0.017</td><td>1.7</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure 15:** Boundary conditions, sparse fracture network configuration (a) and dense fracture network configuration (b) for the problem of infiltration into a fractured dry soil.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 16:** Infiltration into a dry soil involving a sparse fracture distribution. Mesh for the MFE-DFM (a) and MEDFN (b) models and head distribution for the MFE-DFM (c) and MEDFN (d) models.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 17:** Infiltration into a dry soil involving a dense fracture distribution. Mesh for the MFE-DFM (a) and MEDFN (b) models and head distribution for the MFE-DFM (c) and MEDFN (d) models.</figcaption></figure>

The two models are run over the simulation time of 10⁵s for the sparse and the dense fracture configurations. For the sparse fracture distribution configuration, the MFE-DFM model is used with a conforming mesh formed by 11343 triangular elements (Figure 16a) with an element area ratio (between the maximum and minimum area) of 550. The MEDFN model uses 5857 fracture branches. The obtained saturation distributions after 10⁵s with the MEDFN (Figure 16d) and MFE-DFM (Figure 16c) models are in very good agreement. These results demonstrate the validity of the new MEDFN model for the simulation of flow in fractured
unsaturated porous media. Concerning the efficiency of the MEDFN, in comparison to MFE-DFM, the MEDFN model requires only 64s for the whole simulation, whereas the MFE-DFM model requires 475s. Thus, for the sparse fracture distribution configuration, the new MEDFN model is around 7 times more efficient than the MFE-DFM model. Due to the complexity of the fracture network for the dense fracture distribution configuration, the MFE-DFM model requires a highly unstructured mesh (Figure 17a) with a large number of elements (59669 triangles) and a very large area ratio (≈ 2×10⁻⁶) between the largest and smallest area elements. The MEDFN model uses a global mesh formed by 8441 fracture branches. The saturation distributions after 10⁶s, obtained with the two models are plotted in the Figure 17. This figure shows that the results of the MEDFN model (Figure 17b) are in very good agreement with those of the MFE-DFM model (Figure 17a). Concerning efficiency of the two models, the MFE-DFM model requires a large CPU time of 2435s for the whole simulation because of the degradation of the quality of the mesh and the large number of elements, while the MEDFN model requires only 95s for the whole simulation. Thus in this case, the MEDFN model is around 25 times more efficient than the MFE-DFM model. This highlights the great benefit of MEDFN model to simulate flow in unsaturated fractured porous media, especially in the case of complex dense fracture network.

562
564 **5. Conclusion**

565 The simulation of unsaturated flow in fractured reservoirs is a challenging task due to the
566 nonlinearity of the Richards' equation, the high contrast of permeability between the matrix
567 and the fractures and the complexity of the fracture network. The DPM is a widely used
568 model for practical flow simulations in fractured reservoirs, but DPM cannot account for the
569 effect of individual fractures. The DFM model avoids this shortcoming by using an explicit
570 representation of fractures. DFM has been combined with several numerical methods, among
571 which the hybrid formulation of the MFE is the most suitable. The MFE-DFM model assumes
572 the same pressure for the matrix element edge and that of the fracture branch which coincides
573 with that edge. The main limitation of the MFE-DFM model is related to the generation of
574 conforming meshes which can be of poor quality in the case of complex fracture networks
575 involving small distances and/or small angles between the fractures.

576 To avoid these limitations, the MEDFN model is proposed in this work for modeling
577 unsaturated flow in fractured media. The main idea of MEDFN is as follows: instead of
578 replacing the fracture network (which controls the flow in the reservoir) by a continuum (as
579 with the DPM), we replace the matrix (which serves as a sink/source for the fractures) by an
580 equivalent rectangular fracture network. The flow is then solved on a global fracture network
581 formed by the combination of the original fracture network and the added low-permeable
582 fracture network representing the matrix. This yields a very efficient model since the flow
583 problem is only solved on the set of 1D elements forming the global fracture network. The
584 new MEDFN model is combined here with the MOL for the solution of the nonlinear
585 Richards' equation. Numerical experiments are performed to investigate the validity and
586 performance of the MEDFN model for flow simulation in fractured reservoirs under saturated
587 and unsaturated conditions. The results of simulation of saturated flow in unfractured media
show that the flow in the matrix can be accurately represented by the flow across the equivalent fracture network in both steady-state and transient regimes. For saturated fractured media, the MEDFN model yields similar results than the MFE-DFM model for heterogeneous domains involving barriers and fractures as well as for domains involving a complex fracture network. Furthermore, the MEDFN model yields very small errors of the head distribution inside the domain and along fractures when compared to a reference COMSOL solution obtained using a very fine discretization.

For unsaturated conditions, numerical experiments have been performed to simulate water infiltration in a fractured dry soil in the case of sparse fracture distribution as well as in the case of dense fracture distribution. For both situations, the results of the MEDFN model are in very good agreement with those of the MFE-DFM model. Concerning the efficiency of the MEDFN in comparison to MFE-DFM, the developed MEDFN model is found to be 7 times more efficient than the MFE-DFM model for the sparse fracture distribution. For the dense fracture distribution, the MEDFN model is around 25 times more efficient than the MFE-DFM model. These results highlight the great benefit of the new MEDFN model for the simulation of variably-saturated flow in fractured porous media, especially in the case of complex fracture networks.

The MEDFN model has been developed in this work for 2D isotropic fractured domains. MEDFN can be extended to 3D problems by representing the fractures, which are generally 2D planar interfaces, by a fracture network formed by a set of 1D high-conductive linear elements and representing the matrix by an equivalent 3D rectangular fracture network of 1D low-permeable elements. The extension of the MEDFN model to anisotropic media is more challenging and is currently under investigation.
REFERENCES

613

614

615

616

617

618

[1] B. Liang, H. Jiang, J. Li, C. Gong, A systematic study of fracture parameters effect on fracture network permeability based on discrete-fracture model employing Finite Element Analyses, Journal of Natural Gas Science and Engineering 28 (2016) 711-722. https://doi.org/10.1016/j.jngse.2015.12.011.

619

620

621

[2] I. Berre, F. Doster, E. Keilegavlen, Flow in fractured porous media: A review of conceptual models and discretization approaches, (2018). https://doi.org/10.48550/arXiv.1805.05701.

622

623

624

[3] Z. Xu, Y. Yang, The hybrid dimensional representation of permeability tensor: a reinterpretation of the discrete fracture model and its extension on nonconforming meshes, (2021). https://doi.org/10.48550/arXiv.2101.12358.

625

626

627

[4] G. Fu, Y. Yang, A hybrid-mixed finite element method for single-phase Darcy flow in fractured porous media, Advances in Water Resources 161 (2022) 104129. https://doi.org/10.1016/j.advwatres.2022.104129.

628

629

630

[5] G.I. Barenblatt, Iu.P. Zheltov, I.N. Kochina, Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata], Journal of Applied Mathematics and Mechanics 24 (1960) 1286-1303. https://doi.org/10.1016/0021-8928(60)90107-6.

631

632

633

[6] S.. Geiger, M.. Dentz, I.. Neuweiler, A Novel Multirate Dual-Porosity Model for Improved Simulation of Fractured and Multiporosity Reservoirs, SPE Journal 18 (2013) 670-684. https://doi.org/10.2118/148130-PA.

634

635

636

637

[7] G. Aldrich, J.D. Hyman, S. Karra, C.W. Gable, N. Makedonska, H. Viswanathan, J. Woodring, B. Hamann, Analysis and Visualization of Discrete Fracture Networks Using a Flow Topology Graph, IEEE Transactions on Visualization and Computer Graphics 23 (2017) 1896-1909. https://doi.org/10.1109/TVCG.2016.2582174.

638

639

640

641

[8] S. Tabrizinejadas, A. Younes, H. Hoteit, J. Carrayrou, M. Fahs, Coupled fluid flow, solute transport and dissolution processes in discrete fracture networks: An advanced Discontinuous Galerkin model, Advances in Water Resources 180 (2023) 104540. https://doi.org/10.1016/j.advwatres.2023.104540.

642

643

644

[9] V. Martin, J. Jaffré, J.E. Roberts, Modeling Fractures and Barriers as Interfaces for Flow in Porous Media, SIAM J. Sci. Comput. 26 (2005) 1667-1691. https://doi.org/10.1137/S1064827503429363.

645

646

647

[10] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson, A. Tatomir, Benchmarks for single-phase flow in fractured porous media, Advances in Water Resources 111 (2018) 239-258. https://doi.org/10.1016/j.advwatres.2017.10.036.

648

649

650

[11] M. Karimi-Fard, A. Firoozabadi, Numerical Simulation of Water Injection in 2D Fractured Media Using Discrete-Fracture Model, in: OnePetro, 2001. https://doi.org/10.2118/71615-MS.

651

652

653

654

[12] N. Zhang, J. Yao, Z. Huang, Y. Wang, Accurate multiscale finite element method for numerical simulation of two-phase flow in fractured media using discrete-fracture model, Journal of Computational Physics 242 (2013) 420-438. https://doi.org/10.1016/j.jcp.2012.12.006.

655

656

657

[13] J.G. Kim, M.D. Deo, Comparison of the Performance of a Discrete Fracture Multiphase Model With Those Using Conventional Methods, in: OnePetro, 1999. https://doi.org/10.2118/51928-MS.
658

659

660

[14] J.E.P. Monteagudo, A. Firoozabadi, Control-Volume Model for Simulation of Water Injection in Fractured Media: Incorporating Matrix Heterogeneity and Reservoir Wettability Effects, SPE Journal 12 (2007) 355-366. https://doi.org/10.2118/98108-PA.

661

662

663

[15] V. Reichenberger, H. Jakobs, P. Bastian, R. Helmig, A mixed-dimensional finite volume method for two-phase flow in fractured porous media, Advances in Water Resources 29 (2006) 1020-1036. https://doi.org/10.1016/j.advwatres.2005.09.001.

664

665

666

667

[16] R. Zhang, L. Zhang, J. Luo, Z. Yang, M. Xu, Numerical simulation of water flooding in natural fractured reservoirs based on control volume finite element method, Journal of Petroleum Science and Engineering 146 (2016) 1211-1225. https://doi.org/10.1016/j.petrol.2016.08.024.

668

669

670

[17] T.H. Sandve, I. Berre, J.M. Nordbotten, An efficient multi-point flux approximation method for Discrete Fracture-Matrix simulations, Journal of Computational Physics 231 (2012) 3784-3800. https://doi.org/10.1016/j.jcp.2012.01.023.

671

672

673

674

[18] R. Ahmed, M.G. Edwards, S. Lamine, B.A.H. Huisman, M. Pal, Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model, Journal of Computational Physics 284 (2015) 462-489. https://doi.org/10.1016/j.jcp.2014.12.047.

675

676

677

678

[19] W. Fang, C. Liu, J. Li, H. Jiang, J. Pu, H. Gu, X. Qin, A discrete modeling framework for reservoirs with complex fractured media: Theory, validation and case studies, Journal of Petroleum Science and Engineering 170 (2018) 945-957. https://doi.org/10.1016/j.petrol.2017.11.050.

679

680

681

[20] H. Hoteit, A. Firoozabadi, Multicomponent fluid flow by discontinuous Galerkin and mixed methods in unfractured and fractured media, Water Resources Research 41 (2005). https://doi.org/10.1029/2005WR004339.

682

683

684

[21] J. Moortgat, A. Firoozabadi, Higher-order compositional modeling of three-phase flow in 3D fractured porous media based on cross-flow equilibrium, Journal of Computational Physics 250 (2013) 425-445. https://doi.org/10.1016/j.jcp.2013.05.009.

685

686

687

[22] A. Zidane, A. Firoozabadi, An efficient numerical model for multicomponent compressible flow in fractured porous media, Advances in Water Resources 74 (2014) 127-147. https://doi.org/10.1016/j.advwatres.2014.08.010.

688

689

690

691

[23] B. Koohbor, M. Fahs, H. Hoteit, J. Doummar, A. Younes, B. Belfort, An advanced discrete fracture model for variably saturated flow in fractured porous media, Advances in Water Resources 140 (2020) 103602. https://doi.org/10.1016/j.advwatres.2020.103602.

692

693

694

[24] A. Younes, H. Hoteit, R. Helmig, M. Fahs, A robust fully mixed finite element model for flow and transport in unsaturated fractured porous media, Advances in Water Resources 166 (2022) 104259. https://doi.org/10.1016/j.advwatres.2022.104259.

695

696

697

[25] A. Younes, B. Koohbor, M. Fahs, H. Hoteit, An efficient discontinuous Galerkin - mixed finite element model for variable density flow in fractured porous media, Journal of Computational Physics 477 (2023) 111937. https://doi.org/10.1016/j.jcp.2023.111937.

698

699

700

[26] H. Hoteit, A. Firoozabadi, An efficient numerical model for incompressible two-phase flow in fractured media, Advances in Water Resources 31 (2008) 891-905. https://doi.org/10.1016/j.advwatres.2008.02.004.
701 [27] Chavent, G., Jaffré, J., Mathematical Models and Finite Elements for Reservoir Simulation, 1986. https://shop.elsevier.com/books/mathematical-models-and-finite-elements-for-reservoir-simulation/chavent/978-0-444-70099-5 (accessed July 25, 2024).

704 [28] A. Younes, P. Ackerer, F. Delay, Mixed finite elements for solving 2-D diffusion-type equations, Reviews of Geophysics 48 (2010). https://doi.org/10.1029/2008RG000277.

706 [29] A. Moinfar, A. Varavei, K. Sepehrnoori, R.T. Johns, Development of an Efficient Embedded Discrete Fracture Model for 3D Compositional Reservoir Simulation in Fractured Reservoirs, SPE Journal 19 (2013) 289-303. https://doi.org/10.2118/154246-PA.

710 [30] X. Yan, Z. Huang, J. Yao, Y. Li, D. Fan, An efficient embedded discrete fracture model based on mimetic finite difference method, Journal of Petroleum Science and Engineering 145 (2016) 11-21. https://doi.org/10.1016/j.petrol.2016.03.013.

713 [31] Y. Xu, K. Sepehrnoori, Development of an Embedded Discrete Fracture Model for Field-Scale Reservoir Simulation With Complex Corner-Point Grids, SPE Journal 24 (2019) 1552-1575. https://doi.org/10.2118/195572-PA.

716 [32] J.S. de Araujo Cavalcante Filho, M. Shakiba, A. Moinfar, K. Sepehrnoori, Implementation of a Preprocessor for Embedded Discrete Fracture Modeling in an IMPEC Compositional Reservoir Simulator, in: OnePetro, 2015. https://doi.org/10.2118/173289-MS.

720 [33] M. Tene, S.B.M. Bosma, M.S. Al Kobaisi, H. Hajibeygi, Projection-based Embedded Discrete Fracture Model (pEDFM), Advances in Water Resources 105 (2017) 205-216. https://doi.org/10.1016/j.advwatres.2017.05.009.

723 [34] J. Jiang, R.M. Younis, An improved projection-based embedded discrete fracture model (pEDFM) for multiphase flow in fractured reservoirs, Advances in Water Resources 109 (2017) 267-289. https://doi.org/10.1016/j.advwatres.2017.09.017.

726 [35] X. Rao, X. He, K. Du, H. Kwak, A. Yousef, H. Hoteit, A Novel Projection-based Embedded Discrete Fracture Model (pEDFM) for Anisotropic Two-phase Flow Simulation Using Hybrid of Two-point Flux Approximation and Mimetic Finite Difference (TPFA-MFD) Methods, Journal of Computational Physics 499 (2024) 112736. https://doi.org/10.1016/j.jcp.2023.112736.

731 [36] M. Köppel, V. Martin, J. Jaffré, J.E. Roberts, A Lagrange multiplier method for a discrete fracture model for flow in porous media, Comput Geosci 23 (2019) 239-253. https://doi.org/10.1007/s10596-018-9779-8.

734 [37] P. Schädle, P. Zulian, D. Vogler, S.R. Bhopalam, M.G.C. Nestola, A. Ebigbo, R. Krause, M.O. Saar, 3D non-conforming mesh model for flow in fractured porous media using Lagrange multipliers, Computers & Geosciences 132 (2019) 42-55. https://doi.org/10.1016/j.cageo.2019.06.014.

738 [38] A. Fumagalli, A. Scotti, An Efficient XFEM Approximation of Darcy Flows in Arbitrarily Fractured Porous Media, OGSAT 69 (2014) 555-564. https://doi.org/10.2516/ogst/2013192.

741 [39] B. Flemisch, A. Fumagalli, A. Scotti, A Review of the XFEM-Based Approximation of Flow in Fractured Porous Media, in: G. Ventura, E. Benvenuti (Eds.), Advances in Discretization Methods: Discontinuities, Virtual Elements, Fictitious Domain Methods, Springer International Publishing, Cham, 2016: pp. 47-76. https://doi.org/10.1007/978-3-319-41246-7_3.
746

747

748

749

[40] H. Guo, W. Feng, Z. Xu, Y. Yang, Conservative numerical methods for the reinterpreted discrete fracture model on non-conforming meshes and their applications in contaminant transportation in fractured porous media, Advances in Water Resources 153 (2021) 103951. https://doi.org/10.1016/j.advwatres.2021.103951.

750

751

752

753

[41] Z. Xu, Z. Huang, Y. Yang, The hybrid-dimensional Darcy's law: a non-conforming reinterpreted discrete fracture model (RDFM) for single-phase flow in fractured media, Journal of Computational Physics 473 (2023) 111749. https://doi.org/10.1016/j.jcp.2022.111749.

754

755

756

[42] M.Th. van Genuchten, A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal 44 (1980) 892-898. https://doi.org/10.2136/sssaj1980.03615995004400050002x.

757

758

759

[43] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resources Research 12 (1976) 513-522. https://doi.org/10.1029/WR012i003p00513.

760

761

762

763

[44] M. Fahs, A. Younes, F. Lehmann, An easy and efficient combination of the Mixed Finite Element Method and the Method of Lines for the resolution of Richards' Equation, Environmental Modelling & Software 24 (2009) 1122-1126. https://doi.org/10.1016/j.envsoft.2009.02.010.

764

765

766

[45] M.W. Farthing, F.L. Ogden, Numerical Solution of Richards' Equation: A Review of Advances and Challenges, Soil Science Society of America Journal 81 (2017) 1257-1269. https://doi.org/10.2136/sssaj2017.02.0058.

767

768

769

770

[46] A. Younes, B. Koohbor, B. Belfort, P. Ackerer, J. Doummar, M. Fahs, Modeling variable-density flow in saturated-unsaturated porous media: An advanced numerical model, Advances in Water Resources 159 (2022) 104077. https://doi.org/10.1016/j.advwatres.2021.104077.

771

772

773

[47] M.D. Tocci, C.T. Kelley, C.T. Miller, Accurate and economical solution of the pressure-head form of Richards' equation by the method of lines, Advances in Water Resources 20 (1997) 1-14. https://doi.org/10.1016/S0309-1708(96)00008-5.

774

775

776

[48] B. Belfort, F. Ramasomanana, A. Younes, F. Lehmann, An Efficient Lumped Mixed Hybrid Finite Element Formulation for Variably Saturated Groundwater Flow, Vadose Zone Journal 8 (2009) 352-362. https://doi.org/10.2136/vzj2008.0108.

777

778

779

[49] Y. Zha, J. Yang, L. Yin, Y. Zhang, W. Zeng, L. Shi, A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil, Journal of Hydrology 551 (2017) 56-69. https://doi.org/10.1016/j.jhydrol.2017.05.053.

780

781

782