# A Discrete Fracture Network Model for Coupled Variable-Density Flow and Dissolution with Dynamic Fracture Aperture Evolution

Anis A. Younes, Husam Musa Baalousha, Lamia Guellouz, Marwan Fahs

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

## To cite this version:

Anis A. Younes, Husam Musa Baalousha, Lamia Guellouz, Marwan Fahs. A Discrete Fracture Network Model for Coupled Variable-Density Flow and Dissolution with Dynamic Fracture Aperture Evolution. Water, 2025, 17 (13), pp.1904. 10.3390/w17131904. hal-05364965

HAL Id: hal-05364965

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

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Discrete Fracture Network Model for Coupled Variable-Density Flow and Dissolution with Dynamic Fracture Aperture Evolution

Anis Younes <fnref n="1" />*, Husam Musa Baalousha <fnref n="2" />, Lamia Guellouz <fnref n="3" /> and Marwan Fahs <fnref n="1" />,*

## Abstract

Fluid flow and mass transfer processes in some fractured aquifers are negligible in the low-permeability rock matrix and occur mainly in the fracture network. In this work, we consider coupled variable-density flow (VDF) and mass transport with dissolution in discrete fracture networks (DFNs). These three processes are ruled by nonlinear and strongly coupled partial differential equations (PDEs) due to the (i) density variation induced by concentration and (ii) fracture aperture evolution induced by dissolution. In this study, we develop an efficient model to solve the resulting system of nonlinear PDEs. The new model leverages the method of lines (MOL) to combine the robust finite volume (FV) method for spatial discretization with a high-order method for temporal discretization. A suitable upwind scheme is used on the fracture network to eliminate spurious oscillations in the advection-dominated case. The time step size and the order of the time integration are adapted during simulations to reduce the computational burden while preserving accuracy. The developed VDF-DFN model is validated by simulating saltwater intrusion and dissolution in a coastal fractured aquifer. The results of the VDF-DFN model, in the case of a dense fracture network, show excellent agreement with the Henry semi-analytical solution for saltwater intrusion and dissolution in a coastal aquifer. The VDF-DFN model is then employed to investigate coupled flow, mass transfer and dissolution for an injection/extraction well pair problem. This test problem enables an exploration of how dissolution influences the evolution of the fracture aperture, considering both constant and variable dissolution rates.

**Keywords:** discrete fracture network; variable-density flow; dissolution; aperture evolution; finite volume; MOL

## 1. Introduction

In this work, we consider coupled fluid flow, mass transfer and dissolution processes with variable density in fractured aquifers. Variable-density flow in fractured domains

Academic Editor: Chin H Wu

Received: 21 April 2025

Revised: 2 June 2025

Accepted: 23 June 2025

Published: date

**Citation:** Younes, A.; Baalousha, H.M.; Guellouz, L.; Fahs, M. A Discrete Fracture Network Model for Coupled Variable-Density Flow and Dissolution with Dynamic Fracture Aperture Evolution. *Water* **2025**, *17*, x.

https://doi.org/10.3390/xxxxx

**Copyright:** © 2025 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Institut Terre et Environnement de Strasbourg, Université de Strasbourg, Centre National de la Recherche Scientifique (CNRS), Ecole Nationale du Génie et de l'Eau et de l'Environnement de Strasbourg (ENGEES), UMR 7063, 67084 Strasbourg, France</footnote>

<footnote n="2">² Department of Geosciences, College of Petroleum Engineering and Geosciences, King Fahd University of Petroleum and Minerals (KFUPM), Dhahran 31261, Saudi Arabia; baalousha@web.de</footnote>

<footnote n="3">³ Laboratoire Modélisation en Hydraulique et Environnement, Ecole Nationale d'Ingénieur de Tunis, Université Tunis El Manar, Le Belvédère, Tunis 1002, Tunisia; lamia.guellouz@enit.utm.tn</footnote>

<footnote>* Correspondence: younes@unistra.fr (A.Y.); fahs@unistra.fr (M.F.)</footnote>
finds applications in various instances including CO₂ storage [1,2], geothermal energy [3,4] and groundwater management in fractured aquifers [5–9]. Furthermore, when dis-solution occurs, it significantly impacts the transmissivity of the fractures, subsequently influencing the flow and the solute transport phenomena. Variable-density flow with dis-solution has been shown to play a role in various processes encountered in fractured me-dia, such as geothermal energy extraction [10], the management of oil reservoirs [11], CO₂ sequestration [12] and nuclear waste disposal [13]. Further, coupled processes of variable-density flow and dissolution in fractured aquifers can occur as a result of excessive pump-ing, creating voids and sinkholes that lead to land subsidence or ground collapse. This phenomenon is observed across the world and has great economic and social impacts as it can endanger the safety of human lives and properties [14].

Modeling coupled fluid flow, mass transfer and dissolution with variable density in fractured aquifers is challenging due to the high degree of nonlinearities induced by (i) the density variations and (ii) the evolution of the fracture aperture due to dissolution [15]. In the literature, the focus on dissolution processes has predominantly centered on unfractured domains and their role in fostering the creation and propagation of fractures [16,17]. The effect of dissolution on fracture evolution in fractured domains has received limited attention [18,19].

Either explicit or implicit representations of fractures can be used for modeling flow and mass transfer in fractured domains [20]. In the implicit approach, the fractures are taken into account using equivalent single- or multi-porosity models [4]. These models are based on equivalent parameters, which depend on the characteristics of the fractures as well as on the fracture network configuration [21]. The explicit approach, for its part, uses a discrete fracture representation. Among the variants of the explicit approach, the discrete fracture matrix (DFM) model allows for considering flow and mass transfer in both the fracture network and the rock matrix [22–24]. Fractures are represented by thin regions with distinct properties from those of the rock matrix [5]. The DFM model requires a mesh conforming to the fracture network; therefore, because of gridding problems, it is not adapted for situations with a high density of fractures. Recently, methods have emerged using non-conforming meshes to overcome gridding-related challenges [25,26].

Another variant of the explicit approach is the discrete fracture network (DFN) model, suitable for situations where matrix permeability is negligible. In this model, the matrix is considered impermeable, with flow and transport occurring exclusively through interconnected fractures [27,28]. DFN models are particularly well-suited for porous media where the permeability and porosity are entirely derived from fractures [20]. Recently, Tabrizinejadas et al. [29] showed the advantages of DFNs for modelling flow and transport in fractured domains without considering density variations. To the best of our knowledge, coupled variable-density flow and dissolution processes within a DFN framework remain unexplored. This paper aims to fill this gap by developing an efficient numerical model for this highly nonlinear problem to better investigate the effect of density variation and dissolution processes in DFNs.

