# Adaptive mesh refinement for flow and transport problem in heterogeneous porous media

Hanen Amor, Fayssal Benkhaldoun

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

## To cite this version:

Hanen Amor, Fayssal Benkhaldoun. Adaptive mesh refinement for flow and transport problem in heterogeneous porous media. Mathematics and Computers in Simulation, 2026, 239, pp.766 - 789. 10.1016/j.matcom.2025.06.020. hal-05299145

HAL Id: hal-05299145

https://asnr.hal.science/hal-05299145v1

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
<figure><img src="image_3.png" /></figure>

ELSEVIER

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

Original articles

# Adaptive mesh refinement for flow and transport problem in heterogeneous porous media

Hanen Amor<sup>a ID,*</sup>, Fayssal Benkhaldoun<sup>b ID</sup>

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

Check for
updates

## ARTICLE INFO

**Keywords:**

Adaptive mesh refinement
FVFE method
Porous media
Radioactive waste

## ABSTRACT

This paper presents an adaptive mesh refinement strategy for modeling the flow and transport of contaminants in porous media. These phenomena are modeled by a system of equations that combine Darcy's law and the diffusion-convection equation. A combined finite volume and finite element discretization is used to discretize the coupled system. A guaranteed *a posteriori* error estimator, derived from an H(div)-conforming flux reconstruction, is computed and followed by an appropriate method of local refinement and coarsening of mesh elements in two and three dimension to minimize the difference between a discretized solution and a continuous solution in the context of industrial applications. The method is implemented within the MELODIE software developed by the ASNR in order to assess the safety of a deep geological disposal for radioactive waste, and is illustrated through a realistic application of flow and transport in heterogeneous porous media.

## 1. Introduction

Numerical modeling of flow and transport in porous media is essential for addressing environmental engineering applications. The long-term safety of high-level and intermediate long-lived radioactive waste disposal is an important issue in all countries with a significant nuclear program. One solution considered for radioactive waste management in the nuclear industry is disposing it in deep geological underground clay formations chosen for their ability to delay and attenuate potential releases of radionuclides into the biosphere. The disposal is investigated by several countries including France and it typically relies on a multibarrier system to isolate the waste from the biosphere. The multibarrier system typically includes the natural geological barrier provided by the repository host rock and its surroundings, as well as engineered barriers materials (such as backfill, bentonite, concrete, etc.).

The ASNR<fnref n="1" /> has been carrying out independent research and studies for many years in order to simulate the migration of radionuclides from radioactive waste to the biosphere using the MELODIE software<fnref n="2" /> (see e.g. [1,2]). One of the main challenges of such computational modeling is dealing with various geometric scales and the heterogeneity of hydraulic and transport characteristics of the materials in repository components. When assessing confinement capabilities of disposal facilities and of the geological formations, the domain investigated by numerical simulations exhibits significant variability in properties and geometric scales (ranging from centimeters to kilometers), presenting significant computational challenges. Solving these problems involving such phenomena requires high resolution in regions where the numerical solution changes drastically. The numerical method

<footnote><sup>a</sup> ASNR, PSE-ENV/SPDR/UEMIS, France</footnote>

<footnote><sup>b</sup> LAGA, Université Sorbonne Paris Nord, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: hanen.amor@asnr.fr (H. Amor).</footnote>

<footnote n="1">1 Autorité de Sûreté Nucléaire et de Radioprotection - French Authority for Nuclear Safety and Radiation Protection</footnote>

<footnote n="2">2 Modèle d'Evaluation à LOng terme des Déchets Irradiants Enterrés - model for long-term assessment of radioactive waste repositories</footnote>
and mesh should accurately represent the modeled system with the required level of precision to reduce numerical errors while maintaining reasonable computational costs and durations. Therefore, the development of dynamic gridding capability has long been of interest in the community of mathematical modeling of porous media flow and transport. The method must capture regions with higher numerical errors during flow and transport due to higher permeability contrasts, high hydraulic gradients, etc. One solution is to increase the polynomial order of the approximation within elements, but this solution is not easily coded and is costly for practical problems. Another solution is to locally reduce the mesh size. Adaptive mesh refinement and a posteriori error estimation are indispensable for accurately simulating flow and transport in porous media, particularly in high-stakes applications like radioactive waste disposal. Key contributions have advanced the field, addressing multiscale geometries, heterogeneous properties, and computational efficiency. For example, Amaziane et al. (see [3]) highlight finite volume methods tailored for radionuclide transport in highly heterogeneous systems. Ern and Vohralík (see [4,5]) emphasize robust error estimations for heat and diffusion problems, paving the way for adaptive methods that ensure precision while controlling computational costs. Hilhorst et al. [6] explore hybrid finite element/finite volume techniques, demonstrating their effectiveness on non-matching grids. Vohralík [7] extends error estimation frameworks to polytopal meshes, allowing for efficient handling of irregular geometries often encountered in real-world geological systems. Recent works also address thermal and multiphase compositional flow modeling, showing how adaptive mesh refinement can optimize simulations in cases with steep gradients and dynamic interfaces (see [8]). To provide a broader perspective, researchers such as Di Pietro and Vohralík [9] review advances in discretization methods and error analysis for geosciences, presenting a comprehensive outlook on challenges and methodologies. Future directions include integrating machine learning into adaptive algorithms and enhancing parallel computing frameworks to tackle increasingly complex simulations.

The purpose of this paper is to study the performance of specific a posteriori error estimators and to develop a strategy of adaptive mesh refinement for mixed formulations of flow and transport problems in porous media (see e.g. [10,11]). The method will be implemented in the MELODIE software, and numerical simulations will be presented for both academic tests and industrial applications. The chosen estimators are expected to facilitate the management of the local mesh refinement process to minimize numerical errors while maintaining an optimal computational cost. The selected estimators for this study are those developed by Ern, Vohralík et al. (see e.g. [4–6,12]). These estimators were chosen because they satisfy key properties of optimal error estimators, such as guaranteed upper bounds, robustness, and local efficiency. Moreover, their formulation is well-suited for mixed discretizations and can be extended to three-dimensional and nonlinear problems, which aligns with the complexity of the flow and transport problems studied here. In the adaptive mesh refinement process, these estimators are used to identify regions of the domain where the numerical solution requires higher resolution. This process minimizes numerical errors while balancing computational cost.

The paper is organized as follows. The next section provides a brief description of the mathematical equations governing the physics of saturated flow and transport in porous media, along with their discretization using a vertex-centered finite volume method combined with a semi-implicit Euler approach. Following that, the subsequent section outlines the construction of error estimators based on H(div)-conforming flux reconstruction and the refinement process implemented in MELODIE. Results from numerical simulations are presented in the last section. Two tests are investigated for the flow and the transport problems in two dimensions to validate the efficiency and accuracy of the method. The first test focuses on a heterogeneous medium with two different permeability contrasts and a point source, while the second test deals with an example of radionuclide migration released by waste packages in a deep geological repository. This scenario involves a complex heterogeneous layered and faulted porous medium with high permeability contrasts and several source terms. Finally, conclusions and perspectives of this work are provided.

## 2. Mathematical models and numerical methods

The MELODIE software is a numerical tool developed by ASNR for modeling water flow and the transport of radionuclides in saturated and unsaturated porous media within the context of a radioactive waste disposal facility. Indeed, each component of the repository, including the host rock, is represented as an equivalent homogeneous porous medium (see also [13,14]). The physical properties of water, such as density and viscosity, are assumed to remain unaffected by the presence of solutes transported in the water. This approach aligns with studies like Knabner et al. (see [15]), which highlight that changes in water viscosity due to solute interactions, such as with surfactants, are typically negligible in low-concentration scenarios. This assumption is justified by the low concentration of solutes considered in this study, which minimizes their impact on water properties. Additionally, this simplification is commonly used in transport models to balance computational efficiency with the accuracy required for the studied phenomena. The mathematical equations and their discretization for saturated porous media are summarized below.

### 2.1. Conceptual model

To simplify the presentation, a two-dimensional example was considered where the effects of gravity are negligible. The single-phase flow of an incompressible fluid with a dissolved solute in a horizontal porous saturated reservoir $\Omega \subset \mathbb{R}^2$ is governed by Darcy's law and mass conservation (see, for example, [11,16]):

$$ - \vec{\nabla} \cdot (\mathbf{K} \vec{\nabla} P) - q = 0, \quad (1) $$

