# Multiphase flow and reactive transport benchmark for radioactive waste disposal

Javier Samper, Alba Mon, Etienne Ahusborde, Han Yu, Asta Narkuniene,
Milan Hokr, Luis Montenegro, Brahim Amaziane, Mustapha El Ossmani,
Tianfu Xu, et al.

► To cite this version:

Javier Samper, Alba Mon, Etienne Ahusborde, Han Yu, Asta Narkuniene, et al. Multiphase flow and reactive transport benchmark for radioactive waste disposal. Environmental Earth Sciences, 2024, 83 (22), pp.619. 10.1007/s12665-024-11887-6. hal-04771061

HAL Id: hal-04771061

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

Submitted on 7 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Multiphase flow and reactive transport benchmark for radioactive waste disposal

Javier Samper<fnref n="1" /> · Alba Mon<fnref n="1" /> · Etienne Ahusborde<fnref n="2" /> · Han Yu<fnref n="3" /> · Asta Narkuniene<fnref n="4" /> · Milan Hokr<fnref n="5" /> · Luis Montenegro<fnref n="1" /> · Brahim Amaziane<fnref n="2" /> · Mustapha El Ossmani<fnref n="2" />,<fnref n="6" /> · Tianfu Xu<fnref n="3" /> · Yilong Yuan<fnref n="3" /> · Jan Šembera<fnref n="5" /> · Gintautas Poskas<fnref n="4" />

Received: 28 June 2024 / Accepted: 21 September 2024

© The Author(s) 2024

## Abstract

Compacted bentonite is part of the multi-barrier system of radioactive waste repositories. The assessment of the long-term performance of the barrier requires using reactive transport models. Here we present a multiphase flow and reactive transport benchmark for radioactive waste disposal. The numerical model deals with a 1D column of unsaturated bentonite through which water, dry air and CO₂(g) may flow and with the following reactions; aqueous complexation, calcite and gypsum dissolution/precipitation, cation exchange and gas dissolution. INVERSE-FADES-CORE V2, DuMuₓ, TOUGHREACT and iCP were benchmarked with 6 test cases of increasing complexity, starting with conservative tracer transport under variably unsaturated conditions and ending with water flow, gas diffusion, minerals and cation exchange. The solutions of all codes exhibit similar trends. Small discrepancies are found in conservative tracer transport due to differences in hydrodynamic dispersion. Computed CO₂(g) pressures agree when a sufficiently refined grid is used. Small discrepancies in CO₂(g) and pH are found near the no-flow boundary at early times which vanish later. Discrepancies are due differences in the formulations used for gas flow at nearly water-saturated conditions. Computed CO₂(g) pressures show a fluctuation between 10⁻<fnref n="4" /> and 10⁻<fnref n="3" /> years which slows down the in-diffusion of CO₂(g). This fluctuation is associated with chemical reactions involving CO₂. There are discrepancies in solute concentrations due to differences in the Debye–Hückel (DH) formulation. They are overcome when all codes use the same DH formulation. The results of this benchmark will contribute to increase the confidence on multiphase reactive transport models for radioactive waste disposal.

**Keywords** Benchmark · Reactive transport modelling · FEBEX bentonite · Multiphase flow

## Introduction

Reactive Transport Modeling (RTM) is essential for comprehending the interconnected thermal, hydraulic, and chemical (THC) processes crucial for assessing the performance of geological disposal of radioactive waste, as highlighted in previous studies (Bildstein et al. 2019; Claret et al. 2018).

Benchmarks enhance the knowledge of novel challenges by apporting the definition and the resolution of several test cases. Several code benchmark cases for reactive transport models have been proposed. The Groupe de Recherche Mathematical Modeling and Numerical Simulation for Nuclear Waste Management Problems (MoMaS) presented the benchmark of several cases representative of the problems encountered in nuclear waste disposal simulation (Carayrou et al. 2010). According to De Dieuleveult and Erhel (2010), the main challenges of reactive transport simulations deal with solving strongly nonlinear systems with sharp

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

<footnote>Javier Samper
j.samper@udc.es</footnote>

<footnote n="1">¹ Interdisciplinary Center for Chemistry and Biology (CICA), Civil Engineering School & Department, University of A Coruña, Campus Elviña s/n, 15071 A Coruña, Spain</footnote>

<footnote n="2">² E2S UPPA, CNRS, LMAP, Université de Pau et des Pays de l'Adour, Pau, France</footnote>

<footnote n="3">³ College of New Energy and Environment, Jilin University, Changchun 130021, People's Republic of China</footnote>

<footnote n="4">⁴ Nuclear Engineering Laboratory, Lithuanian Energy Institute, Kaunas, Lithuania</footnote>

<footnote n="5">⁵ Faculty of Mechatronics, Informatics and Interdisciplinary Studies, Technical University of Liberec, Liberec, Czech Republic</footnote>

<footnote n="6">⁶ L2M3S-ENSAM, University of Moulay Ismaïl, 50500 Meknes, Morocco</footnote>
fronts and stiff reactions. They reported oscillations due to numerical dispersion which were damped by refining the mesh. In Lagneau and van der Lee (2010), the authors also reported numerical oscillations for a system at thermodynamic equilibrium. Oscillations could be reduced by using a finer grid or by adding a kinetic control. In Alt-Epping et al. (2015), the authors presented a benchmark problem for reactive transport simulators based on a flow-through experiment carried out on saturated bentonite. The results for different codes showed excellent agreement. It was required to modify their thermodynamic databases (consistent activity coefficients and thermodynamic data) to achieve this good agreement. The FluidFlower study considers the injection of CO₂ into brine-saturated formations close to atmospheric conditions (Wapperom et al. 2024; Tian et al. 2024). They establish that the resolution and comparison of unstructured and Cartesian gridding reveals a high sensitivity of the overall distribution of CO₂ throughout the domain to gridding choices. They reported that the grid resolution plays a very important role in FluidFlower modeling. A refined grid is needed for a fully numerically converged solution. In Ahusborde et al. (2024) and De Hoop et al. (2024), the authors presented a benchmark case for CO₂ geological storage with focus on the coupling of two-phase flow and geochemistry. They showed that the results from all methods were comparable, even though there remained some discrepancies.

Previous subsurface environmental simulation benchmarks have stressed the need for benchmarking in the context of multiphase multicomponent reactive transport for radioactive waste disposal (Bildstein et al. 2021).

Here we report a multiphase flow and reactive transport benchmark for radioactive waste disposal. It is based on the conditions of the FEBEX (Full-scale Engineering Barrier Experiment) in situ test (Samper et al. 2018a, 2008; Zheng et al. 2011). The benchmark was performed with 4 multiphase flow and reactive transport codes (INVERSEFADES-CORE V2, DuMuX, TOUGHREACT and iCP) and involves 6 test cases of increasing complexity, including: conservative tracer migration, hydration under variably unsaturated conditions, gas diffusion, calcite and gypsum dissolution/precipitation, and cation exchange reactions.

In “Description of the benchmark: FEBEX in situ test” section, the paper starts with the description of the benchmark, including the mathematical formulation, the description of the FEBEX in situ test, the simplifications adopted for the benchmark and the main features of the model such as physical properties, space and time discretization, thermodynamic database, geochemical properties and initial and boundary conditions. Then, in “Description of the codes” section, the features of the computer codes utilized in the benchmark are presented. Afterwards, in “Comparison and discussion of the results” section, a selection of the benchmark results for test cases TC1 to TC6 are presented, followed by a comparative

analysis and discussion of results obtained by different teams. The paper ends with the main conclusions and perspectives for future work in “Conclusion and perspectives” section.

## Description of the benchmark: FEBEX in situ test

In this section, a comprehensive description of the benchmark is presented. The multiphase mathematical formulation is outlined. Following this, the various proposed cases of increasing complexity are detailed to provide a clear understanding of the scenarios under consideration.

## Multiphase mathematical formulation

This section presents the mathematical formulation of the coupled hydrological-chemical model. This formulation has been extracted from Zheng et al. (2011). Water mass balance is given by Navarro and Alonso (2000):

$$ \frac{\partial (\Phi \rho_l X_l^w S_l + \Phi \rho_g X_g^v S_g)}{\partial t} + \nabla \cdot (\rho_l X_l^w \mathbf{q}_l + \rho_g X_g^v \mathbf{q}_g + \mathbf{j}_v) + r_w = 0, \quad (1) $$

where $\Phi$ (-) is the porosity, $S_\alpha$ (-) is the saturation degree of phase $\alpha$ ($\alpha \in \{l, g\}$). $\rho_l$ and $\rho_g$ (kg/m³) are the bulk densities of the liquid and gaseous phases, respectively, $X_l^w$ (-) is the mass fraction of water in the liquid phase, $X_l^v$ (-) is the mass fraction of the vapor in the gas phase, $r_w^v$ (kg/m³/s) is the sink/source term of liquid water, $\mathbf{q}_\alpha$ (m/s) is the vector of volumetric flux of phase $\alpha$ which is given by Eq. 4, $\mathbf{j}_v$ (kg/m²/s), is the dispersive mass flux of vapor with respect to the mean gas velocity which is given by Eq. 5. The dispersive mass fluxes of air and water with respect to the liquid phase are neglected.

The gas mass balance equation is given by:

$$ \frac{\partial \Phi \rho_g S_g}{\partial t} + \nabla \cdot (\rho_g \mathbf{q}_g) + r_g = 0, \quad (2) $$

where $r_g$ (kg/m³/s) is the sink/source term of gas. The dispersive mass flux of air and the gas with respect to the gaseous phase are neglected.

The air mass balance equation is given by:

$$ \frac{\partial (\Phi \rho_g X_g^a S_g + \Phi \rho_l H S_l)}{\partial t} + \nabla \cdot (\rho_g X_g^a \mathbf{q}_g + \rho_l H \mathbf{q}_l) = 0, \quad (3) $$

where $X_g^a$ (-) is the mass fraction of air in the gaseous phase, $H$ (-) is the Henry's constant.
Our formulation assumes isothermal conditions and that all phases are at local thermal equilibrium and therefore they are all at the same temperature.

The volumetric flux $q_α$ (m/s) for each phase $α$ ($α \in \{l, g\}$), is given by:

$$q_α = -\frac{K_{iα}k_{rα}}{\mu_α}(\nabla p_α + \rho_α g \nabla z), \quad (4)$$

where $p_α$ (Pa) is the pressure of phase $α$, $K_{iα}$ (m²) is the intrinsic permeability tensor of the phase $α$, $k_{rα}$ (-) is the relative permeability of the phase $α$, $μ_α$ (Pa.s) is the viscosity of the phase $α$, $z$ (m) is the elevation and $g$ (m²/s) is the gravitational acceleration.

The dispersive mass flux of vapor $j_v$ (kg/m²/s) is calculated by the Fick's law:

$$j_v = -\rho_g D_v \nabla X_g^v, \quad (5)$$

where $D^v$ (m²/s) is the hydrodynamic dispersion tensor for vapor which includes the effects of mechanical dispersion ($D_{disp}$) and molecular diffusion ($D^e$). The effective molecular diffusion coefficient for the vapor was calculated from Pollock (1986):

$$D_v^e = \frac{5.9 \ 10^{-6} \tau^v (T + 273.15)^{2.3}}{p_g}, \quad (6)$$

where $τ^v$ (-) is the vapor tortuosity factor and $T$ is the temperature expressed in Celsius.

It should be noticed that in the model the gases dissolved in water (air and other gases) are transported by advection and dispersion in the fluid phase together with liquid water.

Solute transport processes include advection, molecular diffusion, and mechanical dispersion. Each of them produces a solute flux per unit surface and unit time. There are as many transport equations as primary chemical species or aqueous components in the system. The primary species are the building blocks of chemical systems of interest, upon which concentrations of secondary species are written through laws of mass action for reactions at thermodynamic equilibrium.

The mass balance equation for the j-th primary species is given by Zheng and Samper (2008):

$$m_l^w \frac{\partial C_j}{\partial t} + \frac{\partial(m_l^w P_j)}{\partial t} + \frac{\partial(m_l^w W_j)}{\partial t} + \frac{\partial(m_l^w G_j)}{\partial t} = L^*(G_j) + r_j(C_j^0 - C_j), \quad (7)$$