In this study, the coupled nonlinear equations governing variable-density flow, solute transport and dissolution are simultaneously solved to prevent any potential errors associated with operator splitting. The numerical model integrates appropriate numerical schemes for both space and time discretization. The method of lines (MOL) is utilized to transform partial differential equations (PDEs) into a set of ordinary differential equations (ODEs) by discretizing spatial derivatives while maintaining time derivatives in their continuous form. The spatial discretization is based on the finite volume (FV) method, recognized for its robustness and cost-effectiveness in handling conservation laws [30]. An upwind scheme is used at the fracture intersection nodes to avoid unphysical oscillations
observed with high advection transport. Temporal discretization of the obtained ODE system is performed with the DASPK [31] time solver, leveraging high-order integration methods and an efficient adaptive time-stepping scheme. DASPK has demonstrated efficiency in solving coupled flow and mass transport equations in both unfractured [32,33] and fractured domains [34]. The newly developed model for variable-density flow and dissolution with a DFN is referred to as VDF-DFN.

This paper is structured as follows. Section 2 presents the mathematical equations ruling the coupled processes of variable-density flow, solute transfer and dissolution. Section 3 provides a comprehensive explanation of the novel VDF-DFN numerical model. Section 4 is devoted to numerical experiments where the model’s validity is first investigated against the semi-analytical solution for a variant of the Henry problem involving saltwater intrusion and dissolution in a fractured coastal aquifer. Then, we apply the VDF-DFN model to simulate coupled fluid flow, mass transfer and dissolution for a horizontal injection/extraction well pair problem. Additionally, we investigate the evolution of the fracture aperture under both constant and variable dissolution rates.

## 2. The Mathematical Model

The mathematical model is constructed under the following assumptions. The fracture network is of codimension one within a two-dimensional bounded domain. Flow and transport in the fractures are one-dimensional, and transport is governed by advection, dispersion and dissolution processes. The fractures are filled with a porous medium through which the flow is governed by Darcy’s law. The fracture permeability is a function of the aperture, and the porosity is constant. The temperature is constant, and mechanical deformations are negligible. The fracture aperture varies solely due to dissolution, which is governed by a first-order dissolution approximation. The fluid density varies linearly with the concentration.

The flow inside a fracture with a variable aperture $b[L]$ due to dissolution is governed by the following mass conservation equation:

$$ \frac{\partial(\phi b)}{\partial t} + \nabla \cdot (b\mathbf{q}) = 0 \qquad (1) $$

where the fluid velocity $\mathbf{q}[L.T^{-1}]$ is obtained using Darcy’s law, which is widely valid for flow in porous media [35,36]. In the case of variable density, it is written as

$$ \mathbf{q} = -\frac{\rho_0 g \kappa}{\mu} \left( \nabla h + \frac{\rho - \rho_0}{\rho_0} \nabla z \right) \qquad (2) $$

where $\phi$ is the porosity [-] of the porous medium in the fracture, $t$ is the time $[T]$, $h = \frac{P}{\rho_0 g} + z$ is the pressure head $[L]$, $P$ is the pressure $[Pa]$, $z$ is the elevation $[L]$, $g$ is the gravity acceleration $[L.T^{-2}]$, $\rho$ is the density of mixed saltwater and freshwater $[M.L^{-3}]$, $\mu$ is the fluid viscosity $[M.L^{-1}.T^{-1}]$, $\rho_0$ is the freshwater density $[M.L^{-3}]$ and $\kappa$ is the permeability $[L^2]$ approximated using the cubic law $\kappa = \frac{b^2}{12}$ [37].

The transport of a contaminant inside a fracture is governed by the following reactive transport equation [38]:
$$ \frac{\partial(\phi bc)}{\partial t} + \nabla \cdot (b\mathbf{q}c) - \nabla \cdot (bD\nabla c) = R(c) \quad (3) $$

where c is the solute concentration $[M.L^{-3}]$, R(c) is the source/sink term $[M.L^{-2}.T^{-1}]$ representing the dissolution process and D $[L^{2}.T^{-1}]$ is the dispersion coefficient in the fracture, given by

$$ D = \alpha_L |\mathbf{q}| + D_m \quad (4) $$

where $D_m$ $[L^2.T^{-1}]$ is the diffusion coefficient and $\alpha_L$ $[L]$ is the longitudinal dispersivity.

Using mass conservation Equation (1), transport Equation (3) becomes

$$ \phi b \frac{\partial c}{\partial t} + b \mathbf{q} \cdot \nabla c - \nabla \cdot (b D \nabla c) = R(c) \quad (5) $$

The flow system (Equations (1) and (2)) and transport Equation (5) are coupled by (i) the dissolution process and (ii) the nonlinearity arising from the variation in density. Assuming a constant temperature and negligible fluid compressibility, the fluid density is assumed to vary linearly with the concentration as follows:

$$ \rho = \rho_0 + (\rho_1 - \rho_0)c \quad (6) $$

where $\rho_1$ is the density of the contaminant fluid.

Because of dissolution, the aperture varies with time as follows:

$$ \rho_r \delta \frac{\partial b}{\partial t} - R(c) = 0 \quad (7) $$

where $\rho_r$ is the rock density $[M.L^{-3}]$ and $\delta$ is the stoichiometric coefficient $[-]$ corresponding to the mass of solute for a unit mass of dissolved rock.

A first-order approximation is used for the reaction term [38,39]:

$$ R(c) = k_c (c_s - c) \quad (8) $$

where $k_c$ is the dissolution coefficient $[L.T^{-1}]$ and $c_s$ $[M.L^{-3}]$ is the concentration at saturation.

Two configurations are examined in this work. In the first configuration, the dissolution coefficient remains constant, while in the second, the rate $k_c$ increases with the fluid velocity, as observed with artificial water circulation systems within hot dry rocks [40].

The system of coupled Equations (1)–(8) governing flow, transport and dissolution must be solved with the corresponding initial and boundary conditions. The flow boundary conditions can be of the Dirichlet or Neumann type, fixing either the head or the flux for boundary fractures. The initial flow condition corresponds to the initial head in all fracture elements. The boundary conditions for transport are of the Dirichlet type, if the concentration is prescribed, or of the Neumann type if the dispersive flux is prescribed. The initial conditions for transport correspond to the initial concentration in the fracture network.

The mathematical model formed by Equations (1)–(8) cannot be solved analytically because of its high nonlinearity. In the next section, an efficient numerical model is developed for an accurate solution to the mathematical model.

### 3. The Numerical Model
We consider a two-dimensional fractured aquifer where the rock matrix is assumed to be completely impermeable. Thus, flow, mass transport and dissolution occur only within the connected fractures. An efficient numerical model is developed on the fracture network to address the system of coupled Equations (1)–(8) governing fluid flow, mass transfer and dissolution with variable density and fracture aperture evolution. The MOL is employed by first discretizing the spatial derivatives to obtain a system of nonlinear ODEs. The latter is then solved using high-order time integration methods. The FV method, known for its robustness and cost-effectiveness in the discretization of conservation laws, is employed for the spatial discretization. The resulting coupled system of ODEs is tackled using the DASPK solver for time integration.

### 3.1. Spatial Discretization of the Flow

Neglecting the variation in the porosity inside the fracture, as compared to the fracture aperture evolution, integration over the fracture element $k$ of continuity Equation (1) yields

$$ \phi^k \ell^k \underset{t}{\stackrel{\sim}{\gamma}} b^k \qquad q_1^k + \bar{q}_2^k = 0 \qquad (9) $$