Here, $\vec{U} = -(\mathbf{K} \vec{\nabla} P)$ represents the Darcy velocity [$m \cdot y^{-1}$], and $P$ is the hydraulic head [$m$], given by $P = \frac{p}{g\rho} + z$, where $p$ is the water hydraulic pressure [$Pa$], $g$ is the gravitational acceleration [$m \cdot y^{-2}$], $\rho$ is the water density [$kg \cdot m^{-3}$], and $z$ is the topographic level [$m$]. $\mathbf{K}$ denotes the permeability tensor of the medium [$m^2$] and $q$ is the flow source term [$y^{-1}$].
The transport of radionuclides in saturated porous media is described by solving the following partial differential equation over a time interval $(0, T_f)$:

$$ \omega R \frac{\partial C}{\partial t} - \vec{V} \cdot (\mathbf{D}\vec{V}C - \vec{U}C) + \lambda \omega RC = S, \quad (2) $$

In this equation, C represents the radionuclide concentration $[Bq \cdot m^{-3}]$, $\omega$ is the porosity, $R^3$ is the retardation factor ($R = K_d \rho_s^{1-\omega} \omega^{\omega}$, where $K_d$ is the distribution coefficient and $\rho_s$ is the dry density of the solid phase), $\lambda$ is the decay constant $[y^{-1}]$, S is the radionuclide source term $[Bq \cdot m^{-3} \cdot y^{-1}]$, and D is the diffusion-dispersion tensor given by:

$$ \mathbf{D} = (\alpha_T |\vec{U}| + \omega d)I + (\alpha_L - \alpha_T) \frac{\vec{U} \otimes \vec{U}}{|\vec{U}|} \quad (3) $$

Here, $\alpha_L$ and $\alpha_T$ are the longitudinal and transverse dispersivities $[m]$, d is the molecular diffusion constant $[m^2 \cdot y^{-1}]$, and $\otimes$ denotes the tensor product. The D expression used in our model is a classical representation of the diffusion-dispersion tensor in porous media. It relies on Fick's law to describe molecular diffusion, which corresponds to the random movement of particles under the effect of concentration gradients and becomes particularly important at low flow velocities. In addition, the model also includes mechanical dispersion, a phenomenon caused by local variations in pore velocity, leading to particle spreading in time and space. Although mechanical dispersion is not directly derived from Fick's law, it is often integrated into the same mathematical framework to provide a complete description of solute transport. This combined approach is widely documented in the literature on transport in porous media, notably in fundamental works on hydrodynamic dispersion equations (e.g., [17]). The system (1)–(2) is completed by appropriate boundary and initial conditions. Both the K and D tensors can be anisotropic, and all other material-dependent parameters can be heterogeneous and radionuclide-dependent.

## 2.2. FVFE method

The finite volume (FV) method is a widely used approach for simulating flow in porous media, primarily due to its local conservation properties and its ability to handle complex geometries. In this framework, the domain is discretized into control volumes, and the governing equations are integrated over each volume. Fluxes are then approximated at the interfaces between volumes. For the flow and transport models considered here, a two-point flux approximation (TPFA) can be employed for the Darcy flow equation. While this approach simplifies flux computation, it may suffer from inaccuracies on non-orthogonal grids. As a result, grid quality plays a critical role in the accuracy of FV-based simulations. For the transport equation, an upwind scheme is typically used to handle advection-dominated regimes, although it can introduce numerical diffusion.

The finite volume–finite element (FVFE) method, while advantageous due to its compatibility with the MELODIE framework and its ability to preserve local mass conservation, also has known limitations. In highly heterogeneous porous media, it may produce non-physical artifacts such as negative concentrations and spurious solute leakage into low-permeability zones. These numerical effects arise from the discretization strategy, which combines explicit upwinding for the convective term with a linear finite element treatment for diffusion. This combination can become overly diffusive or unstable in regions with sharp concentration gradients, unless carefully controlled. Nevertheless, the FVFE method remains appropriate in this context due to its intrinsic modularity, dual-mesh structure, and full integration within MELODIE's existing numerical infrastructure—particularly for multiphysics coupling. Alternative discretization schemes, such as mimetic finite differences or TPFA, were considered during the design phase. However, their integration was found to be less compatible with the modular operator structure and data layout of MELODIE. As a result, our focus has been placed on mitigating numerical artifacts through error-controlled mesh adaptation rather than switching to different schemes. In particular, local negative concentrations are mitigated through careful time-step control and mesh refinement near sharp fronts. The theoretical and practical analyses presented in [18] provide further support for the method's convergence and robustness under appropriate assumptions.

In MELODIE, the spatial system (1)–(2) is discretized using the vertex-centered finite volume method, also known as the control volume finite element method. A Godunov scheme is used for the convection term and a $P_1$-finite element scheme is used for the discretization of the diffusion-dispersion term (see [1,19]). A semi-implicit finite volume scheme is defined, with an explicit approximation of the convection term and an implicit approximation of the diffusion-dispersion term, requiring the construction of a secondary mesh. Indeed, this approach relies on two spatial meshes: a primary mesh, which is a conforming finite element mesh, and a secondary mesh consisting of control volumes centered at the vertices of the primary mesh. Control volumes are constructed around mesh nodes by connecting the midpoints of the edges of a triangle with the triangle's barycenter (see Fig. 1).

The degrees of freedom are located at the vertices of the primary mesh. All model data, including permeability, porosity, sources, and dispersivity, are constant element-wise on the primary mesh. Consequently, the material interfaces align with the edges of the primary mesh. The unknowns of these equation are represented as a $P_1$ finite element functions over the primary mesh. This discretization method is also known as the FVFE (see [1,19,20]). The FVFE method was chosen for its ability to ensure local conservation of flows, a crucial property for transport problems in porous media. It also provides flexibility in handling complex geometries and unstructured meshes, which are often required in geological disposal studies. Additionally, the method ensures

<footnote n="3">³ Retardation factor is a parameter used in modeling to describe the slowdown in the movement of solutes in a porous medium, relative to the flow of water. This slowdown is caused by interactions between solutes and the solid matrix of the medium, such as adsorption (particles bind to solid surfaces) or ionic exchange (chemical exchanges between solutes and solids).</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 1. The primary mesh and the dual mesh.</figcaption></figure>

non-negativity of concentrations under the CFL condition, an essential feature for preserving the physics of the problem. Compared to classical finite volumes or standard finite elements, the FVFE method combines these advantages while offering greater accuracy in complex scenarios.

Considering the partition of the interval $[0, T_f]$ into sub-intervals $[t_{n-1}, t_n]$ such that $0 = t_0 < \dots < t_l < \dots < t_N = T_f$ in the numerical scheme, for each $n$, the FVFE method requires constructing a conforming triangular mesh of the domain $\Omega$ in two dimensions. For each vertex $i \in \mathcal{V}$ of the primary mesh, a control volume $D$ is associated by connecting the midpoints of the edges of adjacent elements with the barycenters of these elements. The mesh obtained by uniting all control volumes $D$, $i \in \mathcal{V}$, forms the dual mesh $D_h$, on which the numerical scheme is entirely built. The FVFE scheme can then be expressed for any $D \in D_h$ as:

$$ - \sum_{\gamma \subset \partial D} \int_{\gamma} \mathbb{K} \nabla P_h \cdot \vec{n}_{\gamma} \, ds = \int_{D} q \, dx, \qquad (4) $$

where $P_h$ denotes the numerical solution of the flow Eq. (1), $\vec{n}_{\gamma}$ is the unit outward normal vector on $\gamma$, $\partial D$ is the set of edges of the dual cell $D$ and $\bar{U}_h = -\mathbb{K} \nabla P_h$ is an approximation of the Darcy velocity.

Denoting $C_h^n$ the numerical solution of the transport Eq. (2) obtained for $t = t^n$, the FVFE scheme can be written $\forall D \in D_h$ with a semi-implicit time stepping, the time discretization is implicit for the diffusion term and it is explicit for the convective term as:

$$ \int_D \omega R \frac{C_h^n - C_h^{n-1}}{\tau_n} \, dx - \sum_{\gamma \subset \partial D} \left\{ \int_{\gamma} \mathbb{D} \nabla C_h^n \cdot \vec{n}_{\gamma} \, ds - \int_{\gamma} \bar{U}_h C_h^{n-1} \cdot \vec{n}_{\gamma} \, ds \right\} + \int_D \lambda \omega R C_h^n \, dx = \int_D S \, dx, \qquad (5) $$

where $\tau_n = t_n - t_{n-1}$ is the current time step. The resulting scheme is designed to be locally conservative and to preserve the non-negativity of radionuclide concentration under a stability CFL condition, although deviations may arise in practice due to factors such as mesh quality (see [10]).

### 3. A posteriori error estimators and adaptive algorithm