where $C_j$ (mol/L) is the total dissolved concentration of the j-th primary species, $m_l^w$ (kg/m³) is the mass of liquid water per unit volume of medium, which is equal to $ρ^l X_l^w θ_l$, where $θ_l = S_l Φ$ is the volumetric watercontent (m³/m³), $P_j$, $W_j$ and $G_j$ are the total precipitated minerals, the total exchanged and total dissolved gas concentrations (mol/L),

respectively, of the j-th primary species, $r_j$ (kg/m³/s) is the sink term, $C_j^0$ (mol/L) is the dissolved concentration of j-th species in the sink term $r_j$, $N_c$ is the number of primary species. $L^*(·)$ defines the transport operator as follows:

$$L^*(\cdot) = \nabla \cdot [m_l^w D^j \cdot \nabla(\cdot)] - m_l^w \mathbf{q}^l \cdot \nabla(\cdot) + (r_e - r_c)(\cdot), \quad (8)$$

where $D^j$ (m²/s) defines the dispersion coefficient including the hydrodynamic or mechanical dispersion and the porewater molecular diffusion while $r_e$ and $r_c$ are the condensation and evaporation rates (kg/m³/s), respectively.

The reactive transport of the f-th gas species was implemented as an additional mass balance equation of the f-th gas species in the gaseous phase. The gas transport equation can be written in compact form as:

$$\rho_g \Phi S_g \frac{\partial C_f}{\partial t} = L_f^*(C_f) + r_f^i(C_f^o - C_f) + R_f, \quad (9)$$

where $C_f$ (mol/kg) is the concentration of the f-th gas species in the gaseous phase, $r_f^i$ (kg/m³/s) is the gas mass flux source entering, $C_f^o$ (mol/kg) is the external concentration of the f-th gas species in the entering gas flux and $R_f$ (mol/m³/s) is the gas chemical reactions term. Finally $L_f^*(·)$ is the following transport operator:

$$L_f^*(\cdot) = \nabla \cdot [\rho_g \Phi S_g D^f \cdot \nabla(\cdot)] - \rho^g \mathbf{q}^g \cdot \nabla(\cdot), \quad (10)$$

where $D^f$ (m²/s) defines the dispersion coefficient including the mechanical dispersion and the molecular diffusion.

The chemical conceptual model for compacted bentonite accounts for the following reactions: aqueous complexation, acid/base, redox, cation exchange, mineral dissolution/precipitation (at equilibrium or in kinetic), and gas dissolution/exsolution. The chemical system is defined in terms of the concentrations of the primary species. The concentrations of the secondary species are computed from the concentrations of the primary species through appropriate mass action laws (Samper et al. 2009). The concentrations of the precipitated, exchanged and adsorbed species are computed using similar equations. A detailed description of the calculations of the chemical reactions can be found in Samper et al. (2009) and Zheng et al. (2011). The Gaines-Thomas convention is used for cation exchange. The total dissolved concentration $C_j$ (mol/L), can be written in an explicit form as a function of the $N_c$ primary species by applying the mass action law:
$$C_j = c_j + \sum_{k=1}^{N_x} v_{kj} x_k \\ = c_j + \sum_{k=1}^{N_x} v_{kj} \left( K_k^{-1} \gamma_k^{-1} \prod_{i=1}^{N_c} c_i^{\gamma_{ki}} \gamma_i^{\gamma_{ki}} \right) \quad (11)$$

where $c_j$ (mol/L) is the dissolved concentration of the primary species j, $N_x$ (-) is the number of secondary species, $x_j$ (mol/L) is the concentration of the secondary species, $K_k$ (-) is the equilibrium constant of the k-th secondary species reaction, $\gamma$ (-) is the thermodynamic activity coefficient and $v_{kj}$ (-) is the stoichiometric coefficient of the j-th primary species on the k-th species.

The total concentration of the precipitated minerals of the j-th primary species, $P_j$ (mol/L) can be written as:

$$P_j = \sum_{m=1}^{N_p} v_{mj}^p p_m, \quad (12)$$

where $p_m$ (mol/L) is the concentration of the m-th mineral phase and $v_{mj}^p$ (-) is the stoichiometric coefficient of the j-th primary species on the m-th mineral.

Under equilibrium conditions, dissolution-precipitation reactions can be described by the mass action law which states that:

$$X_m \lambda_m K_m = \prod_{i=1}^{N_c} c_i^{\gamma_{mi}} \gamma_i^{\gamma_{mi}}, \quad (13)$$

where $X_m$ (-) is the molar fraction of the m-th solid phase; $\lambda_m$ (-) is the thermodynamic activity coefficient ($X_m$ and $\lambda_m$ are taken equal to 1 for pure phases); $c_i$ (mol/L) and $\gamma_i$ (-) are the concentration and activity coefficient of the i-th specie; $v_{mi}$ (-) is the stoichiometric coefficient in the dissolution reaction of the m-th solid phase; and $K_m$ (-) is the corresponding equilibrium constant.

The concentration of the i-th exchanged cation (mol/L) can be obtained from the i-th equivalent fraction $\beta_i$, according to:

$$w_i = \beta_i CEC \rho_s \frac{(1 - \Phi)}{z_i 100 \Phi S_l} \\ = (K_{ij}^*)^{-\nu_i} c_i \gamma_i \left( \frac{\beta_j}{c_j \gamma_j} \right)^{\frac{\nu_i}{\gamma_j}} CEC \rho_s \frac{(1 - \Phi)}{z_i 100 \Phi S_l}, \quad (14)$$

where CEC (meq/100g) is the total cation exchange capacity, $z_i$ (-) is the cation charge and $K_{ij}^*$ (-) is the exchange coefficient or selectivity.

The total concentration of the dissolved gases concentrations of the j-th primary species $G_j$ (mol/kg), can be written as:

$$G_j = \sum_{f=1}^{N_g} v_{ff}^g g_f, \quad (15)$$

where $g_f$ (mol/kg) is the concentration of the f-th gas and $v_{ff}^g$ (-) is the stoichiometric coefficient of the j-th primary species on the f-th gas. For reactions involving aqueous and gas phases, the mass action law states that:

$$p_f \Gamma_f K_f = \prod_{i=1}^{N_c} c_i^{\nu_{fi}} \gamma_i^{\nu_{fi}}, \quad (16)$$

where $p_f$ (Pa) is the partial pressure of the f-th species in the gas phase; $\Gamma_f$ (-) is the activity coefficient and $K_f$ (-) is the equilibrium constant of the reaction, $c_i$ (mol/L) and $\gamma_i$ (-) are the concentration and activity coefficient of the i-th dissolved primary species, $v_{fi}$ (-) is the stoichiometric coefficient of the f-th gas on the i-th specie. The gaseous phase is assumed as an ideal mixture and the fugacity constant, $\Gamma_f$ is equal to 1.

The equilibrium constants for aqueous complexes and minerals depend on temperature. They are calculated with the following expression, which is valid for temperatures between 0 and 300 °C:

$$\log K(T) = \frac{b_1}{T^2} + \frac{b_2}{T} + b_3 \log T + b_4 + b_5 T, \quad (17)$$

where $b_1$ to $b_5$ are coefficients taken in the thermodynamic database of ThermoChimie v11.a (Giffaut et al. 2014), used for aqueous complexes and minerals.

For the activity coefficients of the aqueous species the extended Debye–Hückel formula can be used:

$$\log \gamma_i = \frac{-A z_i^2 (I)^{1/2}}{1 + B a_i (I)^{1/2}} + bI, \quad (18)$$

where $I$ (-) is the ionic strength of the solution; $z_i$ (-) and $a_i$ (-) are the electric charge and the ionic radius in solution of the i-th species, respectively, A and B are constants which depends on the temperature and dielectric constant of water, and b is a constant determined from experimental data (A, B and b were taken from tabulated values from Helgeson and Kirkham (1974)). The value of the ionic strength is calculated as:

$$I = \frac{1}{2} \sum_{i=1}^{N_T} c_i z_i^2, \quad (19)$$

where $N_T$ is the number of ions present in the solution.

The activity of the water can be calculated according to the approximation of Garrels and Christ (1965):

$$a_{H_2O} = 1 - 0.018 \sum_{i=2}^{N_T} c_i. \quad (20)$$
## Description of the FEBEX in situ test reference case

The FEBEX in situ test was performed in a gallery excavated in granite in the underground research laboratory of Grimsel operated by NAGRA in Switzerland. The test included the heating system, the clay barrier and the instrumentation, monitoring and control system (Fig. 1).

The drift was 70.4 m long and 2.27 m in diameter. The test zone was located in the last 17.4 m of the drift where heaters, bentonite and instrumentation were installed. The main elements of the heating system were two heaters, separated horizontally by 1 m, which simulated full-sized canisters. The heaters were placed inside a cylindrical steel liner having a diameter of 0.93 m, which had been installed concentrically with the drift. Each heater was made of carbon steel, measured 4.54 m in length and 0.90 m in diameter, had a wall thickness of 0.10 m and weighed 11 tons. The heaters were designed to maintain a maximum temperature of 100 °C at the liner/bentonite interface. The bentonite barrier was made of blocks of highly compacted bentonite. The test began in February 1997. The first operation period lasted from 1997 to 2002 when heater 1 was switched off and the surrounding area was dismantled. The second operation period started after the emplacement of a shotcrete plug and ended in June 2015 when the entire bentonite barrier was fully dismantled.

The bentonite barrier is hydrated from the external cylindrical surface and water flows towards the internal heater/bentonite interface. The bentonite barrier is initially unsaturated, and it progressively hydrates from the surrounding rock towards the canister. The water content of the bentonite increases near the hydration boundary. Bentonite hydration leads to bentonite swelling. Therefore, the porosity of the bentonite increases during bentonite hydration. Mechanical and swelling processes play an important role in the hydrodynamic and geochemical evolution of the engineered barrier system during its early heating and hydration stage, when bentonite buffer is subjected simultaneously to heating and hydration.

The temperature is fixed at 100 °C on the heater/bentonite interface. Water evaporates near the heater. Vapor diffuses away from the heater and condenses in cooler places. Vapor condensation retards the hydration of the bentonite buffer and affects the concentration of the dissolved species (Villar et al. 2012; Zheng et al. 2011).

The conditions of the FEBEX in situ test call for a coupled non-isothermal multiphase flow conceptual model with: 1) Advection of water in the liquid phase; 2) Advection and diffusion of vapor and other gases in the gaseous phase; 3) Advection and diffusion of air in the liquid and gaseous phases; 4) Convection of heat in the liquid and gaseous phase and; 5) Heat conduction.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** General layout of the FEBEX in situ test indicating the instrumented and sampling sections used in this work. The x coordinates of the sections are referred to the concrete plug on the left</figcaption></figure>
## Simplifications and benchmark test cases

The conceptual model for the multiphase simulations needs to be simplified due to the complexity of the real system. It was agreed to consider a simplified system in the reference case and to increase the complexity of the modelled processes. The main simplification pertains to the temperature, which is assumed to be constant across all test cases, thus modeling isothermal processes. The proposed test cases (TC) are detailed in Table 1.

## Model description

### Physical properties

Table 2 shows the hydrodynamic and transport properties used in the numerical model of the FEBEX in situ test.

### Time and space discretization

The numerical models were performed with two one-dimensional meshes. The first one is uniform (see bottom of Fig. 2). The grid size of the elements, denoted as $Δx$, is equal to 0.01 m throughout the entire domain. The second mesh,