where $\phi^k$ and $\ell$ are the porosity and the length, assumed to be constant for fracture branch $k$. The aperture $b^k$ of $k$ is assumed to vary due to the dissolution process. The fluxes $\bar{q}_1^k$ and $\bar{q}_2^k$ at the two extremity nodes of $k$ are considered positive for outflow (see Figure 1).

<figure><img src="image_3.png" /><figcaption>Figure 1. Flow notation: (a) for the fracture $k$ and (b) for the intersection node $i$.</figcaption></figure>

With $h^k$, $z^k$ and $c^k$ respectively denoting the head, elevation and concentration of fracture branch $k$, and $Tz_i^k$ and $Th_i^k$ respectively denoting the elevation and the head at the extremity node $i$ of $k$, the expression of the flux $\bar{q}_i^k$ is obtained from the discretization of the Darcy Equation (2) as follows:

$$ \bar{q}_i^k = \frac{\rho_0 g (b^k)^2}{12\mu} \ell^k (h^k \quad Th_i^k + r_\rho (z^k - Tz_i^k) c^k) \qquad (10) $$

$$ r_\rho = \frac{\rho_1 - \rho_0}{\rho_0} $$

with Plugging Equation (10) into Equation (9) yields
$$ \phi^k \ell \stackrel{\sim}{\underset{\sim}{\int}}^k \frac{\rho_0 g(b^k)^3}{6\mu\ell} \ell \wedge \ell^k - Th_1^k - Th_2^k = 0 \quad (11) $$

The mass conservation of water at each node $i$ of the intersection of $n$ fractures (see Figure 1b corresponding to $n = 4$) is written as

$$ \sum_{l=1}^{n} \bar{q}_{i}^{l} = Q_{i}^{*} \quad (12) $$

where $Q_i^*$ corresponds to the pumping/injection rate at intersection node $i$.

Combining Equations (10) and (12) gives

$$ \sum_{l=1}^{n} \frac{\rho_0 g(b^l)^3}{6\mu\ell^l} \ell h^l Th_i^l + r_\rho(z^l - Tz_i^l)c^l = Q_i^* \quad (13) $$

The continuity of the head at intersection node $i$ ($Th_i = Th_i^l, \forall l = 1, \dots, n$) yields

$$ Th_i = \frac{\sum_{l=1}^{n} \frac{\rho_0 g(b^l)^3}{6\mu\ell^l} \ell h^l \cdot r_\rho(z^l - Tz_i^l)c^l - Q_i^*}{\sum_{l=1}^{n} \frac{\rho_0 g(b^l)^3}{6\mu\ell^l}} \quad (14) $$

Assuming that the first node of fracture $k$ is shared by $n_1$ fractures and the second node is shared by $n_2$ fractures, the substitution of Equation (14) for $Th_1^k$ and $Th_2^k$ in mass conservation Equation (11) yields the final flow equation to solve:

$$ \phi^k \ell^k \stackrel{\sim}{\underset{\sim}{\int}}^k \frac{\rho_0 g(b^k)^3}{6\mu\ell^k} \left\{ \begin{aligned} & \sum_{l=1}^{n_1} \frac{\rho_0 g(b^l)^3}{6\mu\ell^l} \ell h^l \cdot r_\rho(z^l - Tz_i^l)c^l \cdot Q_i^* \\ & - \sum_{l=1}^{n_2} \frac{\rho_0 g(b^l)^3}{6\mu\ell^l} \ell h^l \cdot r_\rho(z^l - Tz_i^l)c^l \cdot Q_2^* \\ & - \sum_{l=1}^{n_3} \frac{\rho_0 g(b^l)^3}{6\mu\ell^l} \ell h^l \cdot r_\rho(z^l - Tz_i^l)c^l \end{aligned} \right\} = 0 \quad (15) $$

In this equation, the unknowns are the heads, the apertures and the concentrations at fracture branch $k$ and at its sharing fracture branches $l$. The flow system is constructed by writing Equation (15) for all the fracture branches forming the fracture network.

### 3.2. Discretization of Transport and Dissolution

Integrating advection-dispersion transport Equation (5) over $k$ yields

$$ \int_0^k \phi b \frac{\partial c}{\partial t} + \int_0^k b \nabla \cdot (qc) - \int_0^k bc \nabla \cdot q - \int_0^k \nabla \cdot (bD \nabla c) = \int_0^k R(c) \quad (16) $$

which can be rewritten as

$$ \phi^k b^k \ell^k \stackrel{\sim}{\underset{\sim}{\int}}^k c^k q^l \stackrel{\sim}{\underset{\sim}{\int}}^k q^l q^k q^k c^k \ell^k R c^k \quad (17) $$

The total flux $q_m^{t,k}$ at the extremities $m = 1, 2$ of fracture branch $k$ is defined by

$$ q_m^{t,k} = q_m^{adv,k} + q_m^{disp,k} \quad (18) $$
where $q_m^{adv,k}$ and $q_m^{disp,k}$ are, respectively, the advective and dispersive flux at extremity $m$ of $k$.

To avoid unphysical oscillations, an upwind scheme is used for the advective flux as follows:

$$q_m^{adv,k} = \bar{q}_m^k \left[ \underbrace{\lambda_m^k c^k \cdot l_1 \lambda_m^k c_m}_{c_m^{k,*}} \right] \quad (19)$$

where $c^k$ is the concentration at the center of fracture branch $k$, $c_m$ is the concentration at its extremity node $m$, and $\lambda_m^k$ is the upwinding parameter, defined as $\lambda_m^k = 1$ if $(\bar{q}_m^k \ge 0)$ or $\lambda_m^k = 0$ otherwise.

Thus, the advective flux at node $m$ is calculated using the upstream concentration $c_m^{k,*} = \lambda_m^k c^k + (1-\lambda_m^k)c_m$, which corresponds to the concentration of fracture $k$ ($c_m^{k,*} = c^k$) in the case of an outflow ($\bar{q}_m^k \ge 0$) and to the concentration of intersection node $m$ ($c_m^{k,*} = c_m$) in the case of an inflow ($\bar{q}_m^k < 0$).

The mass conservation of the contaminant at intersection node $m$, shared among $n$ fractures, is written as

$$\sum_{l=1}^{n} q_{m}^{t,l} = 0 \quad (20)$$

The dispersive flux at node $m$ of $k$ is approximated by

$$q_m^{disp,k} = \frac{2b^k D_k}{\ell^k} (c^k \quad c_m) \quad (21)$$

where $D_k$ is the dispersion at fracture element $k$, expressed using Equation (4) as

$$D_k = \alpha_L^k |q_k| + D_m^k \quad (22)$$

where $\alpha_L^k$, $|q_k| = \frac{|q_2^k - q_1^k|}{2b^k}$ and $D_m^k$ are, respectively, the dispersivity coefficient, the norm of the velocity and the molecular diffusion in fracture branch $k$.

Substituting Equations (19) and (21) into Equation (20) provides the expression of the concentration at node $m$ as follows:

$$c_m = \frac{\sum_{l=1}^{n} \left( \lambda_m^l \bar{q}_m^l + \frac{2b^l D_l}{\ell^k} \right) c^l}{\sum_{l=1}^{n} \left( -\ell^l D_l \quad \bar{q}_m^l (1-\lambda_m^l) \right)} \quad (23)$$