In this section, an adaptive numerical technique using the finite volume approximation described in the previous section are proposed. The a posteriori error estimators applied in this study are based on the flux equilibration technique, which is well-established for mixed finite volume and finite element discretizations (see [4–7,21–24]). These estimators have been shown to satisfy key theoretical properties, including guaranteed upper bounds, robustness, and local efficiency, making them suitable for guiding adaptive mesh refinement. The novelty of this work lies in their application to industrial-scale problems, involving complex geometries and highly heterogeneous media. While the theoretical derivation of these estimators is not revisited here, their implementation were adapted to account for the challenges posed by the industrial context. The accuracy of these estimators as error indicators in this specific framework is validated through numerical experiments. In the following sub-sections, the adaptive mesh refinement process developed in MELODIE is presented. More specifically, the criterion of the estimator implemented in the software is given, then the refinement and adaptation management strategy are also presented.

#### 3.1. A posteriori error estimators

Before introducing the a posteriori estimators, some definitions and notations are provided. $S_h$ is a (fine) simplicial mesh obtained from the associated dual mesh $D_h$. An illustration of $S_h$ in the two-dimensional case is presented in Fig. 2(a).
Following the flow equation and the FVFE scheme, the approximate flux $-K\nabla P_h$ verifies local conservation on internal dual cells $D$, even though it is not locally mass conservative on the primary mesh $\mathcal{T}_h$. To build an estimator, the main idea is to introduce a suitable equilibrated flux $\mathbf{t}_h \in \mathbf{H}(\text{div}, \Omega)$ corresponds to the flux of flow within the Darcy equation. The space $H(\text{div}, \Omega)$ denotes the space generated by vector-valued functions whose weak divergence has an integrable square:

$$H(\text{div}, \Omega) = \{v \in [L^2(\Omega)]^d : \nabla \cdot v \in L^2(\Omega)\} \text{ d is the dimension.}$$

Using the divergence theorem, the local conservation property is written in the following form:

$$\int_D \nabla \cdot \mathbf{t}_h \, d\mathbf{x} = \sum_{\gamma \subset \partial D} \int_{\gamma} \mathbf{t}_h \cdot \tilde{\mathbf{n}}_\gamma \, d\mathbf{s} = \int_D f \, d\mathbf{x}, \quad \forall D \in D_h^{\text{int}}. \qquad (6)$$

Defining the dual norm by the following formula see [5] :

$$|||u|||_* = \sup_{\varphi \in H_0^1(\Omega) \setminus \{0\}} \frac{\int_{\Omega} K \nabla u \cdot \nabla \varphi \, d\mathbf{x}}{||\varphi||}$$

Then, an upper bound, combining the residual and the flux estimators, of the error measured with this dual norm $|||\cdot|||_*$ is given by

$$\|P - P_h\|_* \le \left\{ \sum_{D \in D_h} (\eta_{R,D} + \eta_{DF,D})^2 \right\}^{1/2} \qquad (7)$$

where the diffusive flux estimator and the residual estimator are given by

$$\eta_{DF,D} := \|K^{\frac{1}{2}} \nabla P_h + K^{-\frac{1}{2}} \mathbf{t}_h\|_D \text{ and } \eta_{R,D} := m_D \|f - \nabla \cdot \mathbf{t}_h\|_D$$

where $m_D$ is a computable constant depending on the diameter of the dual cell $D$ (see, e.g., [5,6]). The norm $\|\cdot\|_D$ refers to a $L^2$ norm defined over the dual cell $D$. For a scalar function $v$, it is typically expressed as:

$$\|v\|_D = (\int_D (v^2))^{1/2}$$

For the transport Eq. (2) and using the equilibrated flux $\theta_h^n \in \mathbf{H}(\text{div}, \Omega)$, the local conservation property is given by the following condition on each control volume $D \in D_h$:

$$\sum_{\gamma \subset \partial D} \int_{\gamma} \theta_h^n \cdot \tilde{\mathbf{n}}_{\gamma} \, d\mathbf{s} = \int_D S \, dx - \int_D \omega R \frac{C_h^n - C_h^{n-1}}{\tau_n} \, dx + \int_D \lambda \omega R C_h^n \, dx. \qquad (8)$$

On $[0, T_f]$, the functions $C_{h,\tau}(t)$ which are a piecewise linear interpolation of $\{C_h^n\}_n$ on each time interval are given by:

$$C_{h,\tau}(t) := C_h^{n-1} + \frac{t-t_{n-1}}{\tau_n} (C_h^n - C_h^{n-1}), \quad \forall t \in (t_{n-1}, t_n], \quad n = 1, \dots, N.$$

The space-time energy norm is defined as

$$\|v\|_X^2 := \int_0^{T_f} \|v\|^2 dt \text{ with } \|v\|^2 := \|D^{\frac{1}{2}} \nabla v\|_{\Omega}^2 + \|(\frac{1}{2} \nabla \cdot \vec{U} - \lambda \omega R)^{\frac{1}{2}} v\|_{\Omega}^2.$$

And then the error estimator for the transport equation is established in the space-time energy norm and formulated as:

$$\|C - C_{h,\tau}\|_X \le \left\{ \sum_{n=1}^{N} \{ (\theta_{sp}^n)^2 + (\eta_{tm}^n)^2 \} \right\}^{1/2} \qquad (9)$$

where the time and space estimators are expressed by

$$\eta_{im}^n := \left\{ \int_{t^n-1}^{t^n} \sum_{D \in D_h^n} \| D^{\frac{1}{2}} \nabla (C_{h,\tau} - C_h^n) - D^{-\frac{1}{2}} \vec{U}_h (C_{h,\tau} - C_h^n) \|_D^2 (t) dt \right\}^{1/2}$$

$$\eta_{sp}^n := \left\{ \tau_n \sum_{D \in D_h^n} m_{D,D}^2 \|S - \omega R \frac{\partial C_{h,\tau}}{\partial t} - \lambda \omega R C_h^n - \nabla \cdot \theta_h^n \|_D^2 \right\}^{1/2} + \left\{ \tau_n \sum_{D \in D_h^n} \|D^{\frac{1}{2}} \nabla C_h^n - D^{-\frac{1}{2}} \vec{U}_h C_h^n + D^{-\frac{1}{2}} \theta_h^n \|_D^2 \right\}^{1/2}$$

with $m_{D,D}$ is a computable constant depending on the diffusion-dispersion tensor $D$, the reaction coefficient $\lambda \omega R$ and the diameter of the dual cell $D$.

The conservative $H(\text{div})$-fluxes $\mathbf{t}_h$ and $\theta_h^n$ verifying respectively (6) and (8) can be derived in the lowest order vector space of Raviart-Thomas-Nédélec RTN$_0(S_h)$ which is a finite-dimensional subspace of $H(\text{div}, \Omega)$. By definition, a vector function from RTN$_0(S_h)$ is piecewise-linear. In order to define such a function, it is sufficient to describe its degrees of freedom associated with the edges (or faces in 3D) of $S_h$. Indeed, considering a dual cell $D \subset D_h$ and $S_D$ the associated local simplicial cell, the fluxes are constructed verifying the local conservation and for all the elements $K \in S_D$ (triangle in 2D and tetrahedron in 3D).
<figure><img src="image_7.png" /><figcaption>Fig. 2. Sub-mesh $S_h$ of the primary mesh $T_h$ and the dual mesh $D_h$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 3. The refinement strategy.</figcaption></figure>

In order to construct such fluxes, the degrees of freedom $u_h$ (discrete approximation of the solution) are defined on the boundary of each element $K \in S_D$. Specifically, the flux across an interface of dual cell is given by $t_h \cdot n_\sigma = -K\nabla u_h \cdot n_\sigma$ where $n_\sigma$ denotes the outward normal vector to the boundary $\sigma$. This formulation ensures that the local conservation principle (6) is satisfied on each element $K \in S_D$. This iterative process propagates across all elements in the domain $S_D$, ensuring that the degrees of freedom on each element are consistent with the conservation principle. The sequential approach ensures compatibility between adjacent elements, avoiding inconsistencies in the fluxes. However, the accuracy and stability of this method depend on the mesh quality and the approximation space used for $u_h$. For more details on the construction of $t_h$ and $\theta_h^\sigma$ see, e.g., [5,6] and the references therein. Drawing from the aforementioned *a posteriori* error estimators, an adaptive mesh refinement procedure, designed to ensure that the relative energy error between the exact and approximate solution will be below a prescribed tolerance $\epsilon$, is presented hereafter.

## 3.2. Adaptive mesh refinement procedure

An overview of the basic time-stepping routine will be given, followed by a detailed description of each time step of the algorithm. The basic adaptive mesh refinement algorithm can be written as the following loop:

solve numerical problem→ calculate the estimator → mark selected elements → refine and/or coarsen each selected element once only
<figure><img src="image_9.png" /><figcaption>Fig. 4. Hierarchical structure used to store the adaptive mesh.</figcaption></figure>

Once the estimators given by (7) and (9) are computed, the next step is to mark the nodes involved in the refinement process. It is important to note that the estimators are defined on the dual cells $D \in D_h$, whereas the refinement process operates on the primary $T_h$. To address this, all elements associated with the dual cell where the estimator value is high (tolerance set by user) are marked for refinement. The refinement strategy for spatial marking can follow either:

• The maximum strategy involves selecting for refinement all elements whose error indicator value exceeds a threshold. This threshold is determined as the maximum value of the error indicator obtained during the simulation, multiplied by a user-defined parameter between 0 and 1. For high values of this parameter (close to 1), the strategy becomes highly selective. Conversely, low values of this parameter (close to 0) lead to the selection of almost all mesh elements.

• The equidistributed Dörfler strategy consists of calculating a threshold defined as the discrete $L^2$ norm of the error indicator multiplied by a reference value. All elements exceeding this threshold are then selected for refinement.

The mesh refinement parameter plays a crucial role in determining which elements of the mesh should be refined, this parameter ranging from 0 to 1 specifies the fraction of the total error to be captured by the selected elements for spatial refinement. Temporal adaptation, on the other hand, is based on the temporal error indicator combined with the CFL condition. Furthermore, the same markers are used for coarsening. The maximum strategy is to select all elements whose error indicator value is below a threshold. This threshold is determined as the maximum error indicator value obtained during simulation, multiplied by a user-defined parameter between 0 and 1. The equidistributed Dörfler strategy consists of calculating a threshold defined as the discrete L2 norm of the error indicator multiplied by a user-defined parameter between 0 and 1. All elements with values of the error indicator below this threshold are then selected for coarsening.

Now that the elements have been marked, they can be refined or coarsened as needed. To achieve this, the refinement algorithm applies the newest vertex bisection method (see [25] and Fig. 3), where elements are bisected into two smaller elements. In two dimensions, each marked triangle from the primary mesh has a specific edge selected for refinement (see Fig. 3(b)), typically the longest edge in the first iteration of the algorithm. This edge, referred to as the refinement edge, is bisected by connecting its midpoint to the opposite vertex, thereby creating two new elements, the aim being to reduce its volume while maintaining a certain regularity in the mesh.

After the marked triangles have been cut in half, usually resulting in a break in conformity, further cuts are propagated to eliminate the orphan nodes. This process, aimed at restoring mesh conformity, is known as completion. Of course, a completion step can produce other orphan nodes and so on. Mitchell [26] and Bansch [27] demonstrate that the completion process will end in a finite number of iterations.

### 3.3. Algorithmic description

Recall that our representation of a triangle is given by element $(t, [0, 1, 2])$, where these indices represent the global indices of the three vertices of triangle $t$. The only constraint on the order of these indices is that the orientation must be counterclockwise. A cyclic permutation of these three indices always represents the same triangle while guaranteeing orientation. This structure was exploited by adopting the following rule: the refinement edge of triangle $t$ is defined as element $(t, [1, 2])$. In other words, when bisecting, element$(t, 0)$ will always be referred to as the newest vertex of triangle $t$.

During the refinement process, orphan nodes are avoided by the refinement algorithm (see Fig. 3(c)). Indeed, bisecting may result in hanging nodes. Neighboring triangles need to be bisected (again following the newest vertex method) to eliminate hanging nodes. This procedure is called completion. For completion, it can observe that if all the edges of the current triangulation are bisected according to the newest vertex method, a conforming triangulation is obtained. Therefore, rather than operating directly on triangles, it is sufficient to bisect a sufficient number of edges to guarantee conformity. This is a recursive process.

To facilitate coarsening and avoid storage of the mesh resulting from refinement, a hierarchical structure based on a binary tree is employed (see [28]). At runtime, the algorithm builds the initial mesh, explicitly establishing connectivity information for each element (e.g., neighbors, opposite vertices), forming so-called macro-elements (see Fig. 4). For child elements, this connectivity
information is inherited from the macro-elements via the hierarchical structure, avoiding the need for explicit storage. During the coarsening process, error indicators are used to mark pairs of triangles for deletion. When such pairs are identified, they are merged, and the hierarchical structure is updated by clearing the relevant level and all its children. This approach ensures efficient management of mesh connectivity and adaptation.

### 3.4. Flow and transport algorithm

The following sequential two-step algorithm was implemented to solve the coupled problem (1)–(2). For an initial conforming triangulation $\mathcal{T}^0$, a prescribed tolerance $\epsilon$, and the coupled parameters $(N_{sp}, N_{tm})$ which correspond to the maximum number of iterations allowed in space and time, the two steps are described in algorithms 1 and 2.

<caption>**Algorithm 1** Resolution of the flow equation</caption>

**Step 1** : compute $P_h$ and $\vec{U}_h$ by solving the Darcy equation Eq. (1) set boundary condition on $P_h$, $\mathcal{T}_h := \mathcal{T}^0$, $n_{sp} = 0$
loop adaptMesh
    solve $P_h = \text{DARCY}(\mathcal{T}_h)$, tol = $\epsilon ||P_h||_*$
    estimate $\eta_D = \eta_{R,D} + \eta_{DF,D}$, $\eta = \{\sum_{D \in D_h} \eta_D^2\}^{1/2}$
    If $(\eta \ge tol$ and $n_{sp} < N_{sp})$ Then
        ! refine the mesh and resolve
        mark($\mathcal{T}_h$, $\eta_D$), $\mathcal{T}_h := \text{refine}(\mathcal{T}_h)$, $n_{sp} = n_{sp} + 1$
        continue loop adaptMesh
    Else exit loop adaptMesh

**Remark.** The function `mark` ($\mathcal{T}_h$, $\eta_D$) introduced in this algorithm defines a set of elements that will be refined. At each refinement stage, the flow is recalculated and then re-evaluated the estimator to determine the regions to be refined.

<caption>**Algorithm 2** Resolution of the transport equation using the flow obtained in the first step</caption>

**Step 2** : compute $C_h^n$ by solving the transport equation Eq. (2). In order to equilibrate the space and time estimators in each iteration of mesh adaptation, according to [6], $\eta_{sp}^n, \eta_{tm}^n < tol := \frac{1}{\sqrt{2}}\epsilon ||C_{h,\tau}||_{X(t^{n-1},t^n)}$ is saved.

set $\mathcal{T}_h^0 := \mathcal{T}_h$ the final mesh obtained by step 1, compute $\vec{U}_h$.
set and interpolate initial, boundary condition $C_h^0$ on $\mathcal{T}_h^0$
init $\tau_1$, set $t^0 = 0$, $t^1 = t^0 + \tau_1$, $\mathcal{T}_h^1 := \mathcal{T}_h^0$, $n = 1$, $n_{tm} = 0$
Loop in time: While ($t^n \le T_f$)
    Loop adapt timestep
        set $n_{sp} = 0$
        Loop adapt mesh
            solve $C_h^n = \text{TRANSPORT}(C_h^{n-1}, \tau_n, \mathcal{T}_h^n)$
            estimate $\eta_{sp}^n, \eta_{tm}^n$, tol := $\frac{1}{\sqrt{2}}\epsilon ||C_{h,\tau}||_{X(t^{n-1},t^n)}$
            If $(\eta_{sp}^n \ge tol$ and $n_{sp} < N_{sp})$ Then
                ! refine/coarsen the mesh and resolve
                mark($\mathcal{T}_h^n$, $\eta_{sp,D}^n$), $\mathcal{T}_h^n := \text{refine\_coarsen}(\mathcal{T}_h^n)$
                interpolate $(C_h^{n-1}, C_h^n, \vec{U}_h)$ on $\mathcal{T}_h^n$, $n_{sp} = n_{sp} + 1$
                continue loop adapt mesh
            Else exit loop adapt mesh
        If $(\eta_{tm}^n \ge tol$ and $n_{tm} < N_{tm})$ Then
            ! refine the time step and resolve
            set $\tau_n = \tau_n/2$, $t^n = t^{n-1} + \tau_n$, $n_{tm} = n_{tm} + 1$
            continue loop adapt timestep
        Else ! go forward and coarsen the time step if needed
            If $(\eta_{tm}^n < 0.5 * tol)$ set $\tau_n = 2 * \tau_n$
            set $n = n + 1$, $t^n = t^{n-1} + \tau_n$, $n_{tm} = 0$
            continue loop in time
        End loop in time