<table><caption>**Table 1** Description of benchmark test cases TC1 to TC6 for the FEBEX in situ test</caption><thead><tr><th></th><th>TC1</th><th>TC2</th><th>TC2B</th><th>TC3</th><th>TC3B</th><th>TC4</th><th>TC5</th><th>TC6</th></tr></thead><tbody><tr><td>Hydration</td><td>Yes</td><td>No</td><td>No</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Aqueous Species</td><td>Tracer</td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup></td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup><br>Ca<sup>2+</sup><br>Na<sup>+</sup><br>HCO<sub>3</sub><sup>-</sup><br>Cl<sup>-</sup><br>O<sub>2</sub>(aq)<br>OH<sup>-</sup><br>CaCl<sup>+</sup><br>CaCO<sub>3</sub><br>Ca(OH)<sup>+</sup><br>NaCl<br>CaCl<sub>2</sub><br>Ca(HCO<sub>3</sub>)<sup>+</sup><br>Na(HCO<sub>3</sub>)<br>Na(OH)<br>Na(CO<sub>3</sub>)<sup>-</sup></td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup><br>Ca<sup>2+</sup><br>Na<sup>+</sup><br>HCO<sub>3</sub><sup>-</sup><br>Cl<sup>-</sup><br>O<sub>2</sub>(aq)<br>OH<sup>-</sup><br>CaCl<sup>+</sup><br>CaCO<sub>3</sub><br>Ca(OH)<sup>+</sup><br>NaCl<br>CaCl<sub>2</sub><br>Ca(HCO<sub>3</sub>)<sup>+</sup><br>Na(HCO<sub>3</sub>)<br>Na(OH)<br>Na(CO<sub>3</sub>)<sup>-</sup><br>H<sub>2</sub>(aq)</td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup><br>Ca<sup>2+</sup><br>Na<sup>+</sup><br>SO<sub>4</sub><sup>2-</sup></td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup><br>Ca<sup>2+</sup><br>Na<sup>+</sup><br>SO<sub>4</sub><sup>2-</sup></td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup><br>Ca<sup>2+</sup><br>Na<sup>+</sup><br>SO<sub>4</sub><sup>2-</sup></td><td>H<sup>+</sup><br>CO<sub>2</sub>(aq)<br>CO<sub>3</sub><sup>2-</sup><br>Ca<sup>2+</sup><br>K<sup>+</sup><br>Mg<sup>2+</sup></td></tr><tr><td>Minerals</td><td>No</td><td>No</td><td>No</td><td>Calcite</td><td>Calcite</td><td>Calcite</td><td>Calcite<br>Gypsum</td><td>Calcite<br>Gypsum</td></tr><tr><td>Gases</td><td>Air</td><td>Air, CO<sub>2(g)</sub></td><td>Air, CO<sub>2(g)</sub></td><td>Air, CO<sub>2(g)</sub></td><td>Air, CO<sub>2(g)</sub></td><td>Air, CO<sub>2(g)</sub></td><td>Air, CO<sub>2(g)</sub></td><td>Air, CO<sub>2(g)</sub></td></tr><tr><td>Cation exchange</td><td>No</td><td>No</td><td>No</td><td>No</td><td>No</td><td>No</td><td>No</td><td>Na<sup>+</sup>, K<sup>+</sup></td></tr><tr><td>Boundary<br>r = 1.14 m</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>P<sub>l</sub> = 10<sup>5</sup> Pa<br>S<sub>l</sub> = 60%<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>S<sub>l</sub> = 60%<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>S<sub>l</sub> = 60%<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>S<sub>l</sub> = 60%<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>P<sub>l</sub> = 10<sup>5</sup> Pa<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>P<sub>l</sub> = 10<sup>5</sup> Pa<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>P<sub>l</sub> = 10<sup>5</sup> Pa<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td><td>P<sub>g</sub> = 10<sup>5</sup> Pa<br>P<sub>l</sub> = 10<sup>5</sup> Pa<br>X<sub>g</sub><sup>CO<sub>2</sub></sup> = 0.1</td></tr></tbody></table>
<table><caption>**Table 2** Hydrodynamic and transport parameters (ENRESA 2006; Mon 2017; Samper et al. 2018a; Zheng and Samper 2008; Zheng et al. 2011, 2010)</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Porosity, Φ (-)</td><td>Φ = 0.4</td></tr><tr><td>Intrinsic permeability for liquid flow, K<sub>il</sub> (m<sup>2</sup>)</td><td>K<sub>il</sub> = k<sub>o</sub> (Φ<sup>3</sup> / (1-Φ)<sup>2</sup>) (1-Φ<sub>o</sub>)<sup>2</sup> / Φ<sub>o</sub><sup>3</sup></td></tr><tr><td>Intrinsic permeability for gas flow K<sub>ig</sub> (m<sup>2</sup>)</td><td>k<sub>o</sub> = 3.75 × 10<sup>-21</sup> m<sup>2</sup></td></tr><tr><td>Relative permeability to liquid, k<sub>rl</sub> (-)</td><td>K<sub>ig</sub> = 5 × 10<sup>-10</sup></td></tr><tr><td>Relative permeability to gas, k<sub>rg</sub> (-)</td><td>k<sub>rl</sub> = S<sub>l</sub><sup>3</sup></td></tr><tr><td>Van Genuchten retention curve, Ψ (kPa)</td><td>k<sub>rg</sub> = (1 - S<sub>l</sub>)<sup>3</sup><br>S<sub>l</sub> = S<sub>r</sub> + (S<sub>r</sub> - S<sub>g</sub>) / (1 + αΨ<sup>1/3</sup>)<sup>b</sup></td></tr><tr><td>Liquid viscosity μ<sub>l</sub> (kg/(m.s))(T in Celsius)</td><td>S<sub>r</sub> = 0.05, S<sub>s</sub> = 0.95</td></tr><tr><td>Gas viscosity μ<sub>g</sub> (kg/(m.s))</td><td>α = 5 × 10<sup>-5</sup> kPa<sup>-1</sup>, b = 0.21</td></tr><tr><td>Vapor tortuosity factor τ<sup>v</sup> (-)</td><td>μ<sub>l</sub> = 661.2 × 10<sup>-3</sup> (T + 44)<sup>-1.562</sup></td></tr><tr><td>Liquid density ρ<sub>l</sub> (kg/m<sup>3</sup>) (T in Celsius)</td><td>μ<sub>g</sub> = 1.76 × 10<sup>-5</sup></td></tr><tr><td>Gas density ρ<sub>g</sub> (kg/m<sup>3</sup>) (T in Celsius)</td><td>τ<sup>v</sup> = 0.10</td></tr><tr><td>Molecular diffusion in water D<sub>o</sub> (m<sup>2</sup>/s) at 25 °C</td><td>ρ<sup>l</sup> = 998.2 e<sup>(5×10<sup>-7</sup>(P<sub>l</sub>-100)-2.1×10<sup>-4</sup>(T-12))</sup></td></tr><tr><td>Air Henry's constant H (-) (T in Celsius)</td><td>ρ<sub>g</sub> = ρ<sub>v</sub> + ρ<sub>a</sub><br>ρ<sub>v</sub> = (e<sup>(0.06374T-0.1634×10<sup>-3</sup>τ<sup>2</sup>)</sup> / 194.4) e<sup>-(2.16677Ψ)/(ρ<sub>l</sub>(T+273.15))</sup></td></tr><tr><td></td><td>ρ<sub>a</sub> = (3.499 P<sub>g</sub> / (T+273.15)) - 1.615ρ<sub>v</sub></td></tr><tr><td></td><td>D<sub>o</sub> = 2 × 10<sup>-11</sup></td></tr><tr><td></td><td>H = ρ<sub>l</sub>(T + 273.15)7.1281183 × 10<sup>-11</sup> e<sup>(1997.32 / (T+273.15))</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Schematic of the meshes</figcaption></figure>

composed of 104 nodes is represented in the upper part of the Fig. 2 and it is non uniform. Δx is equal to 0.01 m in the interval 0.45 m < x < 1 m, Δx is equal to 5 × 10⁻³ m in the interval 1 m < x < 1.1 m; Δx is equal to 2.5 × 10⁻³ m in the interval 1.1 m < x < 1.13 m; and Δx is equal to 6.25 × 10⁻⁴ m at the right boundary (1.13 m < x < 1.14 m). The model domain included the bentonite barrier, which extended from r = 0.45 m to r = 1.14 m. The simulation time horizon covered the entire duration of the in situ test from February 1997 to 2015 (18 years).

## Thermodynamic database

The benchmark cases consider the thermodynamic database ThermoChimie version 11a (Giffaut et al. 2014), available in different formats at https://www.thermochimie-tdb.com/.

## Geochemical properties

The primary species considered in the geochemical model of the benchmark cases are the following: H₂O, H⁺, Na⁺, K⁺, Ca²⁺, Mg²⁺, Cl⁻, SO₄²⁻ and CO₃²⁻. Each test case includes different dissolves species and minerals which are shown in Table 1. Table 3 shows the chemical
<table><caption>**Table 3** Chemical reactions and equilibrium constants for aqueous complexes, minerals, gases and selectivity constants for cation exchange reactions at 25 °C taken from the thermodynamic database ThermoChimie v11.a (Giffaut et al. 2014) to use in the model of the FEBEX in situ test benchmark case</caption><thead><tr><th>Aqueous complexes</th><th>Log K</th></tr></thead><tbody><tr><td>HCO<sub>3</sub><sup>-</sup> ⇌ H<sup>+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 10.3300</td></tr><tr><td>CO<sub>2</sub> (aq) + H<sub>2</sub>O ⇌ 2H<sup>+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 16.6800</td></tr><tr><td>OH<sup>-</sup> + H<sup>+</sup> ⇌ H<sub>2</sub>O</td><td>14.000</td></tr><tr><td>H<sub>2</sub>(aq) + 0.5 O<sub>2</sub>(aq) ⇌ H<sub>2</sub>O</td><td>+ 46.07</td></tr><tr><td>CaCl<sup>+</sup> ⇌ Ca<sup>2+</sup> + Cl<sup>-</sup></td><td>+ 0.2900</td></tr><tr><td>CaCO<sub>3</sub>(aq) ⇌ Ca<sup>2+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 3.2200</td></tr><tr><td>Ca(OH)<sup>+</sup> + H<sup>+</sup> ⇌ Ca<sup>2+</sup> + H<sub>2</sub>O</td><td>+ 12.78</td></tr><tr><td>NaCl<sub>(aq)</sub> ⇌ Na<sup>+</sup> + Cl<sup>-</sup></td><td>+ 0.5000</td></tr><tr><td>CaCl<sub>2</sub>(aq) ⇌ Ca<sup>2+</sup> + 2Cl<sup>-</sup></td><td>+ 0.6400</td></tr><tr><td>Ca(HCO<sub>3</sub>)<sup>+</sup> ⇌ Ca<sup>2+</sup> + H<sup>+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 11.43</td></tr><tr><td>NaHCO<sub>3</sub>(aq) ⇌ Na<sup>+</sup> + H<sup>+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 10.0800</td></tr><tr><td>Na(CO<sub>3</sub>)<sup>-</sup> ⇌ Na<sup>+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 1.2700</td></tr><tr><td>Na(OH) + H<sup>+</sup> ⇌ Na<sup>+</sup> + H<sub>2</sub>O</td><td>+ 14.75</td></tr><tr><td colspan="2"><strong>Minerals</strong></td></tr><tr><td>Calcite ⇌ Ca<sup>2+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 8.4800</td></tr><tr><td>Gypsum ⇌ Ca<sup>2+</sup> + SO<sub>4</sub><sup>2-</sup> + 2H<sub>2</sub>O</td><td>- 4.6100</td></tr><tr><td colspan="2"><strong>Gases</strong></td></tr><tr><td>CO<sub>2</sub>(g) + H<sub>2</sub>O ⇌ 2H<sup>+</sup> + CO<sub>3</sub><sup>2-</sup></td><td>- 18.1500</td></tr><tr><td colspan="2"><strong>Cation exchange</strong></td></tr><tr><td>Na<sup>+</sup> + X-K ⇌ K<sup>+</sup> + X-Na</td><td>0.138</td></tr><tr><td>Na<sup>+</sup> + 0.5 X<sub>2</sub> - Ca ⇌ 0.5 Ca<sup>2+</sup> + X-Na</td><td>0.294</td></tr><tr><td>Na<sup>+</sup> + 0.5 X<sub>2</sub> - Mg ⇌ 0.5 Mg<sup>2+</sup> + X-Na</td><td>0.288</td></tr></tbody></table>

reactions and the equilibrium constants at 25 °C for the secondary aqueous species, minerals and gases used in the numerical model of the FEBEX in situ test benchmark case. Table 3 also depicts the selectivity constants for the cation exchange reactions in the bentonite. The total cation exchange capacity, CEC, is equal to 102.75 meq/100g.

## Initial and boundary conditions