Plugging Equations (19), (21) and (23) into Equation (18), which we substitute into transport Equation (17), yields the final transport equation:
$$ \phi^k b^k \ell^k \stackrel{\sim}{\sim} c^k \quad \begin{array}{c} \left( \phi^k b^k \ell^k \right) \cdot \left( -\lambda_1^k \right) \bar{q}_1^k \left[ \frac{\sum_{l=1}^{n_1} \left( \lambda_1^l \bar{q}_1^l + \frac{2b^l D_l}{\ell^k} \right) c_l^l}{\sum_{l=1}^{n_1} \left( \frac{2b^l}{\ell^l} D_l \frac{\bar{q}^l}{\lambda^l} \right) \lambda^l} \right] \left[ \frac{\sum_{l=1}^{n_2} \left( \lambda_2^l \bar{q}_2^l + \frac{2b^l D_l}{\ell^k} \right) c_l^l}{\sum_{l=1}^{n_2} \left( \frac{2b^l}{\ell^l} D_l \frac{\bar{q}^l}{\lambda^l} \right) \lambda^l} \right] \end{array} \quad (24) \\ + \frac{2b^k D_k}{\ell^k} \left[ \begin{array}{c} \left( \phi^k b^k \ell^k \right) \cdot \left( -\lambda_1^k \right) \bar{q}_1^k \left[ \frac{\sum_{l=1}^{n_1} \left( \lambda_1^l \bar{q}_1^l + \frac{2b^l D_l}{\ell^k} \right) c_l^l}{\sum_{l=1}^{n_1} \left( \frac{2b^l}{\ell^l} D_l \frac{\bar{q}^l}{\lambda^l} \right) \lambda^l} \right] - \frac{\sum_{l=1}^{n_2} \left( \lambda_2^l \bar{q}_2^l + \frac{2b^l D_l}{\ell^k} \right) c_l^l}{\sum_{l=1}^{n_2} \left( \frac{2b^l}{\ell^l} D_l \frac{\bar{q}^l}{\lambda^l} \right) \lambda^l} \end{array} \right] = \ell^k k_c c_s c^k \quad (25) $$

Finally, dissolution Equation (7) is integrated over fracture $k$ and written as follows using Equation (8):

$$ \rho_r \delta \frac{\partial b^k}{\partial t} - k_c (c_s - c^k) = 0 \quad (25) $$

To sum up, the final nonlinear flow-transport-dissolution system with variable density and variable aperture is formed by the coupled Equations (15), (24) and (25). This system has heads $h^k$, concentrations $c^k$ and apertures $b^k$ for all fracture branches $k$ as its unknowns. Hence, the size of the nonlinear flow-transport-dissolution system is triple the number of fractures.

### 3.3. Temporal Discretization

The resulting nonlinear system is integrated with respect to time using high-order methods due to their efficiency for nonlinear systems compared to lowest-order methods [32]. A single system of ODEs is constructed as follows:

$$ G(t, \mathbf{y}, \mathbf{y}') = 0 \quad (26) $$

where the vector $\mathbf{y} = [h^k, c^k, b^k]_{k=1,...,nb\_fractures}$ corresponds to the unknown head, concentration and aperture at all the fracture branches of the network.

For the solution of system (26), we use DASPK, a sophisticated solver for ODE systems and differential-algebraic equation (DAE) systems [31]. DASPK employs the Fixed Leading Coefficient Backward Difference Formulas (FLCBDFs) [41–43] of varying orders up to five. DASPK incorporates an efficient automatic time-stepping scheme. It dynamically adjusts both the order of integration and the time step length to improve efficiency and, at the same time, ensure a small temporal error. The Newton Method is used to linearize system (26), where the Jacobian matrix $J$ is approximated by the finite difference method. To reduce the computational cost, the same Jacobian is utilized for multiple time steps, with updates only occurring when necessary (when convergence fails). The efficiency of the DASPK time solver has been demonstrated for nonlinear coupled processes in porous media [33].

## 4. Numerical Experiments

This section presents numerical experiments conducted using the developed VDF-DFN model for the simulation of two test problems. The first example is a variant of the Henry problem involving saltwater intrusion with dissolution in a fractured coastal aquifer. In this scenario, a constant fracture aperture is assumed, and the results obtained using different fracture networks are compared to the semi-analytical solution established in the literature [44] for saltwater intrusion with dissolution in a coastal aquifer. In the second example, the VDF-DFN model is utilized to simulate coupled fluid flow, mass transfer and dissolution for an injection/extraction well pair problem. Here, we explore
the impact of dissolution on the evolution of the fracture aperture, considering both constant and variable dissolution rates.

## 4.1. The Henry Saltwater Intrusion Problem with Dissolution

The Henry problem [45] serves as a widely recognized illustration of saltwater intrusion within a hypothetical rectangular coastal aquifer 2 m in length and 1 m high (Figure 2). The original problem deals with a confined coastal aquifer where saltwater intrudes from the right boundary due to its higher density. This intrusion persists until equilibrium is attained with the freshwater flowing in from the opposite left vertical boundary, as illustrated in Figure 2. Henry's original semi-analytical solution, developed in 1964, determined steady-state isochores using the Fourier Galerkin (FG) method. The FG method is based on the expansion of the stream function and the concentration with Fourier series. However, the initial FG semi-analytical solution derived by Henry was deemed insufficiently accurate, primarily due to its low number of terms in the Fourier series [46]. Subsequent studies introduced more precise semi-analytical solutions, incorporating a greater number of terms [46,47].

In the following, we consider the semi-analytical solution developed by Younes and Fahs [44], addressing the complex scenario of saltwater intrusion with dissolution and low diffusion. The chosen diffusion coefficient was 10 times smaller than the original coefficient used by Henry. To ensure the stability of the FG solution, a substantial increase in the number of coefficients in the Fourier expansion was necessary, ultimately reaching 6195. Table 1 shows the parameterization and boundary conditions. The computed concentration is relative to that of the sea; hence, it is dimensionless.

<table><caption>**Table 1.** Model parameterizations for the Henry problem with dissolution and low dispersion.</caption><tr><td>Permeability</td><td>k<sub>0</sub> = 1.0204×10<sup>-9</sup> m<sup>2</sup></td></tr><tr><td>Porosity</td><td>φ = 0.35</td></tr><tr><td>Dispersivities</td><td>α<sub>L</sub> = α<sub>T</sub> = 0 m</td></tr><tr><td>Molecular diffusion</td><td>D<sub>m</sub> = 6.6×10<sup>-7</sup> m<sup>2</sup>/s</td></tr><tr><td>Dissolution rate</td><td>k<sub>c</sub> = 10<sup>-6</sup> m/s</td></tr><tr><td>Dimensioless concentration at saturation</td><td>c<sub>s</sub> = 10</td></tr><tr><td>Flow boundary conditions</td><td>- Hydrostatic distribution at right vertical boundary<br/>- Fixed flux at left vertical side: Q<sub>0</sub> = 6.6×10<sup>-5</sup> m<sup>2</sup>/s<br/>- Top and bottom are impermeable boundaries<br/>- Zero concentration c = 0 on left vertical side (ρ<sub>0</sub> = 1000 kg/m<sup>3</sup>).</td></tr><tr><td>Transport boundary conditions</td><td>- Dimensionless concentration c = 1 on right-hand side (ρ<sub>1</sub> = 1025 kg/m<sup>3</sup>).<br/>- Zero normal diffusion along top and bottom</td></tr></table>
In the DFN simulation of this problem, the coastal aquifer is considered fractured with a uniform orthogonal distribution of the fractures, such as in [7]. The network contains (n<sub>f</sub> + 1) horizontal fractures and (2n<sub>f</sub> + 1) vertical fractures (Figure 2). Density-driven flow and dissolution processes exclusively take place within the fracture network. To be equivalent to the unfractured case, the injection of freshwater flux is fixed to Q<sub>0</sub>/(nf + 1) for the intersection nodes that touch the left vertical boundary of the domain. All fractures maintain fixed permeability k<sub>f</sub> = k<sub>0</sub> and constant aperture

$$b = \frac{1}{(nf + 1)}.$$

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Conceptual model of the fractured (n<sub>f</sub> = 5) coastal aquifer for saltwater intrusion with dissolution.</figcaption></figure>

Figure 3a shows the long-time steady-state distribution of the concentration for calculations using a coarse network corresponding to nf = 10. As depicted in Figure 3a, the saltwater intrusion occurs at the right side, moving to the left and ultimately reaching equilibrium. Remarkably, the intrusion is more prominent towards the bottom, predominantly influenced by density effects. Notice that the substantial amount of saltwater in the domain is attributed to both dissolution and intrusion.

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** A solution for the Henry problem on a DFN involving dissolution and low dispersion (a) and comparison with the semi-analytical solution for steady-state conditions (b).</figcaption></figure>

Figure 3b shows a comparison between the FG semi-analytical solution and the VDF-DFN numerical solution obtained using a coarse fracture network corresponding to $nf = 10$ and a dense fracture network corresponding to $nf = 40$. The results of this figure show very good agreement between the semi-analytical solution and the VDF-DFN results obtained on the dense fracture network. This strong agreement serves as a validation for the developed VDF-DFN numerical model for density-driven flow with dissolution.

The developed VDF-DFN can be used with more general fracture networks (not necessarily orthogonal) involving different permeabilities and/or apertures. As an example, Figure 4 shows the concentration distribution on a general fracture network formed by the combination of a first orthogonal network corresponding to $nf = 20$ with a second network formed of a random distribution of fractures. The random fractures are twice as permeable as orthogonal fractures, and their aperture is 5 times smaller than that of regular fractures. The number of fracture elements of the original orthogonal fracture network corresponding to $nf = 20$ is 1660. Adding the random fractures yields a general fracture network formed of 3304 fracture elements. The steady-state results of the Henry saltwater intrusion problem with dissolution on the general fracture network show a greater amount of salt in the domain (Figure 4), reflecting an increase in saltwater intrusion and dissolution induced by the added random fractures.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The concentration distribution (represented as a continuous 2D plot) for the Henry salt-water intrusion problem with dissolution using a general fracture network formed of the uniform orthogonal fracture network corresponding to $nf = 40$, to which we added a random distribution of highly permeable fractures.</figcaption></figure>

4.2. *Fracture Evolution Due to Dissolution in an Injection/Extraction System*
For this test case, we consider a 100m × 100m square fractured horizontal aquifer with a uniform fracture network formed of (n<sub>f</sub> + 1) longitudinal fractures and (n<sub>f</sub> + 1) transverse fractures (Figure 5). The flow field within this aquifer is generated by an extraction and injection well system, with uniform regional flow from the left boundary, where the hydraulic head is imposed to 101 m, to the right boundary, where the head is imposed to 100 m. The separation distance between the extraction and injection wells is 12 m, with the injection well located at (44 m, 50 m) and the pumping well located at (56 m, 50 m). For the transport, a Dirichlet boundary (c = 1) is prescribed at the left side, and a Neumann condition corresponding to zero diffusive flux is prescribed on the right side. Initially, all fractures are saturated with the solute (c = 1). Freshwater (c = 0) intrudes from the injection well. A uniform dense fracture network based on n<sub>f</sub> = 50 is employed.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The boundary conditions and spatial configuration of the injection and pumping wells for the injection/extraction well pair problem.</figcaption></figure>

The permeability of the fractures is ruled by the cubic law, and it changes during the simulation due to the evolution of the aperture induced by dissolution. The simulation was conducted over a time span of T<sub>f</sub> = 10<sup>6</sup> s using the parameters depicted in Table 2. To examine the impact of the dissolution coefficient on the flow behavior and final concentration distribution, simulations were carried out for the following three cases: (a) without dissolution, (b) with a fixed dissolution rate k<sub>c</sub><sup>0</sup> and (c) with a dissolution rate k<sub>c</sub> that increases with the fluid velocity, as suggested in [40]. For this last case, the dissolution coefficient in the fracture is ruled by k<sub>c</sub> = k<sub>c</sub><sup>0</sup>(1 + 10|q<sup>k</sup>|), where |q<sup>k</sup>| = 1/(2b<sup>k</sup>)|q<sub>2</sub><sup>k</sup> - q<sub>1</sub><sup>k</sup>| is the norm of the velocity in fracture branch k. For the studied problem, k<sub>c</sub> varies from k<sub>c</sub><sup>0</sup> (for fractures far from the wells) to 2k<sub>c</sub><sup>0</sup> for fractures with high velocity (near the wells).

<table><caption>**Table 2.** Parameters for the injection/extraction well pair problem.</caption><tr><td>Aperture of Fractures</td><td>b = 1 mm</td></tr></table>
<table><tr><td>Porosity</td><td>$\phi = 0.5$</td></tr><tr><td>Pumping/injection rate</td><td>$Q^* = 10^{-3}$ m²/s</td></tr><tr><td>Dispersivity</td><td>$\alpha_L = 0$ m</td></tr><tr><td>Molecular diffusion</td><td>$D_m = 10^{-8}$ m²/s</td></tr><tr><td>Dissolution rate</td><td>$k_C^0 = 10^{-6}$ m/s</td></tr><tr><td>Dimensionless concentration at saturation</td><td>$c_S = 1$</td></tr><tr><td>Density of rock</td><td>$\rho_r = 2700$ kg/m³</td></tr><tr><td>Stoichiometric coefficient</td><td>$\delta = 0.4$</td></tr></table>

Figure 6 shows the head, concentration and aperture distributions for the three investigated cases. Although these results were obtained on the discrete fracture network, all distributions are represented as continuous 2D plots. In the case without dissolution (case a), where the aperture of all fractures remains unchanged and the corresponding permeability is relatively small according to the cubic law, significant hydraulic effects are observed. This includes a strong depression corresponding to a head of 99.1 m, observed in the pumping well, and a strong overpressure corresponding to a head of 101.9 m in the injection well. As a consequence, a large zone of influence of the injection/extraction well pair is observed (Figure 6(a1)). The concentration distribution shows that the injected clean water occupies a large zone of an almost circular shape around the wells (Figure 6(a2)). Notably, not all the freshwater from the injection well is captured by the pumping well, and a portion of the injected freshwater escapes, reaching the right boundary (Figure 6(a2)).

In case (b), where dissolution occurs at a fixed rate, $k_C^0$, the injected freshwater induces the enlargement of the fractures (Figure 6(b3)). Figure 7 shows the evolution of the aperture of the horizontal fracture located at the center between the injection and extraction wells. This figure shows an almost linear increase in the aperture over time due to dissolution. Close to the injection well, the aperture of the fractures nearly doubles, reaching 1.9 mm. As a result, the permeability of the fractures where dissolution is important experiences a significant increase (due to the cubic law). Hence, the fluid moves more easily from the injection well to the pumping well. The depression observed in the pumping well and the overpressure in the injection well are less important than for case (a). The pressure in the injection well is 100.7 m, while in the pumping well, it reaches 100.2 m (Figure 6(b1)). The zone of influence of the injection/extraction is reduced. The freshwater distribution exhibits an almost circular shape with a smaller diameter around the wells (Figure 6(b2)). It is important to note that the aperture distribution shows more significant enlargement in the vicinity of the injection well compared to the region near the pumping well (Figure 6(b3)). This observed phenomenon aligns with the underlying physics. Indeed, the dissolution is more significant around the injection well, where the concentration is nearly zero, as compared to the area near the extraction well. In the latter case, the concentration of the fluid is non-zero since it has been contaminated by dissolution during its travel from the injection well. In case (c), where the dissolution rate $k_C$ varies with the velocity, the enlargement of the fracture apertures is even more important than in the previous two cases, with values reaching 3.7 mm near the injection well. The aperture of the horizontal fracture located at the center between the injection and extraction wells increases linearly over time, with a more significant slope than in case (b) (see Figure 7). Figure 6(c3) shows that the aperture is significantly enlarged for the fractures between the
injection and pumping wells, where the velocity is significant. This results in an aperture
distribution that takes on an almost tunnel-like shape between the two wells (Figure
6(c3)). The presence of this tunnel with a significantly high aperture (and, consequently,
high permeability) allows a more efficient connection between the injection and pumping
wells. The overall pressure distribution exhibits minimal impacts from the well pair (Fig-
ure 6(c1)), and the zone of influence for the pair of wells is reduced. The concentration
distribution shown in Figure 6(c2) reveals that the freshwater occupies a smaller area with
an ellipsoid shape situated between the two wells.
<table><thead><tr><th></th><th>Head distribution</th><th>Concentration distribution</th><th>Aperture distribution</th></tr></thead><tbody><tr><td>Case (a)</td><td><img src="a1_image.png" alt="Case (a) head distribution"></td><td><img src="a2_image.png" alt="Case (a) concentration distribution"></td><td></td></tr><tr><td>Case (b)</td><td><img src="b1_image.png" alt="Case (b) head distribution"></td><td><img src="b2_image.png" alt="Case (b) concentration distribution"></td><td><img src="b3_image.png" alt="Case (b) aperture distribution"></td></tr></tbody></table>

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

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

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

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

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 6.** The hydraulic head, concentration and aperture distribution for the injection/extraction well pair problem in cases (a1,a2) (without dissolution), (b1–b3) (dissolution with a fixed rate, $k_c^0$) and (c1–c3) (dissolution with a variable rate, $k_c$).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 7.** The aperture evolution of the horizontal fracture located at the center of the domain (between the injection and extraction wells).</figcaption></figure>

**5. Conclusions**

In this work, we considered coupled fluid flow, mass transfer and dissolution in fractured media with variable density and dynamic fracture aperture evolution. We adopted the discrete fracture network (DFN) model, considering that flow, mass transfer and dissolution processes occur only within the connected fractures. The rock matrix was assumed to be completely impermeable. Modeling these three coupled processes in a DFN can present significant computational challenges due to the high nonlinearities stemming from density variations and aperture modification. In this study, we developed an efficient numerical simulator that simultaneously solves flow, transport and dissolution equations with density variations. The spatial discretization uses the robust and cost-effective finite volume method. To prevent unphysical oscillations for advection-dominated cases, we devised an upwind scheme tailored specifically for the fracture network. Upstream concentrations at intersection nodes are computed by considering all connected fracture concentrations with outflows at that node. Integration with respect to time is conducted using sophisticated ODE solvers, which enable the utilization of high-order methods. This allows for dynamic adjustment of the size of the time step and order of integration during the simulation, crucial for optimizing computational efficiency and, at the same time, having good accuracy.

The validity of the developed VDF-DFN model was established by comparing it with the semi-analytical solution of a variant of the Henry problem, which incorporates salt-water intrusion and dissolution in a fractured coastal aquifer. The comparison demonstrated excellent agreement between the semi-analytical solution and the VDF-DFN results obtained on a dense fracture network, thereby validating the accuracy of our approach.

The developed model was then employed to explore the evolution of fracture apertures induced by dissolution in a problem involving both injection and extraction wells. The flow field in this problem is generated by the extraction and injection system, with uniform regional flow. If no dissolution occurs, a large zone of influence of the injection/extraction well pair is observed. In contrast, when dissolution occurs at a fixed rate, the injected freshwater triggers an enlargement of the fractures. The depression observed in the pumping well is less important, and the zone of influence of the injection/extraction wells is reduced. In the final case, where the dissolution rate varies with velocity, the ap-
erture enlargement is even more important compared to the previous two cases. The ap-
erture is significantly enlarged for the fractures between the injection and the pumping
wells, where the velocity is significant. The aperture distribution shows an almost tunnel-
like shape, characterized by high permeability, which connects the two wells. The overall
pressure distribution is weakly affected by the pair of wells, and the concentration distri-
bution reveals a reduced zone with an ellipsoid shape located between the two wells. This
test case points out the benefit of efficient DFN models, which can be suitable and cost-
effective tools for investigating coupled transfer processes in highly fractured aquifers.

**Author Contributions:** All authors contributed to the study analysis. The model was developed by A.Y. and M.F. Conception and numerical simulations were performed by A.Y. The first draft of the manuscript was written by A.Y. and H.M.B. L.G. revised the manuscript and the literature review. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.

**Acknowledgments:** This work was supported by the French National Center of Research (CNRS) and Strasbourg University.

**Conflicts of Interest:** The authors declare no competing interests.

References

1. Jafari Raad, S.M.; Hassanzadeh, H. Onset of Density-Driven Instabilities in Fractured Aquifers. *Phys. Rev. E* **2018**, *97*, 043109.
https://doi.org/10.1103/PhysRevE.97.043109.

2. Kim, M.; Kim, K.; Han, W.S.; Oh, J.; Park, E. Density-Driven Convection in a Fractured Porous Media: Implications for Geological CO2 Storage. *Water Resour. Res.* **2019**, *55*, 5852–5870. https://doi.org/10.1029/2019WR024822.

3. Wang, G.; Ma, X.; Song, X.; Li, G. Modeling Flow and Heat Transfer of Fractured Reservoir: Implications for a Multi-Fracture Enhanced Geothermal System. *J. Clean. Prod.* **2022**, *365*, 132708. https://doi.org/10.1016/j.jclepro.2022.132708.

4. Gao, X.; Li, T.; Zhang, Y.; Kong, X.; Meng, N. A Review of Simulation Models of Heat Extraction for a Geothermal Reservoir in an Enhanced Geothermal System. *Energies* **2022**, *15*, 7148. https://doi.org/10.3390/en15197148.

5. Grillo, A.; Logashenko, D.; Stichel, S.; Wittum, G. Simulation of Density-Driven Flow in Fractured Porous Media. *Adv. Water Resour.* **2010**, *33*, 1494–1507. https://doi.org/10.1016/j.advwatres.2010.08.004.

6. Werner, A.D.; Bakker, M.; Post, V.E.A.; Vandenbohede, A.; Lu, C.; Ataie-Ashtiani, B.; Simmons, C.T.; Barry, D.A. Seawater Intrusion Processes, Investigation and Management: Recent Advances and Future Challenges. *Adv. Water Resour.* **2013**, *51*, 3–26. https://doi.org/10.1016/j.advwatres.2012.03.004.

7. Sebben, M.L.; Werner, A.D.; Graf, T. Seawater Intrusion in Fractured Coastal Aquifers: A Preliminary Numerical Investigation Using a Fractured Henry Problem. *Adv. Water Resour.* **2015**, *85*, 93–108. https://doi.org/10.1016/j.advwatres.2015.09.013.

8. Hosseini, N.; Bajalan, Z.; Khoei, A.R. Numerical Modeling of Density-Driven Solute Transport in Fractured Porous Media with the Extended Finite Element Method. *Adv. Water Resour.* **2020**, *136*, 103453. https://doi.org/10.1016/j.advwatres.2019.103453.

9. Metsias, G.; Hamill, G.A.; Campbell, D.; Straney, R.; Benner, E.M.; Águila, J.F.; McDonnell, M.C.; Ahmed, A.A.; Flynn, R. Laboratory and Numerical Investigation of Saline Intrusion in Fractured Coastal Aquifers. *Adv. Water Resour.* **2021**, *149*, 103866.
https://doi.org/10.1016/j.advwatres.2021.103866.

10. Kim, J.; Sonnenthal, E.; Rutqvist, J. A Sequential Implicit Algorithm of Chemo-Thermo-Poro-Mechanics for Fractured Geothermal Reservoirs. *Comput. Geosci.* **2015**, *76*, 59–71. https://doi.org/10.1016/j.cageo.2014.11.009.

11. Medekova, A.; Jones, G.D. Characterization and Modeling Challenges Associated with Fracture and Karst (Non-Matrix) in the Margin Area of a Carbonate Reservoir. In Proceedings of the SPE Annual Caspian Technical Conference and Exhibition, Astana, Kazakhstan, 12–14 November 2014; SPE: Astana, Kazakhstan, 2014; p. SPE-172275-MS.

12. Tran, M.; Jha, B. Effect of Poroelastic Coupling and Fracture Dynamics on Solute Transport and Geomechanical Stability. *Water Resour. Res.* **2021**, *57*, e2021WR029584. https://doi.org/10.1029/2021WR029584.
13. Natarajan, N.; Suresh Kumar, G. Solute Transport in a Coupled Fracture-Matrix System with Sinusoidal Fracture Geometry. *Int. J. Eng. Sci. Technol.* **2010**, *2*, 1886–1992.

14. Zechner, E.; Konz, M.; Younes, A.; Huggenberger, P. Effects of Tectonic Structures, Salt Solution Mining, and Density-Driven Groundwater Hydraulics on Evaporite Dissolution (Switzerland). *Hydrogeol. J.* **2011**, *19*, 1323–1334. https://doi.org/10.1007/s10040-011-0759-5.

15. Tenthorey, E.; Fitzgerald, J. Feedbacks between Deformation, Hydrothermal Reaction and Permeability Evolution in the Crust: Experimental Insights. *Earth Planet. Sci. Lett.* **2006**, *247*, 117–129. https://doi.org/10.1016/j.epsl.2006.05.005.

16. Deng, H.; Spycher, N. Modeling Reactive Transport Processes in Fractures. *Rev. Mineral. Geochem.* **2019**, *85*, 49–74. https://doi.org/10.2138/rmg.2019.85.3.

17. Ameli, P.; Elkhoury, J.E.; Morris, J.P.; Detwiler, R.L. Fracture Permeability Alteration Due to Chemical and Mechanical Processes: A Coupled High-Resolution Model. *Rock. Mech. Rock. Eng.* **2014**, *47*, 1563–1573. https://doi.org/10.1007/s00603-014-0575-z.

18. Aliouache, M.; Wang, X.; Jourde, H.; Huang, Z.; Yao, J. Incipient Karst Formation in Carbonate Rocks: Influence of Fracture Network Topology. *J. Hydrol.* **2019**, *575*, 824–837. https://doi.org/10.1016/j.jhydrol.2019.05.082.

19. Lopes, J.A.G.; Medeiros, W.E.; La Bruna, V.; De Lima, A.; Bezerra, F.H.R.; Schiozer, D.J. Advancements towards DFKN Modeling: Incorporating Fracture Enlargement Resulting from Karstic Dissolution in Discrete Fracture Networks. *J. Pet. Sci. Eng.* **2022**, *209*, 109944. https://doi.org/10.1016/j.petrol.2021.109944.

20. Berre, I.; Doster, F.; Keilegavlen, E. Flow in Fractured Porous Media: A Review of Conceptual Models and Discretization Approaches. *Transp. Porous Med.* **2019**, *130*, 215–236. https://doi.org/10.1007/s11242-018-1171-6.

21. Liu, R.; Li, B.; Jiang, Y.; Huang, N. Review: Mathematical Expressions for Estimating Equivalent Permeability of Rock Fracture Networks. *Hydrogeol. J.* **2016**, *24*, 1623–1649. https://doi.org/10.1007/s10040-016-1441-8.

22. Kordilla, J.; Sauter, M.; Reimann, T.; Geyer, T. Simulation of Saturated and Unsaturated Flow in Karst Systems at Catchment Scale Using a Double Continuum Approach. *Hydrol. Earth Syst. Sci.* **2012**, *16*, 3909–3923. https://doi.org/10.5194/hess-16-3909-2012.

23. Flemisch, B.; Berre, I.; Boon, W.; Fumagalli, A.; Schwenck, N.; Scotti, A.; Stefansson, I.; Tatomir, A. Benchmarks for Single-Phase Flow in Fractured Porous Media. *Adv. Water Resour.* **2018**, *111*, 239–258. https://doi.org/10.1016/j.advwatres.2017.10.036.

24. Hoteit, H.; Firoozabadi, A. Modeling of Multicomponent Diffusions and Natural Convection in Unfractured and Fractured Media by Discontinuous Galerkin and Mixed Methods. *Numer. Meth Eng.* **2018**, *114*, 535–556. https://doi.org/10.1002/nme.5753.

25. Shakiba, M.; Cavalcante Filho, J.S.D.A.; Sepehrnoori, K. Using Embedded Discrete Fracture Model (EDFM) in Numerical Simulation of Complex Hydraulic Fracture Networks Calibrated by Microseismic Monitoring Data. *J. Nat. Gas. Sci. Eng.* **2018**, *55*, 495–507. https://doi.org/10.1016/j.jngse.2018.04.019.

26. Ţene, M.; Bosma, S.B.M.; Al Kobaisi, M.S.; Hajibeygi, H. Projection-Based Embedded Discrete Fracture Model (pEDFM). *Adv. Water Resour.* **2017**, *105*, 205–216. https://doi.org/10.1016/j.advwatres.2017.05.009.

27. Adler, P.M.; Thovert, J.-F. *Fractures and Fracture Networks*; Theory and Applications of Transport in Porous Media; Springer: Dordrecht, The Netherlands, 1999; Volume 15, ISBN 978-90-481-5192-9.

28. Aldrich, G.; Hyman, J.D.; Karra, S.; Gable, C.W.; Makedonska, N.; Viswanathan, H.; Woodring, J.; Hamann, B. Analysis and Visualization of Discrete Fracture Networks Using a Flow Topology Graph. *IEEE Trans. Visual. Comput. Graph.* **2017**, *23*, 1896–1909. https://doi.org/10.1109/TVCG.2016.2582174.

29. Tabrizinejadas, S.; Younes, A.; Hoteit, H.; Carrayrou, J.; Fahs, M. Coupled Fluid Flow, Solute Transport and Dissolution Processes in Discrete Fracture Networks: An Advanced Discontinuous Galerkin Model. *Adv. Water Resour.* **2023**, *180*, 104540. https://doi.org/10.1016/j.advwatres.2023.104540.

30. Eymard, R.; Gallouët, T.; Herbin, R. Finite Volume Methods. In *Handbook of Numerical Analysis*; Elsevier: Amsterdam, The Netherlands, 2000; Volume 7, pp. 713–1018, ISBN 978-0-444-50350-3.

31. Brown, P.N.; Hindmarsh, A.C.; Petzold, L.R. Using Krylov Methods in the Solution of Large-Scale Differential-Algebraic Systems. *SIAM J. Sci. Comput.* **1994**, *15*, 1467–1488. https://doi.org/10.1137/0915088.

32. Younes, A.; Fahs, M.; Ahmed, S. Solving Density Driven Flow Problems with Efficient Spatial Discretizations and Higher-Order Time Integration Methods. *Adv. Water Resour.* **2009**, *32*, 340–352. https://doi.org/10.1016/j.advwatres.2008.11.003.

33. Younes, A.; Koohbor, B.; Belfort, B.; Ackerer, P.; Doummar, J.; Fahs, M. Modeling Variable-Density Flow in Saturated-Unsaturated Porous Media: An Advanced Numerical Model. *Adv. Water Resour.* **2022**, *159*, 104077. https://doi.org/10.1016/j.advwatres.2021.104077.

34. Younes, A.; Hoteit, H.; Helmig, R.; Fahs, M. A Robust Fully Mixed Finite Element Model for Flow and Transport in Unsaturated Fractured Porous Media. *Adv. Water Resour.* **2022**, *166*, 104259. https://doi.org/10.1016/j.advwatres.2022.104259.
35. Li, Q.; Wu, J.; Li, Q.; Wang, F.; Cheng, Y. Sediment Instability Caused by Gas Production from Hydrate-Bearing Sediment in Northern South China Sea by Horizontal Wellbore: Sensitivity Analysis. *Nat. Resour. Res.* **2025**, *34*, 1667–1699. https://doi.org/10.1007/s11053-025-10478-x.

36. Li, Q.; Liu, J.; Wang, S.; Guo, Y.; Han, X.; Li, Q.; Cheng, Y.; Dong, Z.; Li, X.; Zhang, X. Numerical Insights into Factors Affecting Collapse Behavior of Horizontal Wellbore in Clayey Silt Hydrate-Bearing Sediments and the Accompanying Control Strategy. *Ocean. Eng.* **2024**, *297*, 117029. https://doi.org/10.1016/j.oceaneng.2024.117029.

37. Witherspoon, P.A.; Wang, J.S.Y.; Iwai, K.; Gale, J.E. Validity of Cubic Law for Fluid Flow in a Deformable Rock Fracture. *Water Resour. Res.* **1980**, *16*, 1016–1024. https://doi.org/10.1029/WR016i006p01016.

38. Detwiler, R.L.; Rajaram, H. Predicting Dissolution Patterns in Variable Aperture Fractures: Evaluation of an Enhanced Depth-averaged Computational Model. *Water Resour. Res.* **2007**, *43*, 2006WR005147. https://doi.org/10.1029/2006WR005147.

39. Hanna, R.B.; Rajaram, H. Influence of Aperture Variability on Dissolutional Growth of Fissures in Karst Formations. *Water Resour. Res.* **1998**, *34*, 2843–2853. https://doi.org/10.1029/98WR01528.

40. Suto, Y.; Tanifuji, K.; Watanabe, K.; Hashida, T.; Shoji, T.; Takahashi, H. Dependence of the Rate of Granite Dissolution on Temperature and Fluid Velocity under Simulated Geothermal Reservoir Environments. *J. Geotherm. Res. Soc. Jpn.* **2002**, *24*, 47–56.

41. Tocci, M.D.; Kelley, C.T.; Miller, C.T. Accurate and Economical Solution of the Pressure-Head Form of Richards' Equation by the Method of Lines. *Adv. Water Resour.* **1997**, *20*, 1–14. https://doi.org/10.1016/S0309-1708(96)00008-5.

42. Kavetski, D.; Binning, P.; Sloan, S.W. Adaptive Backward Euler Time Stepping with Truncation Error Control for Numerical Modelling of Unsaturated Fluid Flow. *Numer. Meth Eng.* **2002**, *53*, 1301–1322. https://doi.org/10.1002/nme.329.

43. Farthing, M.W.; Kees, C.E.; Miller, C.T. Mixed Finite Element Methods and Higher Order Temporal Approximations for Variably Saturated Groundwater Flow. *Adv. Water Resour.* **2003**, *26*, 373–394. https://doi.org/10.1016/S0309-1708(02)00187-2.

44. Younes, A.; Fahs, M. A Semi-Analytical Solution for the Reactive Henry Saltwater Intrusion Problem. *Water Air Soil Pollut.* **2013**, *224*, 1779. https://doi.org/10.1007/s11270-013-1779-7.

45. Henry, H.R. Interfaces between Salt Water and Fresh Water in Coastal Aquifers. In *Sea Water in Coastal Aquifers*; Supply Paper 1613-C; US Geological Survey Water: Reston, VA, USA, 1964; C35-70.

46. Ségol, G. *Classic Groundwater Simulations: Proving and Improving Numerical Models*; PTR Prentice Hall: Englewood Cliffs, NJ, USA, 1994; ISBN 978-0-13-137993-0.

47. Younes, A.; Fahs, M. A Semi-Analytical Solution for Saltwater Intrusion with a Very Narrow Transition Zone. *Hydrogeol. J.* **2014**, *22*, 501–506. https://doi.org/10.1007/s10040-014-1102-8.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.