**Remark.** The function `mark` ($\mathcal{T}_h^n$, $\eta_{sp,D}^n$) introduced in this algorithm defines a set of elements that will be refined or coarsened. It should be noted that the local estimators $\eta_D := \eta_{sp,D}^n$ are defined on dual cells $D \in D_h$, while the refinement process will be performed on the primary mesh $\mathcal{T}_h$. One possible approach is to mark all elements $T \in \mathcal{T}_h$ associated with the dual cell $D$ where the value of the corresponding local estimator is significant or negligible compared to a defined strategy in order to refine or coarsen. To resolve the transport equation, the flow is obtained by interpolation from the initial mesh introduced for the calculation.
<table><caption>**Table 1**
Physical parameters for the two porous materials.</caption><thead><tr><th>Media</th><th>K (m.y<sup>-1</sup>)</th><th>ω (%)</th><th>α<sub>L</sub> (m)</th><th>α<sub>T</sub> (m)</th><th>D (m<sup>2</sup>.y<sup>-1</sup>)</th></tr></thead><tbody><tr><td>black</td><td>10<sup>-5</sup></td><td>20</td><td>5</td><td>0.5</td><td>10<sup>-5</sup></td></tr><tr><td>gray</td><td>1.0</td><td>20</td><td>5</td><td>0.5</td><td>10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Heterogeneous domain.</figcaption></figure>

# 4. Numerical results

In this section, the adaptive methodology is illustrated with two example problems to validate the efficiency, accuracy, and portability of the method for industrial applications. These tests include: (1) a heterogeneous medium with two different permeability distributions (values ranging from $10^{-5}$ m.y$^{-1}$ to 1 m.y$^{-1}$), and (2) a realistic heterogeneous fractured case. These tests consider two-dimensional miscible flow problems and the evolution of radionuclide migration from the waste packages through a heterogeneous disposal, where large permeability variations are allowed with different source terms. The used computer for the following simulation is an Intel(R) Xeon 2.40 GHz with 32 GB RAM.

## 4.1. Test problem 1: heterogeneous basic case

For this example, a heterogeneous medium is defined on the domain $\Omega = (0, 1) \times (0, 1)$ with a discontinuous porous material of low permeability equal to $10^{-5}$ m.y$^{-1}$ (black color) inside a high permeability porous material equal to 1 m.y$^{-1}$ (gray color), as shown in Fig. 5. The dispersivity is assumed to be similar for both materials. The values of the physical parameters are given in Table 1. The boundary conditions are illustrated in Fig. 5(a): $p = 0$ m in $\Gamma_{in}$, $p = 1$ m in $\Gamma_{out}$ (red corners), and no-flow conditions are imposed on the other boundaries. A tracer is assumed to be instantaneously released during the first time step of the simulation from the source term location, which corresponds to a tracer injection ($10^3 Bq$) at the lower-left of the domain (red color on $\Gamma_{in}$). The transport equation is solved with an adapted time step size for a time period of 10 years. The primary grid is presented in Fig. 5(b).

In the literature, it is known that the jump in parameters at the inner boundaries may cause problems if the mesh is not adapted to attenuate this jump, especially near the domain's corners. Therefore, in this example, the behavior of the adaptive algorithm near the inner boundaries is studied. For this test problem, the previously described adaptive strategy is used. The error indicators are computed for each dual cell $D$ at each time step and the elements $T$ of the dual cell with the largest error indicator values are refined, and the elements of the dual cell with the smallest error indicator values are coarsened. As expected, for the water flow simulation, the mesh is refined mainly near the interface between the two porous media where a high contrast of permeability is observed (see hydraulic head profiles in Fig. 6(c)).

The Darcy velocities are represented in Figs. 7(a) and 7(b) with the static mesh and in Fig. 7(c) with the refined mesh. Their values increase in the refined mesh area. This phenomenon is observed only in the upper right corner (denoted as $\Gamma_{out}$ in Fig. 5(a)), and may be due to the effect of the Dirichlet boundary condition and the mesh size. Indeed, a difference in the maximum value
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Hydraulic head profiles simulated by MELODIE with different strategy.</figcaption></figure>

is observed between the two static meshes with different discretizations. The maximum value is 0.98 m/y for a simulation carried out with a mesh discretization of $h = 0.1$ (see Fig. 7(a)), while it is 2.11 m/y for a simulation carried out with a much finer mesh discretization of $h = 0.002$ (see Fig. 7(b)). For the simulation obtained using the adaptive mesh strategy, the maximum value of the Darcy velocities is 2.729 m/y.

As can be observed in Figs. 7(b) and 7(c), these two models applied to the same test example can be considered equivalent, even though a small difference between the two simulations of the Darcy velocities can be noted. However, it is clear from Table 2 that the approach with the adaptive mesh required less CPU time than the simulation done with a static mesh with very fine discretization.

Fig. 8 illustrates the evolution of the spatial estimators for the flow equation at the interfaces. This estimator detects the error caused by the difference in permeability between the two materials. As the mesh adaptation iterations progress, the value of the error estimator decreases. This numerical test shows that the proposed adaptive strategy is effective for guiding a space mesh refinement strategy using this class of estimator. Furthermore, it can be observed that a large number of iterations is needed to improve the solution, which can lead to a long simulation time.

To reduce the number of iterations and the CPU time needed to complete the simulation, a multi-level refinement process is considered (see [23,29,30]). Indeed, the basic idea would be to reduce the overall value of the error estimator as quickly as possible
<figure><img src="image_12.png" /><figcaption>Fig. 7. Darcy velocities.</figcaption></figure>

<table><caption>**Table 2**
CPU time for the flow simulation and meshes size.</caption><thead><tr><td>Mesh</td><td>Nodes</td><td>Elements</td><td>cpu time (s)</td></tr></thead><tbody><tr><td>static mesh h = 0.1</td><td>161</td><td>280</td><td>8 × 10<sup>-3</sup></td></tr><tr><td>static mesh h = 0.002</td><td>330 752</td><td>659 502</td><td>3.896 × 10<sup>+3</sup></td></tr><tr><td>adaptative mesh</td><td>2307</td><td>4542</td><td>2.96 × 10<sup>-1</sup></td></tr></tbody></table>

(i.e. with fewer iterations involving expensive system solving). The basic refinement strategy consists in solving the system, then cutting out the marked elements once and repeating this process. The idea of multi-level refinement is to improve approximation in regions of the domain where the error is significant, while minimizing computational cost. In effect, the idea is to add a piece of information that would consist in determining the number of refinements the elements will have to undergo as a function of the magnitude of the estimator value. By introducing this trick, the optimal strategy is attained while reducing computation times indeed, at each step each element of the dual cell can be refined more than once depending on the magnitude of the estimator. The result is presented in Fig. 9. As expected, the number of iterations decreases according to the level of refinement. However, the
<figure><img src="image_13.png" /><figcaption>Fig. 8. Spatial distribution of error estimators at the end of simulation.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 9. Multi-level refinement of the flow equation.</figcaption></figure>

density of the mesh increases. Levels 2 and 3 seem to be a good compromise between the number of iterations of the refinement strategy and the density of the mesh. Level 2 was used for the numerical simulation in this article. This choice was guided by the observation of solution oscillation when the refinement level is equal to 1 and by the need to improve calculation speed in the simulation while maintaining an accurate solution.

For the purpose of testing the developed method, the location of the tracer source was chosen to favor transport along the heterogeneities. The propagation of mesh refinement and the corresponding activities are displayed in Fig. 12. This simulation was also carried out with a static mesh in space and time (see Figs. 10 and 11). The simulation results show that the contaminant moves from the source located in the bottom left corner of the domain to the outside of the domain. A difference is observed between the simulations carried out with coarse and fine static meshes; indeed, this difference is due to the phenomenon of numerical diffusion known in the finite volume scheme.
<figure><img src="image_15.png" /><figcaption>Fig. 10. 2D tracer concentration profiles ($Bq/m^3$) obtained using a static mesh $h = 0.1$.</figcaption></figure>

For the simulation with an adaptive strategy, the mesh refinement follows the displacement of the plume in the corner of the domain at the beginning of the simulation (see Fig. 12), and over time, it tracks the contour lines of the medium with low permeability, as shown in Fig. 13. In fact, diffusion is the dominant transfer mechanism through the low permeability material, whereas both advective and diffusive transfers occur through the high permeability zone. The simulation without a refinement strategy, shown in Fig. 10, does not allow this phenomenon to be observed due, on the one hand, to a too coarse mesh and, on the other hand, to the numerical diffusion of the FVFE scheme. However, a difference in the maximum value of the activity is observed between the simulation performed with a static mesh, where the discretization is $h = 0.002$ (see Fig. 11), and the simulation performed with the adaptive strategy (see Fig. 13). This difference is certainly due to the temporal scheme implemented in this algorithm. Indeed, the latter is based on a semi-implicit Euler approach, which requires knowledge of the concentration over three successive iterations. This scheme likely accumulates errors in this model, as the computation of the various fields on the refined meshes is obtained by interpolation of a non-updated solution for each time iteration.