Bentonite had an initial porosity of 0.40, a volumetric water content of 24%, which corresponds to a gravimetric water content of around 14.4%, a liquid saturation degree of 60% and a suction of 1.17 10⁵ kPa. The gas pressure was set to 100 kPa. The initial temperature was uniform, set at 25 °C, and assumed to remain constant. The initial compositions of the bentonite and the granite pore water, the initial mineral volume fractions

and the initial concentrations of exchanged ions are listed in Tables 4, 5 and 6 which can be found in the “Appendix”.

The temperature and the liquid and gas pressure at the outer boundary (r = 1.14 m) were equal to 25 °C and 100 kPa, respectively (see Fig. 2 and Table 1).

## Description of the codes

Five teams participated in the benchmark with various codes. We are listing the different teams along with their codes, whose main features are described below.

* University of A Coruña (UDC) with INVERSE-FADES-CORE V2;

* University of Pau and the Adour Region (UPPA), with DuMu<sup>X</sup>;

* Jilin University (JU) with TOUGHREACT;

* Technical University of Liberec (TUL) with iCP;

* Lithuanian Energy Institute (LEI) with iCP.

## INVERSE-FADES-CORE V2: UDC

INVERSE-FADES-CORE V2 is a finite element code for modelling non-isothermal multiphase flow, heat transport and multicomponent reactive solute transport under both chemical equilibrium and kinetics conditions (Zheng and Samper 2015; Zheng et al. 2010). The code takes into account the mass balance of water, air, solid and enthalpy; the transport of solids and mechanical equilibrium. The solute transport equation accounts for advection, molecular diffusion, and mechanical dispersion. The reactive gas transport includes additional mass balance equations for the reactive gaseous species in the gaseous phase. The state variables for nonisothermal flow include liquid and gas pressures and temperature. INVERSE-FADES-CORE V2 solves both forward and inverse multiphase flow and multicomponent reactive transport problems in 1-, 2- and 3-D axisymmetric porous and fractured media. The code is the result of integrating the capabilities of a thermo-hydro-mechanical code FADES (Navarro and Alonso 2000), the reactive transport code CORE<sup>2D</sup> (Xu et al. 1999), and the inverse methodology of INVERSE-CORE (Dai and Samper 2004). The code accounts for aqueous complexation, acid-base, redox, mineral dissolution/precipitation, gas dissolution/exsolution, cation exchange, and surface complexation reactions. INVERSE-FADES-CORE V2 solves first the multiphase flow and energy transport equations and then it solves for coupled transport processes and geochemical reactions by adopting the sequential iteration approach. Transport and chemical equilibrium equations are solved separately in an iterative sequential manner. At any given iteration, the chemical
sink/source term is assumed known (or taken equal to the value at the previous iteration) for the purpose of solving the transport equations. This renders the differential equations linear which allows the use of fast numerical methods for solving the linear system of equations. After solving the transport equations, the set of chemical equations is solved nodewise. The nonlinear chemical equations are solved with an iterative Newton-Raphson method. After convergence, the chemical source term is updated. The whole iterative process (transport and chemistry) is repeated until overall convergence is attained. The finite element method is used for spatial discretization, while implicit, explicit and Crank-Nicholson schemes are used for time discretization. INVERSE-FADES-CORE V2 has been used to model laboratory tests and in situ tests (Mon et al. 2023; Samper et al. 2018a, b, 2020a, b, 2008; Zheng and Samper 2008; Zheng et al. 2011, 2010, 2008).

## DuMu<sup>X</sup>: UPPA

DuMu<sup>X</sup> (DUNE for Multi-Phase, Component, Scale, Physics, ... flow and transport in porous media) (Koch et al. 2021) is a free and open-source simulator designed for modeling flow and transport processes in porous media. It is built on the Distributed and Unified Numerics Environment (DUNE) (Bastian et al. 2021), which is an object-oriented software framework written in C++. DUNE handles essential computational tasks such as input/output operations, memory management, grid generation, and parallel computing.

Over the years, UPPA has contributed to DuMu<sup>X</sup> by implementing various numerical schemes for reactive transport modeling. In particular, works such as Ahusborde et al. (2015, 2018) and Ahusborde and El Ossmani (2017) introduced a sequential approach that decomposes the global problem into two sub-problems. The first sub-problem addresses two-phase compositional flow, where only species present in both phases are treated implicitly, and phase exchanges are fully resolved in this step, while other species are treated explicitly. The second sub-problem involves solving a reactive transport problem, using the flow properties (e.g., Darcy velocity, phase saturations, temperature, and phase densities) computed in the first step. In Ahusborde et al. (2015) and Ahusborde and El Ossmani (2017), a sequential iterative approach (SIA) was utilized for the reactive transport sub-problem. However, to mitigate potential time-splitting errors associated with SIA, a global implicit approach was adopted in Ahusborde et al. (2018). Subsequent developments led to the creation of a fully implicit, fully coupled method using a direct substitution approach, as detailed in Ahusborde et al. (2019, 2021). More recently, these methodologies have been extended to study coupled Thermo-Hydro-Chemical processes for both fully

implicit and sequential approaches (Ahusborde et al. 2023). This contribution focuses solely on results derived from the fully implicit approach.

The spatial discretization employs a cell-centered finite volume (FV) scheme. Convective terms are approximated using a fully upwind scheme, while diffusive terms are calculated using a two-point flux approximation (TPFA). This approach is suitable for the orthogonal meshes used in the benchmark geometries. For more complex meshes, multi-point flux approximations (MPFA) could be implemented. The nonlinear system is resolved using a Newton-Raphson algorithm, with the Jacobian matrix approximated by numerical differentiation. The resulting linear systems are solved using the BiConjugate Gradient STABilized (BiCGSTAB) method, preconditioned with an Algebraic Multigrid (AMG) solver. An adaptive time-stepping strategy, based on the iteration count of the Newton method from the previous time step, is employed. A comprehensive description of our methodology is provided in Ahusborde et al. (2021, 2023). This methodology has been validated through multiple test cases, including high-performance computing applications and studies on the geological storage of CO<sub>2</sub> in deep saline aquifers.

## TOUGHREACT: JU

TOUGHREACT is a code for simulating non-isothermal multiphase fluid flow, heat transport and multicomponent reactive solute transport (Xu and Pruess 1998; Xu et al. 2006, 2011). It has been developed by introducing reactive transport into the existing framework of a non-isothermal multi-component fluid and heat flow simulator TOUGH2 (Pruess et al. 1999). A series of subsurface thermo-physical-geochemical processes are considered under various thermo-hydrological and chemical conditions of pressure, temperature, water saturation, and ionic strength. TOUGHREACT can be used in one-, two-, or three dimensional porous and fractured media with physical and chemical heterogeneity. The first version of the TOUGHREACT code was released to the public through the US Department of Energy's Energy Science and Technology Software Center (ESTSC) in August 2004. The mass and energy balance equations in TOUGHREACT are solved by Newton-Raphson iterations and proceeds as in TOUGH2 (Pruess et al. 1999). Space discretization involves an unstructured finite volume scheme (Integral Finite Differences, IFD). Given the chemical transport equations have the same structure as the fluid and heat flow equations, the transport equations can be solved by the same numerical method. TOUGHREACT uses a sequential iteration approach. After solution of the flow equations, the fluid velocities and phase saturations are used for chemical transport modeling. The chemical transport is solved
on a basis of component-by-component. The resulting concentrations obtained from solving transport equations are substituted into the chemical reaction model. The system of mixed equilibrium-kinetic chemical reaction equations is solved on a grid block by grid block basis by Newton-Raphson iteration. Optionally, the chemical transport and reactions are solved iteratively until convergence. The TOUGHREACT code has been widely applied for studies in CO₂ geological sequestration (Gherardi et al. 2007; Hu et al. 2023b, 2019; Zhu et al. 1999), geothermal resources development (Chen et al. 2020; Spycher et al. 2014; Wanner et al. 2014), nuclear waste isolation (Xu et al. 2008), and increasingly for petroleum applications (Hu et al. 2023a).

## iCP: TUL and LEI

iCP v2.1 is an interface developed by Amphos21 (Spain) to couple two standalone simulation programs: the general-purpose finite element platform COMSOL Multiphysics (COMSOL 2022) and the geochemical code PHREEQC (Parkhurst and Appelo 2013). The main goal of the interface is to take advantages and capabilities of mentioned codes, providing a numerical platform that can efficiently simulate a wide number of multiphysics problems coupled with geochemistry. iCP is written in Java and uses the IPhreeqcC++ dynamic library and the COMSOL Java-API (Nardi et al. 2014).

COMSOL Multiphysics (COMSOL 2022) is a commercial finite element code with a set of predefined physics interfaces widely used in various industrial applications. Alternatively, input of user's own partial differential equation is possible through the Partial Differential Equation (PDE) Interface. Then, COMSOL Multiphysics provides the coupling between the individual equations by either a fully coupled or segregated solver. The polynomial order of the finite element base (shape) functions can be chosen by the user. The default settings were used, linear for the transport equations and quadratic for the flow.

PHREEQC (Parkhurst and Appelo 2013) is a free general geochemical code developed by U.S. Geological Survey for simulating chemical reactions and transport processes in water. The program is based on equilibrium chemistry of aqueous solutions interacting with minerals, gases, solid solutions, exchangers, and sorption surfaces, which accounts for the original acronym: pH-REdox-Equilibrium (Parkhurst and Appelo 2013). Although a broader functionality is also included in PHREEQC, the equilibrium reactions were only used in this benchmark. The program uses a modified Newton-Raphson method to solve the simultaneous nonlinear equations (Parkhurst and Appelo 1999).

The coupling between the physical processes solved by COMSOL Multiphysics (water flow, water-phase species

transport, and gas-phase diffusion) and chemical reactions solved by PHREEQC is realised by the operator splitting method, also called Sequential Non Iterative Approach (SNIA) (Nardi et al. 2022). In a single time step controlled by the iCP user setting, the two models run one by one. COMSOL is based on the Finite Element Method and, therefore, the variables (i.e. concentrations) are calculated at the nodes of the mesh. These node concentrations are transferred to PHREEQC at every communication step to perform the chemical calculations. The model can be divided in chemical domains with different geochemical systems (i.e. with different mineralogy). Nodes at the interface of two or more chemical domains are calculated for each domain. The resulting concentrations are arithmetically averaged, assigned to the node and transferred to COMSOL (Nardi et al. 2022). The chemical step is parallelized by grouping sets of nodes and solving them on different PHREEQC instances. These PHREEQC instances are solved on different executable threads by means of the internal thread library of Java (shared memory approach) (Nardi et al. 2022).

iCP interface make additional physics interface available in COMSOL Multiphysics platform named "Molal solute transport by advection-diffusion". This interface was used in transport simulations of a number of chemical species and could be coupled to other interfaces in COMSOL Multiphysics. For modelling of water mass balance in the porous media and CO₂ transport in predefined test cases the physics interfaces "Richard's equation" and "Transport of Diluted Species in Porous Media" were used additionally respectively. Parameter input to the latter interfaces as well as the solver settings and timestepping for iCP requires definition by the user, therefore these inputs may not be unique among the iCP users. Running PHREEQC for the benchmark cases required a non-standard use of its input control. It is made to be run with the full set of species and reactions in the thermodynamic database. Therefore, a reduced database file was created, containing only the species and reactions in the respective Test Case and water. The activity coefficients were not input explicitly but evaluated by PHREEQC via the implemented Debye-Hückel model. The resulting values were checked against the prescribed ones. Contrary to most other reactive transport codes, PHREEQC requires enforcing charge balance at all times, including the initial conditions and in all boundary waters. To compare the results of INVERSE-FADES-CORE V2, TOUGHREACT and DuMuₓ with those of iCP, charge balance was enforced by adding two dummy charged unreactive species. Cl⁻ and Na⁺ were selected as the dummy unreactive species with appropriate concentrations to ensure charge balance. These species are included in the geochemical system, contribute to the

<figure><img src="image_5.png" /></figure>
total ionic strength of the aqueous dissolved system and affect the chemical activities of the dissolved species.

## Comparison and discussion of the results

In this section, we conduct a comparative analysis of results obtained from different numerical codes. To maintain conciseness, we present only a subset of the computed results from participants, with the full set available as supplementary data. In the figures, the teams are identified by the abbreviated name of the code used as follows: UDC is represented as F-CORE, UPPA as DuMu$^x$, JU as T-REACT, LEI as iCP-1, and TUL as iCP-2. All test cases were completed by each code, except for test cases TC2B and TC3B, which were only performed by F-CORE and iCP-1.

### Case TC1. Hydration and conservative tracer.

Case TC1 deals with single-phase variably-saturated water flow and conservative solute transport. Bentonite hydration was simulated by prescribing the liquid pressure equal to 100 kPa at the right boundary ($r = 1.14$ m). A constant temperature of 25 °C was considered in this case. The total time of simulation was 18 years. The computed water content increases with time (Fig. 3a) due to the inflow of water from the right boundary ($r = 1.14$ m). The concentration of the tracer decreases near the right boundary because the concentration of the inflow water ($1.3 \times 10^{-5}$ mol/L) is smaller than the initial concentration ($1.6 \times 10^{-1}$ mol/L) (see Fig. 3b). The tracer concentration front displaces due to advection and spreads due to diffusion and dispersion.

Case TC1 was simulated by 5 teams. The water contents computed by different teams agree for the most part. However, there are small discrepancies of about 0.005 near the left boundary ($r = 0.45$ m) at $t = 18$ years (Fig. 3a). The spatial distribution of the computed tracer concentrations is similar for most teams. However, there are small differences at the tracer front. These differences are possibly related to differences in hydrodynamic dispersion because advection is the dominant transport mechanism. The concentrations computed by F-CORE and iCP-1 coincide in the case of no water flow (only solute diffusion) (see Fig. S1 in Supplementary Material). Model sensitivity runs show that the tracer concentrations are sensitive to changes in the longitudinal dispersity, the diffusion coefficient, and the water flow (see Fig. S2 in Supplementary Material).

### Case TC2 and TC2B. No hydration and gas diffusion

Case TC2 deals with CO₂(g) gas diffusion and aqueous chemical reactions by considering H⁺, CO₂(aq), and CO₃²⁻ as

secondary species in the reactive transport model. The partial pressure of CO₂(g) was prescribed at 0.1 bar at the right boundary ($r = 1.14$ m). CO₂(g) partial pressure in the bentonite increases with time (Fig. 3d). The partial pressure of CO₂(g) becomes uniform after 2 days and equal to the prescribed pressure. Aqueous CO₂(aq) is in equilibrium with gaseous CO₂(g), and therefore, the increase in CO₂(g) pressure leads to an increase in the concentration of dissolved CO₂(aq) and a decrease in pH (Fig. 3c, e).

Case TC2 was simulated by 5 teams. The CO₂(g) partial pressure computed by different teams agree for the most part. However, there are small discrepancies near the left boundary ($r = 0.45$ m) in the computed pH and CO₂(aq) at early times ($t < 0.1$ days). These discrepancies are related to the limited accuracy of the numerical solution provided by the numerical grid. The largest discrepancy in CO₂(aq) is approximately equal to $1.5 \times 10^{-4}$ mol/L at $t = 0.1$ days.

Case TC2B is based on case TC2 and includes additional secondary species to ensure charge balance in the reactive transport model (see Table 1). The computed results for cases TC2B and TC2 are similar (Fig. 3f), except for the computed pH. The computed pH in case TC2B is smaller than that of case TC2. Case TC2B test case was simulated by 2 teams (Fig. 3f). Similar to case TC2, small discrepancies in pH are found in case TC2B near the left boundary at early times.

Computed CO₂(g) partial pressure increases with time until $t = 0.01$ years when the pressure becomes equal to the prescribed pressure at the boundary (Fig. 4a). The computed pH shows a trend opposite to that of CO₂ because it decreases until $t = 0.01$ years to 7.72. The results computed for case TC2 by 5 teams involved agree for the most part (Fig. 4c).

### Case TC3. No hydration, gas diffusion and mineral dissolution/precipitation

Case TC3 deals with gas diffusion and calcite dissolution/precipitation. Water saturation remains constant and equal to 60% while the CO₂(g) partial pressure was prescribed equal to 0.1 bar at the right boundary ($r = 1.14$ m). Case TC3 takes into account calcite dissolution/precipitation. Concentrations of dissolved species decrease near the right boundary because the concentrations prescribed at the right boundary are smaller than the initial concentrations. Similar to case TC2, the partial pressure of CO₂(g) becomes uniform and equal to the prescribed gas pressure at $t = 2$ days (Fig. 5a). The increase in the concentration of CO₂(aq) leads to a decrease in pH (Fig. 5c). In addition, calcite precipitates in a thin band of a few mm near the right boundary (Fig. 5b).

Case TC3 was simulated by 5 teams. The CO₂(g) partial pressure computed by different teams agree for
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Benchmarking of the spatial distribution of computed: **a** Water contents for case TC1 at *t* = 1 and *t* = 18 years; **b** Tracer concentrations for case TC1 at *t* = 1 and *t* = 18 years; **c** pH for cases TC2 at *t* = 0.1, *t* = 0.35 and *t* = 2 days; **d** Partial pressure of CO₂(g) for case</figcaption></figure>

TC2 at *t* = 0.1, *t* = 0.35 and *t* = 2 days; **e** Concentration of dissolved CO₂(aq) for case TC2 at *t* = 0.1, *t* = 0.35 and *t* = 2 days; and **f** pH for cases TC2B at *t* = 0.1, *t* = 0.35 and *t* = 2 days

the most part (Fig. 5a). However, there are some small discrepancies near the left boundary which translate also into small discrepancies in the computed pH at early times. The spatial distribution of the computed calcium concentrations is similar for most teams (Fig. 5d). However, small differences are found at the concentration front similar to those of case TC1. Computed calcite precipitation and pH spatial distribution coincides for most teams.

Case TC3B is similar to case TC3 and incorporates additional secondary species (see Table 1). The computed results in case TC3B show calcite dissolution instead of calcite precipitation in case TC3. This test case was simulated by 2 teams. The computed CO₂(g) partial pressure, dissolved CO₂(aq) concentration, pH and dissolved Ca²⁺ concentrations results for cases TC3B and TC3 are similar (Fig. 5e, f). There are small discrepancies in the computed pH and calcite.

## Case TC4. Hydration, gas diffusion and mineral dissolution/precipitation

Case TC4 deals with liquid flow and gas diffusion and considers the geochemical system of case TC3. The computed water content increases with time due to the water inflow which was described in case TC1 (Fig. 3a). The effective diffusion coefficient of the CO₂(g) decreases at the right boundary due to the water inflow through the right boundary. The diffusion of CO₂(g) is slower than in cases TC2 and TC3 due to the increase of water saturation. The time needed for the partial pressure of CO₂(g) to become uniform in case TC4 is about 20 days (Fig. 6a). This time is much larger than the times required in cases TC2 and TC3. The computed concentration of dissolved Ca²⁺ decreases near the right boundary because the concentration of the boundary water is smaller than the initial concentration. In addition, calcite precipitates in the bentonite and especially in a thin band of a few mm near the right boundary at early times (Fig. 6b). Calcite precipitation consumes dissolved
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Benchmarking of the time evolution of computed: **a** Partial pressure of CO₂(g) for case TC2 at x = 1.12 m, 1 m and 0.6 m; **b** Partial pressure of CO₂(g) for case TC4 at x = 1.12 m and 0.6 m; **c** pH for</figcaption></figure>

case TC2 at the distances x = 1.12 m, 1 m and 0.6 m; and **d** pH for case TC4 at the distances x = 1.12 m, 1 m and 0.6 m

calcium (Fig. 6c). pH decreases in the bentonite while CO₂(g) diffuses into the bentonite. However, pH increases when calcite precipitates near the right boundary. The calcite precipitation front in case TC4 penetrates into the bentonite about 35 cm at the end of the simulation (18 years). The computed pH decreases in the bentonite as the dissolved CO₂(g) diffuses into the bentonite and increases near the right boundary due to calcite precipitation (Fig. 6d).

Case TC4 was simulated by 5 teams. The computed CO₂(g) partial pressures by different teams follow the same trends (Fig. 6a). However, there are some discrepancies in the diffusion profile. These discrepancies are caused by differences in the formulations implemented in different codes for gas flow at the transition from unsaturated to

nearly saturated conditions. The computed CO₂(g) partial pressures are very sensitive to mesh discretization. The gas diffusion profiles computed by different teams with a uniform mesh show significant differences (see Fig. S3 in Supplementary Material). These differences are less relevant when a more refined grid is used.

The concentrations of dissolved Ca²⁺ computed by different teams are similar. However, there are small differences at the concentration front similar to those reported for previous cases. Computed calcite precipitation and pH spatial distribution trends coincides for most teams.

The diffusion of CO₂(g) is fast. Steady and uniform pressures are reached in less than a year. Model results
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Benchmarking of the spatial distribution of computed: **a** Partial pressure of CO₂(g) for TC3 at *t* = 0.35 years and *t* = 0.5 years; **b** Cumulative calcite volume fraction (positive for precipitation) for</figcaption></figure>

TC3 at *t* = 0.5 years; **c** pH for TC3 at *t* = 0.5 years; **d** Dissolved calcium for TC3 at 0.5 years; **e** Cumulative calcite volume fraction (positive for precipitation) for TC3B at *t* = 0.5 years; and **f** pH for TC3B at *t* = 0.5 years

show a fluctuation in the partial pressure between 10⁻⁴ and 10⁻³ years which slows down the diffusion of CO₂(g) (see Fig. 4b). This fluctuation can be seen better when pressures are plotted versus log-time (see Fig. S4 in Supplementary Material). The computed pH increases after 0.01 years (Fig. 4d) due to the pH front associated with calcite precipitation near the right boundary (see Fig. S4 in Supplementary Material).

### Case TC5. Hydration, CO₂(g) diffusion and calcite and gypsum at equilibrium

Case TC5 is an extension of case TC4 which includes calcite and gypsum dissolution/precipitation. Similar to previous cases, the computed water content increases with time due to the water inflow. The computed CO₂(g) partial pressure in the bentonite increases with time. Aqueous CO₂(aq) is in equilibrium with gaseous CO₂(g) and therefore, the increase in CO₂(g) pressure leads to an increase in the concentration

of dissolved CO₂(aq) and a decrease in pH. The computed CO₂(g) partial pressures and dissolved CO₂(aq) diffusion into bentonite follow similar patterns to those of Case TC4 (Fig. 6a). The computed concentrations of dissolved calcium and sulphate are strongly affected by dilution of the inflow water and by calcite and gypsum precipitation/dissolution. The breakthrough curves of dissolved calcium and sulphate are not smooth and show broken lines (Fig. 7a, b). Calcite precipitates while gypsum dissolves in the bentonite near the right boundary. Mineral dissolution/precipitation reactions at the final time take place in a 10 cm thick band (Fig. 7c, d). The computed concentration of dissolved Ca²⁺ decreases near the right boundary because the concentration of the boundary water is smaller than the initial concentration and due to calcite precipitation (Fig. 7a). However, it is counteracted by the release of Ca²⁺ from gypsum dissolution. Computed dissolved sulphate decreases also near the right boundary because its concentration in the boundary water is smaller than the initial concentration (Fig. 7b). On the other
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Case TC4 benchmarking of the spatial distribution of computed: a Dissolved CO₂(g) concentrations at *t* = 0.35 days, 3.5 days and *t* = 20 days; b Cumulative calcite volume fraction (positive for</figcaption></figure>

hand, dissolved sulphate increases near the right boundary due to gypsum dissolution. Computed pH increases near the right boundary where calcite precipitates. Case TC5 was simulated by 5 teams. The curves of computed concentration of dissolved CO₂(aq) follow the same trends. However, there are some discrepancies which are caused by differences in the formulations implemented in different codes for gas flow at the transition from unsaturated to nearly saturated conditions similar to those mentioned in Case TC4. The concentrations of dissolved Ca²⁺ and sulphate computed with different codes are similar. However, there are small differences at the fronts similar to those reported for previous cases. The fronts of pH and total cumulative calcite precipitation and gypsum dissolution are similar for all codes.

precipitation) at *t* = 20 days and *t* = 18 years; c Dissolved calcium at *t* = 20 days and *t* = 18 years; and d pH at *t* = 20 days and *t* = 18 years