Fig. 14 illustrates the evolution of the spatial estimators for the transport equation. As can be seen, this estimator, plotted for the last iteration in time, detects the error caused by the difference in permeability between the two materials at the interfaces, as well as the error following the displacement of the concentration front, which may occur more by diffusion than by convection at the end of this simulation. Moreover, less numerical diffusion is observed in the transport simulation, ensuring the quality of the strategy implemented for these simulations. This numerical test shows that the proposed adaptive strategy for the transport model is effective for guiding a mesh refinement/coarsening strategy using this class of estimator. The CPU time needed for the approach
<figure><img src="image_16.png" /><figcaption>Fig. 11. 2D tracer concentration profiles (Bq/m³) obtained using a static mesh h = 0.002.</figcaption></figure>

<caption>**Table 3**</caption>

CPU time for the transport simulation.

<table><thead><tr><th>mesh</th><th>cpu time (s)</th></tr></thead><tbody><tr><td>static mesh h = 0.1</td><td>4.4 × 10<sup>-2</sup></td></tr><tr><td>static mesh h = 0.002</td><td>2.09 × 10<sup>+3</sup></td></tr><tr><td>adaptative mesh</td><td>1.02 × 10<sup>+3</sup></td></tr></tbody></table>

with the adaptive mesh remains lower than the CPU time for a simulation done with a static mesh with a very small discretization, about half as much (see Table 3).

The maximum and minimum values of the activities in Bq/m³ during the transport simulation are presented in Fig. 15. Notice that for the maximum values of the activities (see Fig. 15(a)), overall the simulation can be considered equivalent even though a small difference appeared at the end of the simulation between the adaptive strategies and the static meshes. This difference can be attributed to the semi-implicit time scheme, which overestimates the activities. For the minimum values of the activities (see Fig. 15(b)), the simulations done with static meshes show that the numerical solution of the concentration is not strictly positive, corresponding to a violation of the minimum principle. The minimum (or maximum) principle is an important property of linear advection-diffusion type problems. In numerical simulations, the violation of the discrete minimum (or maximum) principle leads
<figure><img src="image_17.png" /><figcaption>Fig. 12. Mesh and 2D tracer concentration profiles with adaptive strategy at 3 × 10⁻⁷ year.</figcaption></figure>

<table><caption>**Table 4**
Numerical parameters for transport simulation.</caption><thead><tr><th>Parameters</th><th>Strategy 1</th><th>Strategy 2</th></tr></thead><tbody><tr><td>max iteration on space</td><td>8</td><td>20</td></tr><tr><td>ε</td><td>0.1</td><td>0.05</td></tr><tr><td>mark (τ<sub>h</sub>, η<sub>D</sub>)</td><td>maximum (MS)</td><td>maximum (MS)</td></tr><tr><td>MS refine</td><td>0.75</td><td>0.25</td></tr><tr><td>MS coarsen</td><td>0.55</td><td>0.1</td></tr><tr><td>max iteration on time</td><td>4</td><td>20</td></tr><tr><td>η<sub>m</sub></td><td>0.1</td><td>0.1</td></tr></tbody></table>

to non-physical solutions with numerical artifacts. Modern methods, such as advanced finite volume schemes, have solved many of these problems. As can be seen in Fig. 15(b), combining a semi-implicit scheme in time and a globally refined mesh is not enough to guarantee this principle. This figure also illustrates the simulation realized with the adaptive methodology using two different sets of numerical parameters called strategy 1 and strategy 2 (see Table 4 for more details on these numerical parameters). This numerical test shows overall better behavior with a decrease in the absolute value of the minimum concentration. Indeed, these oscillations are less amplified by adaptive mesh refinement. For strategy 2, these oscillations are almost non-existent, and a well-chosen stopping criterion could be the solution to respect this principle.

It should be noted that some numerical artifacts appear in the simulations, particularly in the form of slight tracer leakage into low-permeability zones and local negative concentrations (Fig. 15(b)). These phenomena are consistent with the known limitations of the FVFE scheme, especially in highly heterogeneous media, where the combination of explicit upwinding and finite element diffusion can lead to non-physical behavior near sharp concentration fronts. Despite these issues, the overall dynamics of the transport process are well captured. Furthermore, the use of an adaptive mesh refinement strategy significantly reduces the extent of such artifacts. As shown in the results, strategy 2 in particular leads to a substantial decrease in the minimum negative concentration, approaching physical admissibility. These observations validate the interest of adaptive control and numerical parameter tuning in improving the robustness and physical relevance of the simulations.

## 4.2. Test problem 2: realistic faulted basin

This test case is chosen to illustrate the capability of the proposed method in capturing discrepancies in geological media characteristics relevant to the geological disposal of radionuclides. The selected scenario uses field data representative of real hydrogeological media. The entire stratigraphic sequence consists of seven hydrogeological media intersected by a sub-vertical
<figure><img src="image_18.png" /><figcaption>Fig. 13. Mesh and 2D tracer concentration profiles with adaptive strategy.</figcaption></figure>

fault. The vertical 2D cross-section dimensions are approximately 18,000 meters long and 2,000 meters deep. Four source terms are considered in this test: two located in a thin layer on the left side of medium 4 at approximately 400 meters depth, and two located in a thicker layer on the right side of medium 4 at approximately 650 meters depth. These two layers (host rock) have the same hydrological properties but differ in thickness due to the presence of the fault. This medium is assumed to represent a sedimentary, low permeability rock with an associated hydraulic conductivity of about $3.15 \times 10^{-3} \text{ m.y}^{-1}$. The hydraulic conductivity of upper and lower layers in contact with the host rock are respectively about $3.15 \times 10^{-1} \text{ m.y}^{-1}$ and $3.15 \times 10^2 \text{ m.y}^{-1}$. The fault is assumed to have a conductive feature with a hydraulic conductivity is $3.15 \times 10^2 \text{ m.y}^{-1}$ (see Fig. 16).

The flow simulation is conducted under the condition of hydraulic heads equal to the topographic level at the model's surface (Dirichlet boundary condition). Hydraulic heads are also prescribed on the left and right sides of the domain. The bottom boundary is assumed to be impermeable. For this test case, iodine-129 (I-129), a long-lived radionuclide that is soluble and non-sorbed, was chosen as an example for release and transport. I-129 is assumed to be instantaneously released from four source locations, each with an initial inventory of $10^9 \text{ Bq}$ distributed equally among them. All physical parameter values chosen for this example are summarized in Table 5. The transport equation is solved using an adapted time step size over a period of 10,000 years. The FVFE method operates
<figure><img src="image_19.png" /><figcaption>Fig. 14. Spatial distribution of error estimators at the end of simulation (10 years).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 15. Evolution of maximum and minimum activities.</figcaption></figure>

on triangular elements and requires the construction of a dual mesh. In this modeling case, the dual mesh is based on centroids. The global mesh consists of 11,614 triangular elements and 5,982 nodes distributed across eight different media.

### 4.2.1. Flow results

The flow simulation is carried out with the initial mesh presented in Fig. 17(c). With this static mesh, the hydraulic head fields are presented in Fig. 17(b), and the Darcy velocities are represented in Fig. 17(c).

As expected, the boundary conditions imposed on the left and right sides of the domain generate a flow towards the fault. Moreover, due to the hydraulic characteristics of the fault, low hydraulic heads are transmitted from the topographic surface to the bottom of the domain through the fault. The heterogeneous flow obtained by this model, due to contrasted hydraulic conductivities,
<figure><img src="image_21.png" /><figcaption>Fig. 16. Realistic case of a faulted layered basin (permeability distribution).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 17. Flow simulation with a static mesh.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>(a) Final mesh</figcaption></figure>

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

<figure><img src="image_25.png" /><figcaption>(b) Associated estimator</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(a) Darcy velocities (m/y)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) A zoom on the junction (m/y)</figcaption></figure>

<caption>Fig. 19. Darcy velocities associated with the flow simulation (adaptive process).</caption>
<figure><img src="image_28.png" /><figcaption>Fig. 20. Adaptative mesh.</figcaption></figure>

<caption>Table 5</caption>

Physical parameters of the different layers and the fault (with the source locations marked in green)