**Case TC6. Hydration, CO₂(g) diffusion, calcite and gypsum at equilibrium and cation exchange**

Case TC6 is based on Case TC5 and includes cation exchange reactions. The concentrations of dissolved species are affected by calcite and gypsum precipitation/ dissolution, dilution from the hydration water and cation exchange reactions. Similar to Case TC5, the concentrations of dissolved species in TC6 decrease near the right boundary (see Fig. S5 in Material Supplementary). The computed pH increases at the locations where calcite precipitates and decreases at a point located 1.7 cm from the right boundary where calcite dissolves (Fig. 8d). Gypsum dissolution takes place during all the simulation time (Fig. 8b) while calcite shows a dissolution/precipitation front near the boundary (Fig. 8a). This front is strongly linked with dissolved Ca gets into the exchange complex at the expense of exchanged Na
<figure><img src="image_10.png" /><figcaption>**Fig. 7** TC5 benchmarking of the spatial distribution of computed: a Dissolved calcium at *t* = 20 days and *t* = 18 years; b Dissolved sulphate at *t* = 20 days and *t* = 18 years; c Cumulative calcite volume</figcaption></figure>

fraction (positive for precipitation) at *t* = 18 years; and d Cumulative gypsum volume fraction (positive for precipitation) at *t* = 18 years

(X-Na). In fact, the concentration of exchanged Ca increases and that of exchanged Na decreases. The concentrations of exchanged K and Mg show small changes (Fig. 8c).

Case TC6 was simulated with 5 codes. The computed concentrations of dissolved species of different codes follow similar trends. However, there are small discrepancies in the fronts which are caused by differences in the formulations implemented in different codes for gas flow at the transition from unsaturated to nearly saturated conditions similar to those seen in previous cases. The fronts of pH and total cumulative calcite precipitation and gypsum dissolution are similar for all codes (Fig. 8).

## Discussions

One of the main challenges of the benchmark had to do with ensuring that all codes and teams solved the benchmark test cases for the same flow and transport conditions, reactions

and initial and boundary conditions. To help ensuring the consistency among codes, the benchmark was designed in test cases of increasing complexity. Consistency was first warranted for conservative solute transport (Case TC1) and then for gas diffusion (TC2). Dirichlet conditions were used for gas partial pressures and concentrations of dissolved species at the inflow boundary to ensure consistency in boundary conditions. Some of the discrepancies were related to differences among codes regarding: 1) Numerical methods for spatial discretization (finite elements, finite differences and finite volumes) which may lead to oscillations and numerical dispersion, 2) Numerical schemes for coupled transport and chemistry, and 3) Formulation and parameters for activity coefficients.

The computed concentrations of dissolved species are very sensitive to water flow, solute dispersion and diffusion. Discrepancies are found in conservative tracer concentrations in Case TC1 due to differences in hydrodynamic
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Case TC6 benchmarking of the spatial distribution of computed: **a** Cumulative calcite volume fraction (positive for precipitation and negative for dissolution) at *t* = 5 years; the inner figure shows a ampliation of the cumulative calcite precipitation for 1 m < *x* < 1.14</figcaption></figure>

m; **b** Cumulative gypsum volume fraction (positive for precipitation and negative for dissolution) at *t* = 5 years; **c** Concentrations of exchange cations at *t* = 5 years; and **d** pH at *t* = 5 years

dispersion because concentration profiles fully agree when water flow is disregarded. Discrepancies in the concentration of the tracer in Case TC1 are similar to those of a “dummy” tracer in test Case TC5 which is defined as the difference between the concentrations of Ca²⁺ and SO₄²⁻ (Fig. S6).

The CO₂(g) partial pressures computed by different codes in Case TC2 follow the same trends. However, there are discrepancies in gas diffusion profiles at early times near the no-flow boundary. These discrepancies, which vanish at later times, are attributed to differences in the formulations implemented in different codes for gas flow at the transition from unsaturated to nearly saturated conditions. The

discrepancies in CO₂(g) pressures lead to small discrepancies in computed pH.

Some codes such as iCP require the preservation of charge balance. This requirement led to the definition of Cases TC2B and TC3B, in which a more complete suite of ions and aqueous complexes were considered to ensure and preserve charge balance. Although the aqueous chemical systems of Cases TC2B and TC3B are similar to those of Cases TC2 and TC3, there are important differences in the computed pH and calcite dissolution/precipitation between Cases TC2B and TC3B and Cases TC2 and TC3. Computed pH for Cases TC2B and TC3B is smaller than those of Cases
TC2 and TC3. In addition, calcite dissolves in Cases TC2B and TC3B while it precipitates in Cases TC2 and TC3. Charge balance in Cases TC5 and TC6 was preserved by adding Na⁺ and Cl⁻ to the chemical system.

The CO₂(g) partial pressures computed in Case TC4 are very sensitive to mesh discretization. The gas diffusion profiles computed with different codes with a uniform mesh show significant differences. The gas diffusion coefficient decreases significantly when water flows into the column. Most codes exhibit large discretization errors for modeling gas diffusion when the mesh is not sufficiently refined. Several discretization schemes were tested. Grid refinement near the inflow boundary was performed until convergence of the numerical solution was reached. A mesh-independent solution was achieved by refining the mesh near the inflow boundary with a transition zone of gradually increasing grid size. All solutions agree when the refined grid is used. The optimum refined grid of Case TC4 was used also for the rest of the test cases. The refined grid helped to overcome the discrepancies in Cases TC5 and TC6.

All codes use the same thermodynamic data base for solving the benchmark test cases. However, there are differences among codes about the parameters and the equation used to calculate the activity coefficients of dissolved species. Activity coefficients are assumed to be constant in DuMuₓ. iCP and INVERSE-FADES-CORE V2 use the extended Debye–Hückel (Eq. (18)) while TOUGHREACT uses the Helgeson–Kirkham–Flowers (HKF) extended Debye–Hückel equation (Appelo et al. 2014). The computed concentrations of dissolved species are very sensitive to the values of the activity coefficients, especially in the most reactive test cases (TC5 and TC6) in which the activity coefficients vary with time significantly. Figure S7 shows the concentrations of dissolved Ca²⁺ computed with several options for the calculation of the activity coefficients. One can see clearly that the concentrations of dissolved Ca²⁺ calculated with constant activities deviate from those calculated with time-varying activities. Tables S1 to S4, provided as in supplementary material, list the initial and final concentrations and activity coefficients calculated by different codes for Cases TC5 and TC6. Although the initial activities of the codes are similar, the final activities show significant differences because activity coefficients vary with time near the boundary (see Fig. S8).

## Conclusion and perspectives

A multiphase flow and reactive transport benchmark for radioactive waste disposal has been presented. The benchmark model is based on a simplification of THCM model of the FEBEX in situ test presented by Samper et al. (2018a). The main features of the numerical model include: (1) 1D model;

(2) Unsaturated bentonite; (3) Gases (dry air and CO₂(g)); (4) Minerals (calcite and gypsum); and (5) Geochemical reactions (aqueous complexation, mineral dissolution/precipitation, cation exchange and gas dissolution/ex-solution). The following multiphase flow and reactive transport codes were used in the benchmark: INVERSE-FADES-CORE V2, DuMuₓ, TOUGHREACT and iCP. A total of 6 test cases were studied with increasing complexity in flow and chemistry conditions ranging from conservative tracer migration with hydration under variably unsaturated conditions to multiphase flow with calcite and gypsum dissolution/precipitation and cation exchange reactions. All codes provide results with similar trends. Some discrepancies are found in conservative tracer concentrations due to differences in hydrodynamic dispersion. The computed CO₂(g) partial pressures by different codes follow the same trends. However, there are some discrepancies in the gas diffusion profiles at early times near the left boundary (no flow boundary). These discrepancies which vanish at later times are caused by differences in the formulations implemented in different codes for gas flow at the transition from unsaturated to nearly saturated conditions. The discrepancies in CO₂(g) pressures translate into small discrepancies in computed pH. The computed CO₂(g) partial pressures are very sensitive to mesh discretization. The gas diffusion profiles computed with different codes with a uniform mesh show significant differences. All solutions tend to agree when a more refined grid is used. The concentrations of dissolved Ca computed with different codes are similar. However, there are small differences at the concentration fronts. All codes agree in the spatial distribution of pH and calcite precipitation. The diffusion of CO₂(g) is fast. Steady and uniform pressures are reached in less than a year. Computed CO₂(g) partial pressures show a fluctuation between 10⁻⁴ and 10⁻³ years which slows down the diffusion of CO₂(g). This fluctuation is associated with chemical reactions involving CO₂ and is shown more clearly when pressures are plotted versus log-time. The computed pH increases after 0.01 years due to the pH front associated with calcite precipitation near the right boundary. There are some discrepancies in the computed concentrations of the dissolves species which are due to differences in the chemical activities of dissolved species. Some codes use different formulations of the Debye–Hückel and others assume the activities remain constant in time. Special attention was paid to ensure that all codes use the same activities. The differences in the codes caused by differences in activities have a significant impact on the fronts of dissolved Ca and sulfate.

The results of this multiphase flow code benchmark for radioactive waste disposal enhance the confidence on the use of multiphase reactive transport models. Future work should be done to extend the benchmarking to non-isothermal conditions, including vapor generation, vapor and gas transport and metallic canister corrosion.
# Appendix

the initial mineral volume fractions, and the initial concentrations of exchanged ions for all test cases.

This appendix includes Tables 4, 5 and 6, which detail the initial compositions of the bentonite and granite pore water,

<table><caption>**Table 4** Equilibrated initial pore water composition and activities, initial gas partial pressure in the bentonite and granite for TC1, TC2 and TC2B</caption><thead><tr><th rowspan="2">TC1</th><th colspan="2">Bentonite</th><th colspan="2">Granite (r = 1.14 m)</th></tr><tr><th>Concentration</th><th>Activity coefficient</th><th>Concentration</th><th>Activity coefficient</th></tr></thead><tbody><tr><td>Cl<sup>-</sup></td><td>1.6 × 10<sup>-1</sup> mol/L</td><td></td><td>1.3 × 10<sup>-5</sup> mol/L</td><td></td></tr><tr><td colspan="5">TC2</td></tr><tr><td>H<sup>+</sup></td><td>1.799 × 10<sup>-10</sup> mol/L</td><td>0.9921</td><td>–</td><td>–</td></tr><tr><td>CO<sub>2</sub>(aq)</td><td>3.388 × 10<sup>-8</sup> mol/L</td><td>1.0000</td><td>–</td><td>–</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>2.295 × 10<sup>-5</sup> mol/L</td><td>–</td><td>–</td><td>–</td></tr><tr><td colspan="5">TC2B</td></tr><tr><td>H<sup>+</sup></td><td>1.847 × 10<sup>-10</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>8.583 × 10<sup>-5</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>Cl<sup>-</sup></td><td>10 × 10<sup>-10</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>Ca<sup>2+</sup></td><td>9.656 × 10<sup>-11</sup> mol/L</td><td>0.935</td><td>–</td><td>–</td></tr><tr><td>Na<sup>+</sup></td><td>1.91 × 10<sup>-4</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>CO<sub>2</sub>(aq)</td><td>3.43 × 10<sup>-8</sup> mol/L</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>CaCl<sup>+</sup></td><td>4.63 × 10<sup>-21</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>CaCO<sub>3</sub></td><td>3.256 × 10<sup>-12</sup> mol/L</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>Ca(OH)<sup>+</sup></td><td>8.389 × 10<sup>-14</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>NaCl</td><td>5.84 × 10<sup>-15</sup> mol/L</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>CaCl<sub>2</sub></td><td>2.00 × 10<sup>-31</sup> mol/L</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>2.325 × 10<sup>-5</sup> mol/L</td><td>0.935</td><td>–</td><td>–</td></tr><tr><td>Ca(HCO<sub>3</sub>)<sup>+</sup></td><td>9.75 × 10<sup>-14</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>Na(HCO<sub>3</sub>)</td><td>8.91 × 10<sup>-9</sup> mol/L</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>Na(CO<sub>3</sub>)<sup>-</sup></td><td>7.731 × 10<sup>-5</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>Na(OH)</td><td>1.839 × 10<sup>-9</sup> mol/L</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>OH<sup>-</sup></td><td>5.601 × 10<sup>-5</sup> mol/L</td><td>0.983</td><td>–</td><td>–</td></tr><tr><td>CO<sub>2</sub>(g)</td><td>10<sup>-6</sup> bar</td><td></td><td>0.1 bar</td><td></td></tr></tbody></table>
<table><caption>**Table 5** Equilibrated initial pore water composition and activities, initial mineral volume fractions and gas partial pressure in the bentonite and granite for TC3, TC3B and TC4</caption><thead><tr><th rowspan="2">TC3 and TC4</th><th colspan="2">Bentonite</th><th colspan="2">Granite (r = 1.14 m)</th></tr><tr><th>Concentration</th><th>Activity coefficient</th><th>Concentration</th><th>Activity coefficient</th></tr></thead><tbody><tr><td>H<sup>+</sup></td><td>1.77 × 10<sup>-9</sup> mol/L</td><td>0.8122</td><td>3.201 × 10<sup>-8</sup> mol/L</td><td>0.9819</td></tr><tr><td>CO<sub>2</sub>(aq)</td><td>3.374 × 10<sup>-8</sup> mol/L</td><td>1.0042</td><td>3.388 × 10<sup>-3</sup> mol/L</td><td>1.0000</td></tr><tr><td>Ca<sup>2+</sup></td><td>2.195 × 10<sup>-2</sup> mol/L</td><td>0.4405</td><td>4.799 × 10<sup>-5</sup> mol/L</td><td>0.9295</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>7.701 × 10<sup>-7</sup> mol/L</td><td>0.4447</td><td>7.707 × 10<sup>-5</sup> mol/L</td><td>0.9296</td></tr><tr><td>CO<sub>2(g)</sub></td><td>10<sup>-6</sup> bar</td><td></td><td>0.1 bar</td><td></td></tr><tr><td>Calcite</td><td>0.01 vf</td><td></td><td>—</td><td></td></tr><tr><td colspan="5">TC3B</td></tr><tr><td>H<sup>+</sup></td><td>1.641 × 10<sup>-9</sup> mol/L</td><td>0.782</td><td>3.615 × 10<sup>-8</sup> mol/L</td><td>0.938</td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>1.41 × 10<sup>-5</sup> mol/L</td><td>0.775</td><td>3.22 × 10<sup>-3</sup> mol/L</td><td>0.937</td></tr><tr><td>Cl<sup>-</sup></td><td>4.37 × 10<sup>-2</sup> mol/L</td><td>0.775</td><td>9.99 × 10<sup>-11</sup> mol/L</td><td>0.937</td></tr><tr><td>O<sub>2</sub>(aq)</td><td>1.96 × 10<sup>-31</sup> mol/L</td><td>1.006</td><td>6.41 × 10<sup>-40</sup> mol/L</td><td>1.000</td></tr><tr><td>Na<sup>+</sup></td><td>9.92 × 10<sup>-11</sup> mol/L</td><td>0.782</td><td>3.16 × 10<sup>-3</sup> mol/L</td><td>0.938</td></tr><tr><td>Ca<sup>2+</sup></td><td>2.17 × 10<sup>-2</sup> mol/L</td><td>0.382</td><td>4.63 × 10<sup>-5</sup> mol/L</td><td>0.775</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>1.03 × 10<sup>-6</sup> mol/L</td><td>0.386</td><td>5.37 × 10<sup>-6</sup> mol/L</td><td>0.775</td></tr><tr><td>CO<sub>2</sub>(aq)</td><td>3.13 × 10<sup>-8</sup> mol/L</td><td>1.006</td><td>2.29 × 10<sup>-4</sup> mol/L</td><td>1.00</td></tr><tr><td>OH<sup>-</sup></td><td>1.004 × 10<sup>-5</sup> mol/L</td><td>0.775</td><td>3.15 × 10<sup>-7</sup> mol/L</td><td>0.937</td></tr><tr><td>H<sub>2</sub>(aq)</td><td>1.90 × 10<sup>-31</sup> mol/L</td><td>1.006</td><td>3.36 × 10<sup>-27</sup> mol/L</td><td>1.00</td></tr><tr><td>CaCl<sup>+</sup></td><td>1.84 × 10<sup>-4</sup> mol/L</td><td>0.782</td><td>1.84 × 10<sup>-15</sup> mol/L</td><td>0.938</td></tr><tr><td>CaCO<sub>3</sub></td><td>5.46 × 10<sup>-6</sup> mol/L</td><td>1.0006</td><td>2.48 × 10<sup>-7</sup> mol/L</td><td>1.00</td></tr><tr><td>Ca(OH)<sup>+</sup></td><td>1.368 × 10<sup>-6</sup> mol/L</td><td>0.782</td><td>1.87 × 10<sup>-10</sup> mol/L</td><td>0.938</td></tr><tr><td>NaCl</td><td>8.26 × 10<sup>-13</sup> mol/L</td><td>1.006</td><td>8.76 × 10<sup>-14</sup> mol/L</td><td>1.00</td></tr><tr><td>CaCl<sub>2</sub></td><td>2.16 × 10<sup>-6</sup> mol/L</td><td>1.006</td><td>7.20 × 10<sup>-26</sup> mol/L</td><td>1.00</td></tr><tr><td>Ca(HCO<sub>3</sub>)<sup>+</sup></td><td>1.46 × 10<sup>-6</sup> mol/L</td><td>0.782</td><td>1.45 × 10<sup>-6</sup> mol/L</td><td>0.938</td></tr><tr><td>Na(HCO<sub>3</sub>)</td><td>4.75 × 10<sup>-16</sup> mol/L</td><td>1.006</td><td>5.03 × 10<sup>-6</sup> mol/L</td><td>1.00</td></tr><tr><td>Na(CO<sub>3</sub>)<sup>-</sup></td><td>7.44 × 10<sup>-16</sup> mol/L</td><td>0.775</td><td>2.45 × 10<sup>-7</sup> mol/L</td><td>0.937</td></tr><tr><td>Na(OH)</td><td>1.06 × 10<sup>-16</sup> mol/L</td><td>1.006</td><td>1.55 × 10<sup>-10</sup> mol/L</td><td>1.00</td></tr><tr><td>CO<sub>2(g)</sub></td><td>10<sup>-6</sup> bar</td><td></td><td>0.1 bar</td><td></td></tr><tr><td>Calcite</td><td>0.01 vf</td><td></td><td>—</td><td></td></tr></tbody></table>
<table><caption>**Table 6** Equilibrated initial pore water composition and activities, initial mineral volume fractions (vf) and gas partial pressure in the bentonite and granite for TC5 and TC6</caption><thead><tr><th rowspan="2">TC5</th><th colspan="2">Bentonite</th><th colspan="2">Granite (r = 1.14 m)</th></tr><tr><th>Concentration</th><th>Activity coefficient</th><th>Concentration</th><th>Activity coefficient</th></tr></thead><tbody><tr><td>H<sup>+</sup></td><td>1.030 × 10<sup>-9</sup> mol/L</td><td>0.8111</td><td>3.331 × 10<sup>-8</sup> mol/L</td><td>0.9745</td></tr><tr><td>CO<sub>2</sub>(aq)</td><td>3.369 × 10<sup>-8</sup> mol/L</td><td>1.0055</td><td>3.392 × 10<sup>-3</sup> mol/L</td><td>1.0000</td></tr><tr><td>Ca<sup>2+</sup></td><td>7.124 × 10<sup>-3</sup> mol/L</td><td>0.4590</td><td>4.799 × 10<sup>-5</sup> mol/L</td><td>0.9028</td></tr><tr><td>SO<sub>4</sub><sup>2-</sup></td><td>1.715 × 10<sup>-2</sup> mol/L</td><td>0.4386</td><td>7.900 × 10<sup>-5</sup> mol/L</td><td>0.9021</td></tr><tr><td>Na<sup>+</sup></td><td>2.01 × 10<sup>-2</sup> mol/L</td><td>0.8311</td><td>2.19 × 10<sup>-4</sup> mol/L</td><td>0.9745</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>2.309 × 10<sup>-6</sup> mol/L</td><td>0.4386</td><td>7.465 × 10<sup>-5</sup> mol/L</td><td>0.9021</td></tr><tr><td>CO<sub>2</sub>(g)</td><td>10<sup>-6</sup> bar</td><td></td><td>0.1 bar</td><td></td></tr><tr><td>Calcite</td><td>0.01 vf</td><td></td><td>—</td><td></td></tr><tr><th colspan="5">TC6</th></tr><tr><td>H<sup>+</sup></td><td>1.183 × 10<sup>-9</sup> mol/L</td><td>0.7065</td><td>3.443 × 10<sup>-8</sup> mol/L</td><td>0.9707</td></tr><tr><td>CO<sub>2</sub>(aq)</td><td>3.288 × 10<sup>-8</sup> mol/L</td><td>1.0306</td><td>3.392 × 10<sup>-3</sup> mol/L</td><td>1.0001</td></tr><tr><td>Ca<sup>2+</sup></td><td>1.163 × 10<sup>-2</sup> mol/L</td><td>0.2831</td><td>4.879 × 10<sup>-5</sup> mol/L</td><td>0.8892</td></tr><tr><td>Mg<sup>2+</sup></td><td>3.124 × 10<sup>-2</sup> mol/L</td><td>0.2831</td><td>1.300 × 10<sup>-6</sup> mol/L</td><td>0.8892</td></tr><tr><td>Na<sup>+</sup></td><td>1.578 × 10<sup>-1</sup> mol/L</td><td>0.7065</td><td>3.800 × 10<sup>-4</sup> mol/L</td><td>0.9707</td></tr><tr><td>K<sup>+</sup></td><td>2.086 × 10<sup>-3</sup> mol/L</td><td>0.7065</td><td>7.800 × 10<sup>-6</sup> mol/L</td><td>0.9707</td></tr><tr><td>SO<sub>4</sub><sup>2-</sup></td><td>3.039 × 10<sup>-2</sup> mol/L</td><td>0.249</td><td>7.900 × 10<sup>-5</sup> mol/L</td><td>0.8884</td></tr><tr><td>Cl<sup>-</sup></td><td>1.848 × 10<sup>-1</sup> mol/L</td><td>0.6918</td><td>1.784 × 10<sup>-4</sup> mol/L</td><td>0.9706</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>4.038 × 10<sup>-6</sup> mol/L</td><td>0.249</td><td>7.151 × 10<sup>-5</sup> mol/L</td><td>0.8884</td></tr><tr><td>CO<sub>2</sub>(g)</td><td>10<sup>-6</sup> bar</td><td></td><td>0.1 bar</td><td></td></tr><tr><td>Calcite</td><td>0.01 vf</td><td></td><td>—</td><td></td></tr><tr><td>Gypsum</td><td>0.0015 vf</td><td></td><td>—</td><td></td></tr><tr><td>Exchanged Ca<sup>2+</sup></td><td>19.65 meq/100g</td><td></td><td>—</td><td></td></tr><tr><td>Exchanged Mg<sup>2+</sup></td><td>54.87 meq/100g</td><td></td><td>—</td><td></td></tr><tr><td>Exchanged Na<sup>+</sup></td><td>25.77 meq/100g</td><td></td><td>—</td><td></td></tr><tr><td>Exchanged K<sup>+</sup></td><td>2.46 meq/100g</td><td></td><td>—</td><td></td></tr></tbody></table>

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s12665-024-11887-6.

**Acknowledgements** This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 847593 (WP DONUT). The work of UDC was funded by ENRESA, the Spanish Ministry of Science and Innovation (PID2023-153202OB-I00) and the Galician Regional Government (Grant ED431C2021/54). The work of M. El Ossmani has been partly supported by the Carnot ISIFoR Institute and "la Région Nouvelle-Aquitaine", France. The work of TUL was partly funded by SÚRAO under the project "Research support for the safety assessment of the technical solution of the deep geological repository" (SO2021-053) and under the contract SO2020-017. These supports are gratefully acknowledged. Funding for open access charge: Universidade da Coruña/CISUG.

**Funding** Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. H2020 Euratom (Grant No. 84759), Spanish Ministry of Science and Innovation (PID2023-153202OB-I00), Galician Regional Government (Grant ED431C2021/54), Lithuanian State budget, Carnot ISIFoR Institute