<table><thead><tr><th>Media</th><th>K (m.y<sup>-1</sup>)</th><th>ω (%)</th><th>α<sub>L</sub> (m)</th><th>α<sub>T</sub> (m)</th><th>D (m<sup>2</sup>.y<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>3.15 × 10<sup>-3</sup></td><td>30</td><td>1.0</td><td>1.0</td><td>10<sup>-2</sup></td></tr><tr><td>2</td><td>3.15 × 10<sup>-3</sup></td><td>15</td><td>0.5</td><td>0.5</td><td>5 × 10<sup>-2</sup></td></tr><tr><td>3</td><td>3.15 × 10<sup>-2</sup></td><td>17</td><td>10.0</td><td>1.0</td><td>5 × 10<sup>-2</sup></td></tr><tr><td>4</td><td>3.15 × 10<sup>-3</sup></td><td>30</td><td>8.0</td><td>4.0</td><td>1</td></tr><tr><td>5</td><td>3.15 × 10<sup>-1</sup></td><td>20</td><td>5.0</td><td>0.5</td><td>1</td></tr><tr><td>6</td><td>3.15 × 10<sup>-3</sup></td><td>20</td><td>12.0</td><td>6.0</td><td>5 × 10<sup>-1</sup></td></tr><tr><td>7</td><td>3.15 × 10<sup>-1</sup></td><td>15</td><td>0.5</td><td>0.5</td><td>5 × 10<sup>-1</sup></td></tr><tr><td>fault (8)</td><td>3.15 × 10<sup>-2</sup></td><td>20</td><td>12.0</td><td>6.0</td><td>5 × 10<sup>-2</sup></td></tr></tbody></table>

is consistent with the higher velocities observed in underlying aquifers and the fault, and the lower velocities present in other media. The behavior of the final mesh and its corresponding a posteriori error estimate obtained from the last simulation are presented in Fig. 18(b). The fixed mesh used for this simulation is shown in Fig. 17(a). The error detected by the estimator is located around the interface of layer 3, which has high permeability, and also around the interface of the fault. This result demonstrates that the indicator implemented in this work is efficient for industrial applications with a complicated model and justifies the use of the mesh refinement algorithm presented in algorithm 1.

The results of the Darcy velocities associated with the flow simulation using an adaptive process are summarized in Fig. 19(a). As expected, significant Darcy velocities are obtained in this example around the underlying aquifers of the host rock and fault. This result is a consequence of a very refined mesh at the junction between media with high hydraulic conductivity and the fault. It can
<figure><img src="image_29.png" /><figcaption>Fig. 21. I-129 concentration 2D profiles obtained at 30 years (log scale) $Bq/m^3$.</figcaption></figure>

<table><caption>**Table 6**

Mesh evolution.</caption><thead><tr><td>Simulation</td><td>Triangles</td><td>Nodes</td></tr></thead><tbody><tr><td>static mesh</td><td>11,614</td><td>5,982</td></tr><tr><td>mesh at 0 year</td><td>60,804</td><td>30,679</td></tr><tr><td>mesh at 30 years</td><td>62,908</td><td>31,731</td></tr><tr><td>mesh at 500 years</td><td>88,753</td><td>44,826</td></tr><tr><td>mesh at 10,000 years</td><td>1,361,747</td><td>683,053</td></tr></tbody></table>

also be associated with minor local perturbation at this junction. Indeed, in highly heterogeneous zones (such as faults), pressure gradients or velocities can intensify, leading to behaviors that depart from Darcy's hypotheses.

Furthermore, the number of elements obtained during the simulation (presented in Fig. 18(a)) is significantly reduced using the adaptive refinement strategy compared to using a very fine static mesh. On average, the number of cells is reduced by more than 90%, which is a significant gain. Indeed, to obtain the same quality of results, a uniformly refined mesh with 2,973,184 elements would be required.

### 4.2.2. Transport results

The propagation of the mesh refinement and the corresponding concentrations at 30, 500, and 10,000 years are displayed in Figs. 20, 21, 22 and 23, respectively. Details on mesh evolution are provided in Table 6. The simulation results show that most of the contaminant moves from the different source locations to the fault. In these figures, the refinement follows the contaminant plume as time evolves, and the de-refinement process is carried out in areas free from the contaminant.

At the beginning of the simulation, for the static mesh the contaminant plumes seems to reach the surrounding aquifers very fast. However, the contaminant plumes in the right and left part of the domain had to be diffusive and concentrated in the host rock due to its higher thickness combined and its low hydraulic conductivity, which delays migration to the surrounding aquifers.
<figure><img src="image_30.png" /><figcaption>Fig. 22. I-129 concentration 2D profiles obtained at 500 years (log scale) $Bq/m^3$.</figcaption></figure>

Thus, the propagation phenomenon surrounding aquifers observed in Fig. 21(a) is due to the non-refined mesh and the numerical diffusion of the FVFE scheme. Fig. 21(b) shows the evolution of the contaminant plumes with an adaptive mesh refinement strategy. It can be remarked that the refinement follows the source term locations as time evolves. The high concentrations observed in the left part of the domain are due to the better mesh obtained with the refinement process. Finally, due to its location, the leftmost source term seems to migrate into the aquifer with low permeability and join the fracture more quickly.

At the beginning of the simulation, with the static mesh, the contaminant plumes seem to reach the surrounding aquifers very quickly. However, the contaminant plumes in the right and left parts of the domain should be diffusive and concentrated in the host rock due to its greater thickness combined with its low hydraulic conductivity, which delays migration to the surrounding aquifers. Thus, the propagation phenomenon observed in the surrounding aquifers in Fig. 21(a) is due to the non-refined mesh and the numerical diffusion of the FVFE scheme. Fig. 21(b) shows the evolution of the contaminant plumes with an adaptive mesh refinement strategy. It can be noted that the refinement follows the source term locations as time evolves. The high concentrations observed in the left part of the domain are due to the better mesh obtained with the refinement process. Finally, due to its location, the leftmost source term seems to migrate into the low permeability aquifer and join the fracture more quickly.

After 500 years, with the static mesh (Fig. 22(a)), the radionuclide moves upward by diffusion and advection, and the fault behaves as a drain. Indeed, the radionuclide is transported along the underlying aquifer and trapped by the fault when it reaches the topographic surface. Due to the coarse discretization, more numerical diffusion in the concentration front displacement is observed. For the refined mesh (Fig. 22(b)), this phenomenon is not observed for all the source terms. In fact, due to the refined mesh strategy, the contaminant plumes are concentrated around the initial locations of the source terms. Moreover, these contaminants do not propagate towards the fault with the numerical diffusion phenomenon.

After 10,000 years, the contaminant plumes progressively move towards the outlets at the topographic surface, following the direction of flow with the refinement process (Fig. 23(b)). However, with the static mesh, the contaminant plumes seem to disperse more into other aquifers (Fig. 23(a)).

As observed in this realistic case, the use of a static mesh leads to noticeable numerical diffusion, especially in the early stages of the simulation. The contaminant plumes tend to spread too quickly into the surrounding aquifers, and their shapes are smoothed,
<figure><img src="image_31.png" /><figcaption>Fig. 23. I-129 concentration 2D profiles obtained at 10,000 years (log scale) Bq/m³.</figcaption></figure>

which is a known artifact of the FVFE scheme when applied to coarse, non-adaptive meshes. These inaccuracies are particularly pronounced near the fault and in low-permeability regions, where strong contrasts in hydraulic properties exist. By contrast, the adaptive mesh refinement strategy significantly improves the resolution of the concentration fronts, reduces non-physical dispersion, and ensures better consistency with the expected behavior of transport in low-permeability host rocks. These observations confirm the limitations and strengths of the FVFE method in heterogeneous settings and reinforce the importance of coupling it with adaptive mesh control to mitigate numerical artifacts.

**5. Conclusions**

The reliability of the method implemented in the MELODIE software to calculate flow and transport is supported by several calculation cases that are presented and discussed. Moreover, this paper proposes the implementation of an adaptive mesh refinement algorithm for modeling flow and transport in porous media, specifically in relation to the simulation of radionuclide migration around geological disposals. The usual system of equations is discretized using a vertex-centered finite volume method combined with a finite element method.

The adaptive mesh refinement strategy implemented in MELODIE was developed with a posteriori error estimators based on $H(\text{div})$-conforming flux reconstruction. This methodology allows a considerable improvement in computational efficiency, as illustrated by the two numerical examples in two dimensions presented in this study. These examples demonstrate that the proposed adaptive mesh refinement is appropriate for solving the flow and transport of radionuclides in heterogeneous porous media in a robust and accurate manner for various physically relevant phenomena. As a future perspective, it propose exploring comparisons with methods like TPFA in order to eliminate violations of the maximum principle in the numerical scheme implemented.

Additionally, the use of mesh adaptation for the simulation of these problems has significantly increased the accuracy and precision of the results while achieving a very significant reduction in simulation time. The objective of future work will be to extend this approach to more realistic three-dimensional problems and to clarify the stopping criteria for adaptive mesh refinement implemented in the MELODIE software.
CRediT authorship contribution statement

Hanen Amor: Writing – original draft, Software, Methodology, Conceptualization. Fayssal Benkhaldoun: Conceptualization.

Acknowledgments

We thank Pascal Omnes (CEA) and Emmanuel Audusse (University Paris 13) for their constructive suggestions during the
development of this work.

References

[1] G. Mathieu, M. Dymitrowska, M. Bourgeois, Modeling of radionuclide transport through repository components using finite volume finite element and multidomain methods, Phys. Chem. the Earth, Parts A/ B 33 (2008) http://dx.doi.org/10.1016/j.pce.2008.10.041.

[2] H. Amor, M. Bourgeois, G. Mathieu, Benchmark 3D: a linear finite element solver : 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids, Springer Berl. Heidelb. : Finite Vol. Complex Appl. VI Probl. Perspect. 4 (2011) 895–930.

[3] B. Amaziane, M. Bourgeois, M. El Fatini, Adaptive Mesh Refinement for a Finite Volume Method for Flow and Transport of Radionuclides in Heterogeneous Porous Media, Oil Gas Sci. Technol. - Rev. D' IFP Energies Nouv. 69 (4) (2014) 687–699, http://dx.doi.org/10.2516/ogst/2013176, URL https://hal.archives-ouvertes.fr/hal-01933397.

[4] A. Ern, M. Vohralík, A posteriori error estimation based on potential and flux reconstruction for the heat equation, SIAM J. Numer. Anal. 48 (2010) http://dx.doi.org/10.1137/090759008.

[5] M. Vohralík, Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients, J. Sci. Comput. 46 (2011) http://dx.doi.org/10.1007/s10915-010-9410-1.

[6] D. Hilhorst, M. Vohralík, A posteriori error estimates for combined finite volume–finite element discretizations of reactive transport equations on nonmatching grids, Comput. Methods Appl. Mech. Engrg. 200 (2011) 597–613, http://dx.doi.org/10.1016/j.cma.2010.08.017.

[7] S. Yousef, M. Vohralík, A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows, Comput. Methods Appl. Mech. Engrg. 331 (2017) http://dx.doi.org/10.1016/j.cma.2017.11.027.

[8] D.A. Di Pietro, M. Vohralík, S. Yousef, A posteriori error estimates with application of adaptive mesh refinement for thermal multiphase compositional flows in porous media, Comput. Math. Appl. (2013) URL https://apisemanticscholar.org/CorpusID:15662357.

[9] D. Di Pietro, M. Vohralík, A review of recent advances in discretization methods, a posteriori error analysis, and adaptive algorithms for numerical modeling in geosciences, Oil Gas Sci. Technol. 69 (2014) 701–729, http://dx.doi.org/10.2516/ogst/2013158.

[10] B. Amaziane, M. Ossmani, C. Serres, Numerical modeling of the flow and transport of radionuclides in heterogeneous porous media, Comput. Geosci. 12 (2008) 437–449, http://dx.doi.org/10.1007/s10596-008-9083-0.

[11] G. Chavent, J. Jaffré, Mathematical Models and Finite Elements in Reservoir Simulation, North Holland, 1986.

[12] D. Di Pietro, E. Flauraud, M. Vohralík, S. Yousef, A posteriori error estimates, stopping criteria, and adaptivity for multiphase compositional Darcy flows in porous media, J. Comput. Phys. 276 (2014) 163–187, http://dx.doi.org/10.1016/j.jcp.2014.06.061.

[13] A. Bourgeat, M. Kern, S. Schumacher, J. Talandier, The COUPLEX test cases: Nuclear waste disposal simulation: Simulation of transport around a nuclear waste disposal site: The COUPLEX test cases, in: A. Bourgeat, M. Kern (Eds.), Comput. Geosci. 8 (2004) http://dx.doi.org/10.1023/B:COMG.0000035073.03009.5d.

[14] Oecd, Safety of Geological Disposal of High-level and Long-lived Radioactive Waste in France: An International Peer Review of the "Dossier 2005 Argile" Concerning Disposal in the Callovo-Oxfordian Formation, OECD Pap. 6 (7) (2006) 1–76, URL https://ideas.repec.org/a/oec/packaa/5l9hnjnvhj20.html.

[15] P. Knabner, S. Bitterlich, R.I. Teran, A. Prechtel, E. Schneid, Influence of surfactants on spreading of contaminants and soil remediation, Math. — Key Technol. Futur. (2003) 152–161.

[16] L. Amir, M. Kern, A global method for coupling transport with chemistry in heterogeneous porous media, Comput. Geosci. 14 (3) (2010) 465–481, http://dx.doi.org/10.1007/s10596-009-9162-x, URL https://hal.inria.fr/inria-00442243.

[17] J. Bear, Dynamics of fluids in porous media., J. Fluid Mech. 61 (1) (1973) 206–208, http://dx.doi.org/10.1017/S0022112073210662.

[18] R. Eymard, T. Gallouët, R. Herbin, J.-C. Latché, A convergent Finite Element-Finite Volume scheme for the compressible Stokes problem, Part II – the isentropic case, Math. Comp. 79 (270) (2010) 649–675, http://dx.doi.org/10.1090/S0025-5718-09-02310-2, URL https://hal.science/hal-00359140.

[19] E. Ahusborde, B. Amaziane, M. El Ossmani, A fully implicit finite volume scheme for single phase flow with reactive transport in porous media, in: Reactive Transport Modeling Workshop, Paris, France, 2018.

[20] R. Eymard, T. Gallouet, R. Herbin, J.-C. Latche, A convergent finite element-finite volume scheme for the compressible Stokes problem, part II – the isentropic case, Math. Comp. 79 (270) (2010) 649–675, http://dx.doi.org/10.1090/S0025-5718-09-02310-2.

[21] I. Ben Gharbia, J. Dabaghi, V. Martin, M. Vohralík, A posteriori error estimates and adaptive stopping criteria for a compositional two-phase flow with nonlinear complementarity constraints, Comput. Geosci. (2019) http://dx.doi.org/10.1007/s10596-019-09909-5, URL https://hal.science/hal-01919067.

[22] A. Ern, I. Smears, M. Vohralík, Equilibrated flux a posteriori error estimates in L²(H¹)-norms for high-order discretizations of parabolic problems, IMA J. Numer. Anal. 39 (3) (2018) 1158–1179, http://dx.doi.org/10.1093/imanum/dry035.

[23] P. Daniel, A. Ern, M. Vohralík, An adaptive hp-refinement strategy with inexact solvers and computable guaranteed bound on the error reduction factor, Comput. Methods Appl. Mech. Engrg. 359 (2020) 112607, http://dx.doi.org/10.1016/j.camwa.2018.05.034.

[24] H. Amor, F. Benkhaldoun, T. Ghoudi, I. Kissami, M. Seaid, New criteria for mesh adaptation in finite volume simulation of planar ionization wavefront propagation, in: Finite Volumes for Complex Applications VIII - Hyperbolic, Elliptic and Parabolic Problems, Springer International Publishing, 2017, pp. 547–555.

[25] I. Kossaczky, A recursive approach to local mesh refinement in two and three dimensions, J. Comput. Appl. Math. 55 (3) (1994) 275–288, http://dx.doi.org/10.1016/0377-0427(94)90034-5, URL http://dx.doi.org/10.1016/0377-0427(94)90034-5.

[26] W.F. Mitchell, Unified multilevel adaptive finite element methods for elliptic problems. URL https://www.osti.gov/biblio/7230871.

[27] E. Bänsch, Local mesh refinement in 2 and 3 dimensions, IMPACT Comput. Sci. Eng. 3 (1991) 181–191.

[28] A. Schmidt, K. Siebert, Design of adaptive finite element software: The finite element toolbox ALBERTA, Springer, 2005.

[29] I. Kissami, C. Cérin, F. Benkhaldoun, G. Scarella, Towards Parallel CFD Computation for the ADAPT Framework, vol. 10048, 2016, pp. 374–387.

[30] E. Biotteau, A. Combescure, A. Gravouil, A. Lubrecht, Stratégie multigrille et méthode de raffinement automatique à précision contrôlée pour la modélisation des structures élastoplastiques 3D en dynamique transitoire non-linéaire, in: 10e Colloque National en Calcul des Structures, Giens, France, URL https://hal.archives-ouvertes.fr/hal-00731723.