and "la Région Nouvelle-Aquitaine", SÚRAO (contract SO2020-017), The National Natural Science Foundation of China (Grant Nos. No. 42130303, No. 42372283).

**Data availability** The authors confirm that the data supporting the findings of this study are available within the article, its supplementary materials and additional files which are openly available in ZENODO at https://doi.org/10.5281/zenodo.13903270.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
## References

Ahusborde E, El Ossmani M (2017) A sequential approach for numerical simulation of two-phase multicomponent flow with reactive transport in porous media. Math Comput Simul 137:71–89

Ahusborde E, Kern M, Vostrikov V (2015) Numerical simulation of two-phase multicomponent flow with reactive transport in porous media: application to geological sequestration of CO₂. ESAIM Proc 50:21–39

Ahusborde E, Amaziane B, El Ossmani M (2018) Improvement of numerical approximation of coupled two-phase multicomponent flow with reactive geochemical transport in porous media. Oil Gas Sci Technol -Rev, IFP Energies nouvelles 73:73

Ahusborde E, Amaziane B, El Ossmani M, Id Moulay M (2019) Numerical modeling and simulation of fully coupled processes of reactive multiphase flow in porous media. J Math Study 52:359–377

Ahusborde E, Amaziane B, Id Moulay M (2021) High performance computing of 3D reactive multiphase flow in porous media: application to geological storage of CO₂. Comput Geosci 25:2131–2147

Ahusborde E, Amaziane B, Crocololo F, Pillardou N (2023) Numerical simulation of a thermal-hydraulic-chemical multiphase flow model for CO₂ sequestration in saline aquifers. Math Geosci. https://doi.org/10.1007/s11004-023-10093-7

Ahusborde E, Amaziane B, de Hoop S et al (2024) A benchmark study on reactive two-phase flow in porous media: Part II—results and discussion. Comput Geosci. https://doi.org/10.1007/s10596-024-10269-y

Alt-Epping P, Tournassat C, Rasouli P et al (2015) Benchmark reactive transport simulations of a column experiment in compacted bentonite with multispecies diffusion and explicit treatment of electrostatic effects. Comput Geosci 19:535–550

Appelo CAJ, Parkhurst DL, Post VEA (2014) Equations for calculating hydrogeochemical reactions of minerals and gases such as CO₂ at high pressures and temperatures. Geochim Cosmochim Acta 125:49–67

Bastian P, Blatt M, Dedner A, Dreier NA, Engwer C, Fritze R, Gräser C, Grüninger C, Kempf D, Klöfkorn R, Ohlberger M, Sander O (2021) The DUNE framework: basic concepts and recent developments. Comput Math Appl 81:75–112

Bildstein O, Claret F, Frugier P (2019) RTM for waste repositories. Rev Mineral Geochem 85:419–457

Bildstein O, Claret F, Lagneau V (2021) Guest editorial to the special issue: subsurface environmental simulation benchmarks. Comput Geosci 25:1281–1283

Carrayrou J, Kern M, Knabner P (2010) Reactive transport benchmark of MoMaS. Comput Geosci 14:385–392

Chen J, Xu T, Jiang Z, Feng B, Liang X (2020) Reducing formation damage by artificially controlling the fluid-rock chemical interaction in a double-well geothermal heat production system. Renew Energy 149:455–467

Claret F, Marty N, Tournassat C (2018) Modeling the long-term stability of multi-barrier systems for nuclear waste disposal in geological clay formations, reactive transport modeling: applications in subsurface energy and environmental problems. Wiley, Hoboken, pp 395–451

COMSOL (2022) COMSOL multiphysics, version 6.0. www.comsol.com

Dai Z, Samper J (2004) Inverse problem of multicomponent reactive chemical transport in porous media: formulation and applications. Water Resour Res 40:W07407

De Dieuleveult C, Erhel J (2010) A global approach to reactive transport: application to the MoMaS benchmark. Comput Geosci 14:451–464

De Hoop S, Voskov D, Ahusborde E et al (2024) A benchmark study on reactive two-phase flow in porous media: Part I—model description. Comput Geosci 28:175–189

ENRESA (2006) FEBEX: updated final report 1994–2004. ENRESA Technical Publication PT 05-0/2006

Garrels RM, Christ CL (1965) Solutions, minerals and equilibria. Freeman Cooper, San Francisco

Gherardi F, Xu T, Pruess K (2007) Numerical modeling of self-limiting and self-enhancing caprock alteration induced by CO₂ storage in a depleted gas reservoir. Chem Geol 244:103–129

Giffaut E, Grivé M, Blanc P, Vieillard P, Colás E, Gailhanou H, Gaboreau S, Marty N, Madé B, Duro L (2014) ANDRA thermodynamic database for performance assessment: ThermoChimie. Appl Geochem 49:225–236

Helgeson HC, Kirkham DH (1974) Theoretical prediction of the thermodynamic behaviour of aqueous electrolytes at high pressures and temperatures: II. Debye–Hückel parameters for activity coefficients and relative partial molal properties. Am J Sci 274:1199–1261

Hu T, Xu T, Wang F, Yang Y, Wang F, Tian H, Yang Z (2019) The evolution of water chemical characteristics and their indicative function in CO₂-enhanced water recovery. Int J Greenh Gas Control 88:403–415

Hu T, Wang Y, Xu T, Tian H, Rui Z, Zhao Y, Liu T (2023a) A novel technique for quantifying the fate of CO₂ injected in oil reservoir for geological utilization and storage. Energy Rep 9:5350–5361

Hu Y, Gan Q, Hurst A, Elsworth D (2023b) Investigation of coupled hydro-mechanical modelling of hydraulic fracture propagation and interaction with natural fractures. Int J Rock Mech Min Sci 169:105418

Koch T, Gläser D, Weishaupt K, Ackermann S, Beck M, Becker B, Burbulla S, Class H, Coltman E, Emmert S, Fetzer T, Grüninger C, Heck K, Hommel J, Kurz T, Lipp M, Mohammadi F, Scherrer S, Schneider M, Seitz G, Stadler L, Utz M, Weinhardt F, Flemisch B (2021) DuMuₓ 3—an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. Comput Math Appl 81:423–443

Lagneau V, van der Lee J (2010) HYTEC results of the MoMas reactive transport benchmark. Comput Geosci 14:435–449

Mon A (2017) Coupled thermo-hydro-chemical-mechanical models for the bentonite barrier in a radioactive waste repository. Ph.D. Dissertation, University of A Coruña, Spain

Mon A, Samper J, Montenegro L, Turrero MJ, Torres E, Cuevas J, Fernández R, De Windt L (2023) Reactive transport models of the geochemical interactions at the iron/bentonite interface in laboratory corrosion tests. Appl Clay Sci 242:106981

Nardi A, Idiart A, Trinchero P, de Vries LM, Molinero J (2014) Interface COMSOL-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry. Comput Geosci 69:10–21

Nardi A et al (2022) iCP interface COMSOL-PHREEQC. User's guide, v\_2.1, 116 p

Navarro V, Alonso EE (2000) Modeling swelling soils for disposal barriers. Comput Geotech 27:19–43

Parkhurst DL, Appelo CAJ (1999) User's guide to PHREEQC (Version 2). Water Resour Investig Rep 99:4259

Parkhurst DL, Appelo CAJ (2013) Description of input and examples for PHREEQC version 3-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations: U.S. Geological Survey Techniques and Methods, book 6, chap. A43, 497 p

Pollock DW (1986) Simulation of fluid flow and energy transport processes associated with high level radioactive waste disposal in unsaturated alluvium. Water Resour Res 22:765–775
Pruess K, Oldenburf CM, Moridis GJ (1999) TOUGH2 User's Guide Version 2. Office of Scientific & Technical Information Technical Reports LBNL-29400. Berkeley, California

Samper J, Zheng L, Montenegro L, Fernández AM, Rivas P (2008) Testing coupled thermo-hydro-chemical models of compacted bentonite after dismantling the FEBEX in situ test. *Appl Geochem* 23:1186–1201

Samper J, Xu T, Yang C (2009) A sequential partly iterative approach for multicomponent reactive transport with CORE2D. Comput Geosci 13:301–316

Samper J, Mon A, Montenegro L (2018a) A revisited thermal, hydro-dynamic, chemical and mechanical model of compacted bentonite for the entire duration of the FEBEX in situ test. *Appl Clay Sci* 160:58–70

Samper J, Mon A, Montenegro L, Naves A, Cuevas J, Fernández R, Turrero MJ, Torres E (2018b) Coupled THCM model of heating and hydration concrete-bentonite column test. *Appl Geochem* 94:67–81

Samper J, Mon A, Montenegro L (2020a) A coupled THMC model of the geochemical interactions of concrete and bentonite after 13 years of FEBEX plug operation. *Appl Geochem* 121:104687

Samper J, Mon A, Montenegro L (2020b) Reactive transport model of lab tests of lime-mortar, compacted bentonite and magnetite geochemical interactions. *Appl Geochem* 119:104633

Spycher N, Peiffer L, Sonnenthal EL, Saldi G, Reed MH, Kennedy BM (2014) Integrated multicomponent solute geothermometry. *Geothermics* 51:113–123

Tian X, Wapperom M, Gunning J, Jackson S, Wilkins A, Green C, Ennis-King J, Voskov D (2024) A history matching study for the FluidFlower benchmark project. *Transp Porous Media* 151:1113–1139

Villar MV, Martín PL, Bárcena I, García-Siñeriz JL, Gómez-Espina R, Lloret A (2012) Long-term experimental evidences of saturation of compacted bentonite under repository conditions. *Eng Geol* 149:57–69

Wanner C, Peiffer L, Sonnenthal E, Spycher N, Iovenitti J, Kennedy BM (2014) Reactive transport modeling of the Dixie Valley geothermal area: insights on flow and geothermometry. *Geothermics* 51:130–141

Wapperom M, Tian X, Novikov A et al (2024) FluidFlower benchmark: lessons learned from the perspective of subsurface simulation. *Transp Porous Med* 151:1033–1052

Xu T, Pruess K (1998) Coupled modeling of non-isothermal multiphase flow, solute transport and reactive chemistry in porous and fractured media: 1. Model development and validation.

Lawrence Berkeley National Laboratory Report LBNL-42050, Berkeley, California

Xu T, Samper J, Ayora C, Manzano M, Custodio E (1999) Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems. *J Hydrol* 214:144–164

Xu T, Sonnenthal E, Spycher N, Pruess K (2006) TOUGHREACT—a simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media: applications to geothermal injectivity and CO₂ geological sequestration. *Comput Geosci* 32:145–165

Xu T, Senger R, Finsterle S (2008) Corrosion-induced gas generation in a nuclear waste repository: reactive geochemistry and multiphase flow effects. *Appl Geochem* 23:3423–3433

Xu T, Spycher N, Sonnenthal E, Zhang G, Zheng L, Pruess K (2011) TOUGHREACT Version 2.0: a simulator for subsurface reactive transport under non-isothermal multiphase flow conditions. *Comput Geosci* 37:763–774

Zheng L, Samper J (2008) A coupled THMC model of FEBEX mock-up test. *Phys Chem Earth* 33:S486–S498

Zheng L, Samper J (2015) Dual-continuum multicomponent reactive transport with nth-order solute transfer terms for structured porous media. *Comput Geosci* 19:709–726

Zheng L, Samper J, Montenegro L, Mayor JC (2008) Flow and reactive transport model of a ventilation experiment in Opallinus clay. *Phys Chem Earth* 33:S486–S498

Zheng L, Samper J, Montenegro L, Fernández AM (2010) A coupled THMC model of a heating and hydration laboratory experiment in unsaturated compacted FEBEX bentonite. *J Hydrol* 386:80–94

Zheng L, Samper J, Montenegro L (2011) A coupled THC model of the FEBEX in situ test with bentonite swelling and chemical and thermal osmosis. *J Contam Hydrol* 126:45–60

Zhu H, Xu T, Tian H, Feng G, Yang Z, Zhou B (1999) Understanding of long-term CO₂-brine-rock geochemical reactions using numerical modeling and natural analogue study. *Geofluids* 2019:1426061

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.