# Improvement of numerical approximation of coupled multiphase multicomponent flow with reactive geochemical transport in porous media

Etienne Ahusborde, Brahim Amaziane, Mustapha El Ossmani

► To cite this version:

Etienne Ahusborde, Brahim Amaziane, Mustapha El Ossmani. Improvement of numerical approximation of coupled multiphase multicomponent flow with reactive geochemical transport in porous media. Oil & Gas Science and Technology - Revue d'IFP Energies nouvelles, 2018, 73, pp.73. 10.2516/ogst/2018033. hal-01952153

HAL Id: hal-01952153

https://hal.science/hal-01952153

Submitted on 11 Dec 2018

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
Numerical methods and HPC

A. Anciaux-Sedrakian and Q. H. Tran (Guest editors)

REGULAR ARTICLE

OPEN ACCESS

# Improvement of numerical approximation of coupled multiphase multicomponent flow with reactive geochemical transport in porous media

Etienne Ahusborde<fnref n="1" />, Brahim Amaziane<fnref n="1" />,\* and Mustapha El Ossmani<fnref n="1" />,<fnref n="2" />

Received: 20 February 2018 / Accepted: 6 July 2018

**Abstract.** In this paper, we consider a parallel finite volume algorithm for modeling complex processes in porous media that include multiphase flow and geochemical interactions. Coupled flow and reactive transport phenomena often occur in a wide range of subsurface systems such as hydrocarbon reservoir production, groundwater management, carbon dioxide sequestration, nuclear waste repository or geothermal energy production. This work aims to develop and implement a parallel code coupling approach for non-isothermal multiphase multicomponent flow and reactive transport simulation in the framework of the parallel open-source platform DuMu<sup>X</sup>. Modeling such problems leads to a highly nonlinear coupled system of degenerate partial differential equations to algebraic or ordinary differential equations requiring special numerical treatment. We propose a sequential fully implicit scheme solving firstly a multiphase compositional flow problem and then a Direct Substitution Approach (DSA) is used to solve the reactive transport problem. Both subsystems are discretized by a fully implicit cell-centred finite volume scheme and then an efficient sequential coupling has been implemented in DuMu<sup>X</sup>. We focus on the stability and robustness of the coupling process and the numerical benefits of the DSA approach. Parallelization is carried out using the DUNE parallel library package based on MPI providing high parallel efficiency and allowing simulations with several tens of millions of degrees of freedom to be carried out, ideal for large-scale field applications involving multicomponent chemistry. As we deal with complex codes, we have tested and demonstrated the correctness of the implemented software by benchmarking, including the MoMaS reactive transport benchmark, and comparison to existing simulations in the literature. The accuracy and effectiveness of the approach is demonstrated through 2D and 3D numerical simulations. Parallel scalability is investigated for 3D simulations with different grid resolutions. Numerical results for long-term fate of injected CO<sub>2</sub> for geological storage are presented. The numerical results have demonstrated that this approach yields physically realistic flow fields in highly heterogeneous media and showed that this approach performs significantly better than the Sequential Iterative Approach (SIA).

## 1 Introduction

Multiphase multicomponent flow with chemical reactions in porous media play a significant role for many applications in geological and reservoir engineering processes. We can mention for instance the sequestration of CO₂ in saline aquifers, the geological storage of nuclear waste or the prevention of groundwater pollution and the contaminant remediation. Numerical models have been increasingly used for this purpose, a trend that will continue because more sophisticated models and codes are being developed and computer costs keep decreasing. Significant efforts and

attempts have been made during recent years toward the development of such tools. There is an extensive literature on this subject. We will not attempt a literature review here but will merely mention a few references. We refer, for instance, to the books [1, 2] and the references therein. A recent review of reactive codes for subsurface environmental simulation can be viewed in [3].

Carbon Capture and Storage (CCS) is a promising way to mitigate the effects of global warming. Assessing the viability of geological storage must rely on numerical simulations due to the long time scales involved. Several physical and geochemical trapping mechanisms must be combined to ensure a high containment rate, and geochemical trapping becomes increasingly important over longer

<footnote n="1">¹ CNRS / Univ Pau & Pays Adour / E2S UPPA, Laboratoire de Mathématiques et de leurs applications de PAU, Fédération IPRA, UMR 5142, 64000 Pau, France</footnote>

<footnote n="2">² University Moulay Ismaïl, EMMACS-ENSAM, Marjane II, 50000 Meknès, Morocco</footnote>

<footnote>\* Corresponding author: brahim.amaziane@univ-pau.fr</footnote>
time scales [4]: carbon dissolution in water occurs over hundreds of years, and formation of carbonate minerals over millions of years, see for instance [5]. Many references can be found for the numerical approximation of such phenomena, see e.g., [6, 7], and the references therein. Two main issues concern the reactive transport in the framework of geological storage of CO₂: firstly the simulation of the different trapping mechanisms and secondly the sustainability of the storage through the durability of well cements.

For the numerical simulation of the geological sequestration, the studies can be classified into two categories according on whether we consider a single (liquid) or two-phase flow (supercritical/gas-liquid). Under certain simplifying conditions, single-phase flow is considered for instance in [8] and [9]. In [8], in the framework of SHPCO₂ Project [10], the gas phase is assumed to be immobile and therefore gaseous carbon dioxide is considered as a fixed species neglecting the two-phase flow effects. In [9], an initial amount of supercritical CO₂ is converted into a source term of liquid CO₂ and then the authors study the transport of the dissolved CO₂ and the precipitation/dissolution process of minerals. In [11], the authors employ a single-phase reactive flow to model the leaking of CO₂-saturated brine in a fractured pathway once supercritical CO₂ is totally dissolved. CO₂ is generally injected in its supercritical form. This injection may induce important pressure build-up that can damage the reservoir or induce fracturing and seismic events. Moreover, the supercritical CO₂ that is less dense than the brine present in the aquifer, will migrate vertically firstly and then along the top of the aquifer and finally it builds up under the cover rock inducing a risk of leakage through faults. In [12], the authors propose an alternative strategy that consists in injecting dissolved CO₂ to circumvent the above-mentioned risks and increase the security of its geological sequestration. In [13], a study of this process and its interactions with the carbonate reservoir through geochemical reactions is proposed. Even if in [14], the authors show that modified one-phase flow models can predict pressure build-up far from the injection as well as complex two-phase flow models, most of the studies deal with two-phase reactive flows. The following is a non-exhaustive list of references: [15–21], the main difference between these references being the complexity of the geochemical system.

During its storage, liquid or supercritical CO₂ is injected through a well composed of cases of cement. Once the injection is terminated, the well is closed by a cement plug. By consequent, the integrity of the disposal and its impermeability can be strongly influenced by the chemical reactivity of cement (see [22] or [23] where a comparison between numerical and experimental results is given).

The long-term safety of the disposal of nuclear waste is an important issue in all countries with a significant nuclear program. Repositories for the disposal of high-level and long-lived radioactive waste generally rely on a multi-barrier system to isolate the waste from the biosphere. The multi-barrier system typically comprises the natural geological barrier provided by the repository host rock and its surroundings, and an engineered barrier system. As for the geological storage of CO₂, the integrity and the sustainability

of the disposal must be ensured for thousands of years. In [24–28], the authors study the corrosion or/and the alteration of carbon steel, compacted bentonite and concrete that constitute the engineered multi-barrier system. In the event of a leak of radionuclide, the chemical interactions between the leaked radionuclide and the surrounding media must be taken into account to predict at best their migration (see for instance [29–32]). Many references related to the Hanford site in the USA can be found (see for instance [33, 34]). In this case, cesium has been released into the environment accidentally due to leaking of High Level Waste (HLW) storage tanks and the retardation effect arising from adsorption is studied.

Groundwater contamination is a crucial issue in the management of the water quality. Reactive transport modeling and simulations allow to better understand the physical, chemical and biological processes involved in contaminant transport and remediation, see for instance [35–37].

Multiphase multicomponent reactive flows are modeled by a mass balance law for each phase, Darcy-Muskat's law, capillary pressure law, solubility laws and equations of state. Coupling between flow and chemistry occurs through reactions rates. In the case of equilibrium reactions, these rates are unknown and are commonly eliminated through linear transformations [38, 39] and replaced by mass actions laws that are algebraic equations relating the activities of concerned species. For kinetic reactions, the rates are nonlinear functions of concentrations [40] and involve ordinary differential equations. By consequence, the problem is modeled by a system of partial differential equations (describing a multiphase compositional flow) coupled with algebraic or ordinary differential equations related to chemical reactions.

The numerical methods for solving this system can be divided into three dominant algorithms: the Global Implicit (GIA), the Sequential Iterative (SIA) and Sequential Non-Iterative (SNIA) approaches [41, 42]. In the GIA, the nonlinear system gathering all equations is solved at each time step. For the sequential solution approaches, flow and reactive transport (or possibly, flow, transport and chemistry) are solved sequentially at each time step. The difference between the SIA and SNIA lies on the fact that for the SIA, the procedure is present in an iterative loop. Sequential approaches are also named operator-splitting approaches.

No comparative study exists yet to quantify the accuracy loss for sequential approaches in comparison with GIA, but their gain in implementation and saving in computing time is fairly obvious. Nonetheless, sequential approaches introduce operator splitting errors [43, 44] and restrictions on the time step are mandatory to ensure mass conservation for instance. In [42], the authors describe the GIA as “research tools for one-dimensional investigations” due to their complexity and their high computational requirements. Thanks to the advance of high-performance computing in the last decades, these restrictions are no longer relevant. The Groupement Mathematical Modeling and Numerical Simulation for Nuclear Waste Management has proposed in [45] a benchmark to test numerical methods
used to deal with reactive transport problem in porous media. In this framework several sequential and implicit algorithms have been compared. In [46] and [47], the authors propose respectively a SNIA and a SIA. The other participants [48–51] deal with various global implicit algorithms. More precisely, in [48], the authors propose a method where the chemical problem is eliminated locally, leading to a nonlinear system where the transport and chemistry subsystems remain separated. In [49, 52], the problem is written in the form of Differential Algebraic Equations (DAE). In [50], the authors use a reduction technique introduced in [53, 54] that aims to reduce the number of coupled nonlinear differential equations drastically. Finally in [51], a Direct Substitution Approach (DSA) consisting in substituting the equations of chemistry directly in the equations of transport is employed. In [55], the results provided by the different teams are compared with a good agreement.

This work deals with a numerical method to perform numerical simulation of two-phase multicomponent flow with reactive transport in porous media. The major difficulties related to this model are in the nonlinear degenerate structure of the equations, as well as in the strong coupling between the flow and reactive transport equations. The chemical processes involve inter-phase mass transfer as well as a host of chemical reactions, including dissolution, ion exchange, adsorption, precipitation, and oxidation/ reduction. The aim of this paper is to extend our previous results by using a GIA for solving the reactive transport problem. We develop an efficient parallel algorithm implementing the discrete model based on a finite volume scheme. Parallel computing offers an opportunity for building detailed models for such problems and in providing the capability of solving larger, more realistic and practical problems faster.

As in [56, 57], we propose a sequential approach that splits the original problem into two sub-problems. The first sub-problem computes an implicit two-phase compositional flow where only species present in both phases are taken into account. Exchanges between phases are totally solved in this step and the contribution of the other species is treated explicitly. The second sub-problem calculates a reactive transport problem where flow properties (Darcy velocity for each phase, saturation of each phase, temperature, density, etc.) are given by the first step. In [56, 57], a SIA has been implemented for the reactive transport sub-problem. To improve the robustness of the scheme and the accuracy loss due to the time-splitting involved by the SIA, in this paper, we switched to GIA for the reactive transport subsystem. More precisely, we used a DSA [42].

The rest of the paper is organized as follows. In Section 2, we describe the governing equations for a two-phase multicomponent flow with reactive transport. Our methodology consisting in decoupling the original problem is explained in Section 3 where the finite-volume discretization of each sub-problem is also detailed. In Section 4, a description of the implementation of our strategy in the DuMuX framework [58, 59] is given and some numerical results are presented. The developed model is effectively verified by published results. We present two examples

involving reactive transport modeling to demonstrate capabilities of the code. The first one aims to validate the geochemistry module. The second one reports the results of the numerical investigations of long-term fate of injected CO₂ for geological sequestration including parallel scaling results. For this test, a comparison between the DSA and the SIA for solving the reactive transport sub-problem is given and proves that the DSA is sensibly less CPU consuming than the SIA. In Section 5, some concluding remarks are presented.

## 2 Governing equations

In this section, the governing equations that describe two-phase multicomponent flow in a porous medium coupled with chemical reactions are presented.

### 2.1 Geochemical model

Before describing the chemical reactions and distinguishing those at equilibrium from the kinetic one, we specify some useful notations to lighten the writing of equations. Indeed, we note by $I$ the set of all the chemical components involved in the chemical reactions. One distinguishes among these components the primary and secondary ones that are noted respectively $I_p$ and $I_s$ such that $I = I_p \cup I_s$. We note by $I_{se}$ the set of components of $I_s$ that are in equilibrium reactions and by $I_{sk}$ those in kinetic reactions ($I_s = I_{se} \cup I_{sk}$).

In order to use a sequential approach to decouple two-phase flow and reactive transport, we introduce two subsets $I_{2\phi}$ and $I_{rt}$. $I_{2\phi}$ contains all the species present in both phases while $I_{rt}$ consists of the remaining ones such that $I = I_{2\phi} \cup I_{rt}$. So $N_c = \text{card}\{I\}$ is the total number of chemical components involved in $N_r = \text{card}\{I_s\}$ chemical reactions ($N_r = N_e + N_k$), with $N_e = \text{card}\{I_{se}\}$ is the number of reactions at equilibrium and $N_k = \text{card}\{I_{sk}\}$ is the number of kinetic reactions. The chemical system writes:

$$ \sum_{j=1}^{N_c} v_{ij}A_j = 0, \quad i = 1, \dots, N_r, $$

where $v_{ij}$ is the stoichiometric coefficient of the species $A_j$ in the reaction $i$.

As an illustration example, we present in Table 1 a chemical system used in [16] in a scenario of CO₂ injection in a deep saline aquifer.

Ano, Cal and Kao refer respectively to Anorthite, Calcite and Kaolinite. In this case, $N_c = 14$ and $N_r = 8$. The first five reactions are at equilibrium while the last three ones are kinetic. Below we give the various corresponding sets above-mentioned:

$$ \begin{aligned} I_s &= [\text{H}_2\text{O}_{(g)} \text{CO}_{2(g)} \text{HCO}_3^- \text{CO}_3^{2-} \text{OH}^- \text{Ano Cal Kao}], \\ I_p &= [\text{H}_2\text{O}_{(l)} \text{CO}_{2(l)} \text{H}^+ \text{Ca}^{2+} \text{Al}^{3+} \text{SiO}_{2(l)}], \\ I_{2\phi} &= [\text{H}_2\text{O}_{(l)} \text{CO}_{2(l)} \text{H}_2\text{O}_{(g)} \text{CO}_{2(g)}], \\ I_{rt} &= [\text{H}^+ \text{Ca}^{2+} \text{Al}^{3+} \text{SiO}_{2(l)} \text{HCO}_3^- \text{CO}_3^{2-} \text{OH}^- \text{Ano Cal Kao}], \\ I_{se} &= [\text{H}_2\text{O}_{(g)} \text{CO}_{2(g)} \text{HCO}_3^- \text{CO}_3^{2-} \text{OH}^-], \\ I_{sk} &= [\text{Ano Cal Kao}]. \end{aligned} $$
<table><caption>**Table 1.** Chemical reactions.</caption><thead><tr><th colspan="2">Reactions</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sub>(g)</sub></td><td>= H<sub>2</sub>O<sub>(l)</sub></td></tr><tr><td>CO<sub>2</sub>(g)</td><td>= CO<sub>2</sub>(l)</td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>= CO<sub>2</sub>(l) + H<sub>2</sub>O<sub>(l)</sub> - H<sup>+</sup></td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>= CO<sub>2</sub>(l) + H<sub>2</sub>O<sub>(l)</sub> - 2H<sup>+</sup></td></tr><tr><td>OH<sup>-</sup></td><td>= H<sub>2</sub>O<sub>(l)</sub> - H<sup>+</sup></td></tr><tr><td>Ano</td><td>= 4H<sub>2</sub>O<sub>(l)</sub> + Ca<sup>2+</sup> + 2Al<sup>3+</sup> + 2SiO<sub>2</sub>(l) - 8H<sup>+</sup></td></tr><tr><td>Cal</td><td>= Ca<sup>2+</sup> + HCO<sub>3</sub><sup>-</sup> - H<sup>+</sup></td></tr><tr><td>Kao</td><td>= 5H<sub>2</sub>O<sub>(l)</sub> + 2Al<sup>3+</sup> + 2SiO<sub>2</sub>(l) - 6H<sup>+</sup></td></tr></tbody></table>

Each equilibrium reaction gives rise to an algebraic relation called mass action law that links the activities of the species involved in the reaction. The mass action law writes as follows:

$$a_{\alpha}^{j} = K_{j} \prod_{i \in I_{p}} (a_{\alpha}^{i})^{v_{ji}}, \quad j \in I_{sc}, \qquad (1)$$

where $a_{\alpha}^j$ is the activity of species $j$ in its phase $\alpha$, $K_j$ is the equilibrium constant of reaction $j$. The activity of water and solid species are set equal to 1. The activity of aqueous species is often written in terms of molality:

$a_1^j = \gamma_1^j \frac{m_1^j}{m_0^j}$, where $\gamma_1^j$ is the activity coefficient for species $j$ in the aqueous phase, $m_1^j$ is the molality of species $j$ [mol kg$^{-1}$] in the aqueous phase, and $m_0$ is standard molality often taken as 1 mol kg$^{-1}$. For each aqueous species $j$, the molality $m_1^j$ and the mole fraction of species $j$ in aqueous phase $x_1^j$ are related by $m_1^j = \frac{x_1^j}{M_{H_2O} x_1^{H_2O}}$, where

$M_{H_2O}$ is the molar mass of water and $x_1^{H_2O}$ is the molar fraction of water in the aqueous phase. The activity of gaseous species is often represented by fugacity:

Each kinetic reaction leads to an ordinary differential equation:

$$\frac{dc_s^j}{dt} = -r_i, \quad i \in I_{sk}, \qquad (2)$$

where $c_s^j$ denotes the concentration of mineral $i$ and $r_i$ is the kinetic reaction rate depending on the activities of the species present in the reaction. Expression for $r_i$ can be found for instance in [40].

## 2.2 Mathematical model for two-phase multicomponent flow with reactive transport

In the sequel, the index $\alpha \in \{l, g, s\}$ (l for liquid, g for gas and s for solid) refers to the phase, while the superscript $i$ refers to the component. We consider the mass balance

equation (see for instance [60]) for primary and secondary species in its phase:

$$\begin{aligned} \frac{\partial}{\partial t}(\theta_{\alpha} \rho_{\text{mol}, \alpha} x_{\alpha}^{i}) + \nabla \cdot (\rho_{\text{mol}, \alpha} x_{\alpha}^{i} \vec{q}_{\alpha}) - \nabla \cdot (\rho_{\text{mol}, \alpha} D_{\alpha} \nabla x_{\alpha}^{i}) \\ = \sum_{j \in I_s} v_{ji} r_j, \quad i \in I_p, \end{aligned} \qquad (3)$$

$$\begin{aligned} \frac{\partial}{\partial t}(\theta_{\alpha} \rho_{\text{mol}, \alpha} x_{\alpha}^{i}) + \nabla \cdot (\rho_{\text{mol}, \alpha} x_{\alpha}^{i} \vec{q}_{\alpha}) - \nabla \cdot (\rho_{\text{mol}, \alpha} D_{\alpha} \nabla x_{\alpha}^{i}) \\ = -r_{i}, \quad i \in I_{s}, \end{aligned} \qquad (4)$$

where $\theta_{\alpha}$ [-] denotes the volumetric content of phase $\alpha$ ($\theta_{\alpha} = \phi S_{\alpha}$, $\phi$ [-] being the porosity of the medium and $S_{\alpha}$ [-] the saturation of phase $\alpha$ if $\alpha \in \{l, g\}$ and $\theta_{s} = 1 - \phi$), $\rho_{\text{mol}, \alpha}$ [mol m$^{-3}$] is the molar density of phase $\alpha$, $x_{\alpha}^i$ is the molar fraction of species $i$ in phase $\alpha$, $D_{\alpha}$ [m$^2$ s$^{-1}$] denotes the diffusivity of the phase $\alpha$:

$$D_{\alpha} = \theta_{\alpha} \left( D_{\text{m}} \| + d_{\text{L}} |\vec{q}_{\alpha}| \| + (d_{\text{L}} - d_{\text{T}}) \frac{\vec{q}_{\alpha} \vec{q}_{\alpha}^{T}}{|\vec{q}_{\alpha}|} \right).$$

For the sake of simplicity, the diffusion coefficient is assumed to be independent of the chemical species $i$ and depends only on the phase $\alpha$. $D_{\text{m}}$ [m$^2$ s$^{-1}$] is the molecular diffusion, $d_{\text{L}}$ [m] and $d_{\text{T}}$ [m] are the magnitudes of longitudinal and transverse dispersion respectively and $\vec{q}_{\alpha}$ [m s$^{-1}$] is the Darcy velocity of phase $\alpha$, expressed as follows:

$$\vec{q}_{\alpha} = - \frac{k_{r\alpha}(S_l)}{\mu_{\alpha}} \mathbb{K}(\nabla P_{\alpha} - \rho_{\alpha} \vec{g}), \qquad (5)$$

where $k_{r\alpha}$ [-] denotes the relative permeability of phase $\alpha$, $\mu_{\alpha}$ [Pa s] is the dynamic viscosity of phase $\alpha$, $\mathbb{K}$ [m$^2$] is the absolute permeability tensor, $P_{\alpha}$ [Pa] is the pressure of phase $\alpha$, $\rho_{\alpha}$ [kg m$^{-3}$] is the mass density of phase $\alpha$ and $\vec{g}$ [m s$^{-2}$] is the gravitational acceleration.

The phase pressures are connected by the capillary pressure law:

$$P_c(S_l) = P_g - P_l, \qquad (6)$$

Finally, $r_j$ [mol m$^{-3}$ s$^{-1}$] is the rate of reaction $j$ (it can be equilibrium or kinetic).

We choose to reformulate equations (3) and (4) in term of molar concentration $c_{\alpha}^i = \rho_{\text{mol}, \alpha} x_{\alpha}^i$. By neglecting the gradient of molar density $\nabla \rho_{\text{mol}, \alpha}$ resulting from the diffusive flux, equations (3) and (4) write:

$$\frac{\partial}{\partial t}(\theta_{\alpha} c_{\alpha}^{i}) + \nabla \cdot (c_{\alpha}^{i} \vec{q}_{\alpha}) - \nabla \cdot (D_{\alpha} \nabla c_{\alpha}^{i}) = \sum_{j \in I_s} v_{ji} r_j, \quad i \in I_p, \qquad (7)$$

$$\frac{\partial}{\partial t}(\theta_{\alpha} c_{\alpha}^{i}) + \nabla \cdot (c_{\alpha}^{i} \vec{q}_{\alpha}) - \nabla \cdot (D_{\alpha} \nabla c_{\alpha}^{i}) = -r_{i}, \quad i \in I_{s}. \qquad (8)$$

For sake of simplicity, we introduce the advection-diffusion operator:

$$L_{\alpha}(c) = \nabla \cdot (c \vec{q}_{\alpha}) - \nabla \cdot (D_{\alpha} \nabla c). \qquad (9)$$
We can remark that, since $D_s = 0$ and $\vec{q}_s = \vec{0}$, so that $L_s(\cdot) = 0$.

In order to eliminate the reaction rates $r_j$ in equations (7), we make linear combinations between equations (8) with each equation (7). This introduces $N_p$ new conservation laws that write:

$$ \frac{\partial}{\partial t} \left( \theta_{\alpha} c_{\alpha}^{j} + \sum_{j \in I_s} \theta_{\alpha} v_{ji} c_{\alpha}^{j} \right) + L_{\alpha}(c_{\alpha}^{j}) + \sum_{j \in I_s} L_{\alpha}(v_{ji} c_{\alpha}^{j}) = 0, \quad i \in I_p. \quad (10) $$

To retrieve the same number of equations as there are unknowns, the $N_s$ equations (8) are replaced by $N_e$ mass actions laws defined by (1) corresponding to the equilibrium reactions and $N_k$ ordinary differential equations corresponding to the kinetic reactions given by (2).

## 3 Numerical formulation

As in [56, 57], we adopt a sequential strategy that consists in simplifying the original problem (3)–(4) by a splitting into two sub-problems. The first sub-problem is dedicated to a two-phase compositional flow with only the species belonging to liquid and gas phases and we assume that all the exchanges between the phases are computed in this step. The second sub-problem is a reactive transport problem that computes the contribution of the other species.

### 3.1 Formulation of the problem

In this subsection, we present the formulation of a geochemical model and the governing equations modeling two-phase multicomponent flow with reactive transport in porous media that will be used in the sequel.

The system of equations for the two-phase compositional flow is given by:

$$ \frac{\partial}{\partial t} \left( \theta_l c_l^i + \theta_g c_g^i \right) + L_l(c_l^i) + L_g(c_g^i) + \sum_{j \in I_s / I_{2\phi}} v_{ji} \left[ \frac{\partial}{\partial t} (\theta_\alpha c_\alpha^j) + L_\alpha(c_\alpha^j) \right] = 0, \quad i \in I_p \cap I_{2\phi}, \quad (11) $$

$$ a_g^i = K_i a_i^i, \quad \bar{i} \in I_{\text{se}} \cap I_{2\phi}. \quad (12) $$

Equation (11) corresponds to equation (10) for the species present in liquid and gas phases where $i$ and $\bar{i}$ represent respectively the species in the liquid and the gas phase that are in chemical equilibrium and the phase $\alpha$ denotes the phase containing the species $j$. Equation (12) is the solubility law where $a_i^i$ and $a_g^i$ are the activities of the species in liquid and gas phases. It traduces the phase equilibrium.

The system of equations describing the reactive transport problem writes:

$$ \frac{\partial}{\partial t} \left( \theta_{\alpha} c_{\alpha}^{i} + \sum_{j \in I_s} \theta_{\alpha} v_{ji} c_{\alpha}^{j} \right) + L_{\alpha}(c_{\alpha}^{i}) + \sum_{j \in I_s} L_{\alpha}(v_{ji} c_{\alpha}^{j}) = 0, \quad i \in I_p \cap I_{rt}, \quad (13) $$

$$ a_{\alpha}^{j} = K_{j} \prod_{i \in I_{p}} (a_{\alpha}^{i})^{v_{ji}}, \quad j \in I_{\text{se}} \setminus I_{2\phi}, \quad (14) $$

$$ \frac{\mathrm{d}c_s^i}{\mathrm{d}t} = -K_i^s A_i^s \left( 1 - K_i \prod_{j \in I_p} (a_j^i)^{v_{ji}} \right), \quad i \in I_{sk}, \quad (15) $$

where $K_i^s$ and $A_i^s$ are respectively the kinetic-rate constant $[\text{mol m}^{-2} \text{s}^{-1}]$ and the reactive surface area $[\text{m}^2 \text{m}^{-3}]$ of mineral $i$.

Equation (13) corresponds to equation (10) for the primary species present only in one phase whereas equation (14) is the mass action law involved by equilibrium reactions. Finally, equation (15) corresponds to the ordinary differential equations involved by kinetic reactions.

Subsystem (11)–(12) and (13)–(15) are solved sequentially. In equation (13), the Darcy velocities, the phase saturations, the molar densities and the activities or species $i \in I_{2\phi}$ are provided after the computation of the subsystem (11)–(12). In this latter, subsystem (13)–(15) provides explicitly the contribution of species $j \in I_s / I_{2\phi}$ as well as a possible update of porosity thanks to the computation of the mineral concentrations.

### 3.2 Finite volume discretization of two-phase compositional flow and reactive transport

The spatial discretization of the coupled system (11)–(15) employs a conservative Finite Volume (FV) method based on a fully upwinding scheme to treat the convective terms and a conforming finite element scheme with piecewise linear elements for the diffusive terms. The time discretization is done by an implicit Euler method.

Here, we choose a cell-centred FV method. It consists in integrating the equations (11)–(15) on a control volume $V_k$ (see Fig. 1) and evaluating the fluxes at the interface $\gamma_{kl}$ between two adjacent elements $V_k$ and $V_l$.

– Spatial discretization of two-phase compositional flow:

$$ \int_{V_k} \frac{\partial}{\partial t} \left\{ \theta_l c_l^i + \theta_g c_g^i \right\} dV + \sum_{l \in V(k)} \int_{\gamma_{kl}} \left\{ c_l^i \vec{q}_l + c_g^i \vec{q}_g \right\} \cdot \vec{n}_{kl} d\gamma \\ - \sum_{l \in V(k)} \int_{\gamma_{kl}} \left\{ D_l \nabla c_l^i + D_g \nabla c_g^i \right\} \cdot \vec{n}_{kl} d\gamma \\ = - \sum_{j \in I_s / I_{2\phi}} v_{ji} \left( \int_{V_k} \frac{\partial}{\partial t} (\theta_\alpha c_\alpha^j) dV \right. \\ \left. + \sum_{l \in V(k)} \int_{\gamma_{kl}} \left\{ c_\alpha^j \vec{q}_\alpha - D_\alpha \nabla c_\alpha^j \right\} \cdot \vec{n}_{kl} d\gamma \right), \quad i \in I_p \cap I_{2\phi}, \quad (16) $$

$$ \{a^i\}_k = K_i \{a^i\}_k, \quad \bar{i} \in I_{\text{se}} \cap I_{2\phi}, \quad (17) $$

$$ \vec{q}_{\alpha} = -K \lambda_{\alpha}(S_l)(\nabla P_{\alpha} - \rho_{\alpha} \vec{g}), \quad (18) $$

where $\lambda_{\alpha}(S_l) = \frac{k_{\alpha}(S_l)}{\mu_{\alpha}}$ is the phase mobility,
<figure><img src="image_2.png" /><figcaption>Fig. 1. Discretization by the cell-centered finite volume method.</figcaption></figure>

- Spatial discretization of reactive transport:

$$ 
\begin{aligned}
& \int_{V_k} \frac{\partial}{\partial t} \left\{ \theta_\alpha c_\alpha^i + \sum_{j \in I_s} v_{ji} \theta_\alpha c_\alpha^j \right\} dV + \\
& \sum_{l \in V(k)} \int_{\gamma_{kl}} \left\{ c_\alpha^i \vec{q}_\alpha + \sum_{j \in I_s} v_{ji} c_\alpha^j \vec{q}_\alpha \right\} \cdot \vec{n}_{kl} d\gamma \\
& - \sum_{l \in V(k)} \int_{\gamma_{kl}} \left\{ D_\alpha \nabla c_\alpha^i + \sum_{j \in I_s} v_{ji} D_\alpha \nabla c_\alpha^j \right\} \cdot \vec{n}_{kl} d\gamma \\
& = 0, \quad i \in I_p \cap I_{rt},
\end{aligned}
\quad (19) $$

$$ \{a_{\alpha}^{j}\}_{k} = K_{j} \prod_{i \in I_{p}} \{(a_{\alpha}^{i})^{v_{ji}}\}_{k}, j \in I_{se}/I_{2\phi}, \quad (20) $$

$$ \frac{d\{c_s^i\}_k}{dt} + K_i^s A_i^s (1 - K_i \prod_{j \in I_p} \{(a_l^j)^{v_{ji}}\}_k) = 0, \quad i \in I_{sk}, \quad (21) $$

where $\vec{n}_{kl}$ denotes the unit outer normal to $\gamma_{kl}$, $V(k)$ is the set of adjacent elements of $V_k$ and $\{\cdot\}_k$ represents the average values on $V_k$. By using the implicit Euler scheme for the time discretization and due to the fact that the primary unknowns ($S_\alpha$, $P_\alpha$, $c_\alpha^i$) and the physical parameters are constant on each element $V_k$, the cell-centred FV schemes corresponding to the discretization of the flow and reactive transport are given by:

- Cell-centred FV scheme of two-phase compositional flow:

$$ 
\begin{aligned}
& \frac{|V_k|}{\Delta t^m} \left( \left\{ \theta_\alpha c_1^i + \theta_\alpha c_2^j \right\}_k^{n+1} - \left\{ \theta_\alpha c_1^i + \theta_\alpha c_2^j \right\}_k^n \right) \\
& + \sum_{l \in V(k)} |\gamma_{kl}| \left( \left\{ c_1^i \right\}_{kl}^{n+1} \left\{ \vec{q}_1 \right\}_{kl}^{n+1} + \left\{ c_2^j \right\}_{kl}^{n+1} \left\{ \vec{q}_2 \right\}_{kl}^{n+1} \right) \cdot \vec{n}_{kl} \\
& - \sum_{l \in V(k)} |\gamma_{kl}| \left( \left\{ D_1 \right\}_{kl}^{n+1} \left\{ \nabla c_1^i \right\}_{kl}^{n+1} + \left\{ D_2 \right\}_{kl}^{n+1} \left\{ \nabla c_2^j \right\}_{kl}^{n+1} \right) \cdot \vec{n}_{kl} \\
& = - \sum_{j \in I_s / I_{2\phi}} v_{ji} \frac{|V_k|}{\Delta t^m} \left( \left\{ \theta_\alpha c_\alpha^j \right\}_k^n - \left\{ \theta_\alpha c_\alpha^i \right\}_k^{n-1} \right) \\
& - \sum_{j \in I_s / I_{2\phi}} v_{ji} \sum_{l \in V(k)} |\gamma_{kl}| \left\{ c_\alpha^j \right\}_{kl}^n \left\{ \vec{q}_\alpha \right\}_{kl}^{n+1} \cdot \vec{n}_{kl} \\
& + \sum_{j \in I_s / I_{2\phi}} v_{ji} \sum_{l \in V(k)} |\gamma_{kl}| \left\{ D_\alpha \right\}_{kl}^{n+1} \left\{ \nabla c_\alpha^j \right\}_{kl}^n \cdot \vec{n}_{kl}, \quad i \in I_p \cap I_{2\phi},
\end{aligned}
\quad (22) $$

$$ \{a_g^{\bar{i}}\}_k^{n+1} = K_i \{a_i^j\}_k^{n+1}, \quad \bar{i} \in I_{se} \cap I_{2\phi}, \quad (23) $$

$$ \{\vec{q}_\alpha\}_{kl}^{n+1} = -\{\mathbb{K}\}_{kl} \{\lambda_\alpha(S_l)\}_{kl}^{n+1} \left( \{\nabla P_\alpha\}_{kl}^{n+1} - \{\rho_\alpha\}_{kl}^{n+1} \vec{g} \right), \quad (24) $$

- Cell-centred FV scheme for reactive transport:

$$ 
\begin{aligned}
& \frac{|V_k|}{\Delta t^m} \left( \left\{ \theta_\alpha c_\alpha^i + \sum_{j \in I_s} v_{ji} \theta_\alpha c_\alpha^j \right\}_k^{n+1} - \left\{ \theta_\alpha c_\alpha^i + \sum_{j \in I_s} v_{ji} \theta_\alpha c_\alpha^j \right\}_k^n \right) \\
& + \sum_{l \in V(k)} |\gamma_{kl}| \left( \left\{ c_\alpha^i \right\}_{kl}^{n+1} \left\{ \vec{q}_\alpha \right\}_{kl}^{n+1} + \sum_{j \in I_s} v_{ji} \left\{ c_\alpha^j \right\}_{kl}^{n+1} \left\{ \vec{q}_\alpha \right\}_{kl}^{n+1} \right) \cdot \vec{n}_{kl} \\
& - \sum_{l \in V(k)} |\gamma_{kl}| \left( \left\{ D_\alpha \right\}_{kl}^{n+1} \left\{ \nabla c_\alpha^i \right\}_{kl}^{n+1} \right. \\
& \qquad \left. + \sum_{j \in I_s} v_{ji} \left\{ D_\alpha \right\}_{kl}^{n+1} \left\{ \nabla c_\alpha^j \right\}_{kl}^{n+1} \right) \cdot \vec{n}_{kl} = 0, \quad i \in I_p \cap I_{rt},
\end{aligned}
\quad (25) $$

$$ \{a_{\alpha}^{j}\}_{k}^{n+1} = K_{j} \prod_{i \in I_{p}} \{(a_{\alpha}^{i})^{v_{ji}}\}_{k}^{n+1}, \quad j \in I_{se}/I_{2\phi}, \quad (26) $$

$$ \{c_s^i\}_k^{n+1} = \{c_s^i\}_k^n - \Delta t^m K_i^s A_i^s \left( 1 - K_i \prod_{j \in I_p} \{(a_l^j)^{v_{ij}}\}_k^{n+1} \right), \quad i \in I_{sk}. \quad (27) $$

Now to define the finite volume scheme it is enough to approach the convective and the diffusive fluxes on the interfaces $\gamma_{kl}$. For this purpose, a fully upwinding scheme is used to calculate the numerical flux for the convective term. More precisely, the quantities ($S_\alpha$, $c_\alpha^i$, $P_\alpha$, $\lambda_\alpha$) are evaluated implicitly and upstream at the interface $\gamma_{kl}$ between two adjacent elements as:

$$ \{\cdot\}_{kl}^{n+1} = \begin{cases} \{\cdot\}_{k}^{n+1} & \text{if } \{\vec{q}_\alpha\}_{kl}^{n+1} \cdot \vec{n}_{kl} > 0, \\ \{\cdot\}_{l}^{n+1} & \text{else.} \end{cases} \quad (28) $$

The gradient operators on the interfaces $\gamma_{kl}$ are calculated by a $\mathbb{P}_1/\mathbb{Q}_1$ finite element method with piecewise linear elements. An harmonic average of the values between two adjacent elements is used to calculate the absolute permeability $\{\mathbb{K}\}_{kl}$ and the diffusion coefficients $\{D_\alpha\}_{kl}^{n+1}$ at the interface $\gamma_{kl}$. $\{\rho_\alpha\}_{kl}^{n+1}$ is computed as the arithmetic average of two elements $V_k$ and $V_l$.
# 4 Numerical simulations

## 4.1 Implementation

The physical models described above are built in DuMu<sup>X</sup> [58, 59] framework, which handles general input/output, memory management, grid generation, parallelism, etc. The code is an object-oriented software written in C++ and equipped with efficient solvers and massively parallel computation capability.

The first sub-problem consists in solving a convection-diffusion system modeling a two-phase compositional flow. To tackle this system, we have used a module implemented in DuMu<sup>X</sup> called 2pnc (two-phase, n-component), where n denotes the number of components present in both phases. The approach is fully implicit. The spatial discretization is performed by the cell-centred finite volume approach described in Section (3.2) by equations (22)–(24). The nonlinear system is solved by a Newton method and a preconditioned BiConjugate Gradient STABilized (BiCGSTAB) method is used to solve the linear system. The control of the time-step is based on the number of iterations required by the Newton method to achieve convergence for the last time iteration. The time-step is reduced, if the number of iterations exceeds a specified threshold, whereas it is increased if the method converges within less iterations. The contribution of the other species is treated explicitly and put as a source term in right hand side of the discretized mass conservation laws (22). We have also implemented the change of porosity linked to the concentrations of minerals.

The second sub-problem consists in solving a single-phase reactive transport problem. In [57], we have developed in the DuMu<sup>X</sup> framework a one-phase multicomponent transport module named 1pmc-react. In this context, we have implemented a new 1pmc module (one-phase, m-component) that is coupled with the calculations of the chemical problem computed with GSL [61] through a SIA. In this step, the properties of the flow (velocity, saturation, density, concentration of species present in both phases) are given by the two-phase compositional flow. Here, we propose a GIA for this reactive transport sub-problem. More precisely, we used a DSA [42]. The DSA consists in integrating directly the mass action laws (26) in the discretized conservation laws (25). Again, the spatial discretization is carried out using the finite volume method described in Section (3.2) by equations (25)–(27).

A parallel simulator has been developed based on these approximation techniques and is being used to study the sequestration of CO<sub>2</sub> and gas migration in a nuclear waste repository. Non-isothermal reactive transport simulations based on 1D, 2D and 3D unstructured grids can be performed. Geochemical reactions can be treated as equilibrium or kinetic. The entire framework preserves its flexibility for further numerical developments.

## 4.2 Numerical results

To validate our strategy, several tests have been performed. In [62], we studied the migration of hydrogen produced by the corrosion reaction in deep geological

radioactive waste repository and obtained some results coinciding with those presented in [27]. In this contribution, we focus on two tests cases. The first one aims to validate the implementation of the B-dot model [63] for the activity, chemical equilibrium calculations, mineral dissolution and precipitation reaction and the chosen numerical scheme to deal with geochemistry. The second test [16] deals with a scenario of injection of CO<sub>2</sub> in a deep saline aquifer where parallel computations are performed and their performances are detailed.

### 4.2.1 Validation of geochemistry calculations

We consider a test case presented in [18] which aims to validate the implementation of the activity, chemical equilibrium calculations, mineral dissolution and precipitation reaction and the chosen numerical scheme to deal with geochemistry. There is no flow, no injection and no production. In [18], the authors compared their results with those obtained by the geochemistry software “The Geochemist’s Workbench® (GWB)” [64]. GWB models only the geochemistry but not the transport of components in the porous medium. Consequently, as in [18], we used a one dimensional grid for an example considering only chemistry with no space dependence. The test deals with the dissolution of the calcite at 25 °C with the reactions presented in Table 2. Reactions (1)–(6) are in equilibrium while (7) is a kinetic reaction whose parameters are displayed in Table 3.

Initial molalities [mol kg<sup>-1</sup>] are exhibited in Table 4. The activity of each aqueous species is written in terms of molality as:

$$d_i^i = \gamma^i m_i^i, \qquad (29)$$

where γ<sup>i</sup> is the activity coefficient for species i in the aqueous phase. Several models of activity coefficient can be used. Here, the B-dot model is considered [63]:

$$\gamma^i = - \frac{Az_i \sqrt{I}}{1 + \dot{A}_i B \sqrt{I}} + \dot{B}I, \qquad (30)$$

where A, B and Ḃ are constants depending on temperature, z<sub>i</sub> is the ion electrical charge for species i and Ȧ<sub>i</sub> is the ion size of species i. Finally, I is the ionic strength of the aqueous phase and expresses as I = 0.5 Σ<sub>i</sub> m<sub>i</sub>z<sub>i</sub><sup>2</sup>.

Tables 4 and 5 exhibit the parameters for the B-dot model. Initially, the system contains 10 g of calcite and 100 kg of H<sub>2</sub>O.

Figures 2 and 3 show respectively the evolution of the molality of the different components and the quantity of dissolved calcite as a function of the time. Solid lines represent our results and markers are used for results obtained in [18]. We can observe a very good agreement between both calculations.

### 4.2.2 Application to geological CO₂ sequestration

We now apply our simulator to a 3D coupled two-phase flow and reactive transport problem to compare our results with those of [16]. In this work the authors propose to use
<table><caption>**Table 2.** Chemical reactions.</caption><thead><tr><th>No.</th><th>Reactions</th><th>log<sub>10</sub> (K)</th></tr></thead><tbody><tr><td>(1)</td><td>H<sub>2</sub>O = OH<sup>-</sup> + H<sup>+</sup></td><td>-13.9868</td></tr><tr><td>(2)</td><td>CO<sub>2(l)</sub> + H<sub>2</sub>O = H<sup>+</sup> + HCO<sub>3</sub><sup>-</sup></td><td>-6.366</td></tr><tr><td>(3)</td><td>HCO<sub>3</sub><sup>-</sup> = CO<sub>3</sub><sup>2-</sup> + H<sup>+</sup></td><td>-10.3439</td></tr><tr><td>(4)</td><td>Ca<sup>2+</sup> + HCO<sub>3</sub><sup>-</sup> = CaHCO<sub>3</sub><sup>+</sup></td><td>1.22</td></tr><tr><td>(5)</td><td>Ca<sup>2+</sup> + Cl<sup>-</sup> = CaCl<sup>+</sup></td><td>0.7</td></tr><tr><td>(6)</td><td>HCO<sub>3</sub><sup>-</sup> + Ca<sup>2+</sup> = CaCO<sub>3(l)</sub> + H<sup>+</sup></td><td>-7.128</td></tr><tr><td>(7)</td><td>Calcite + H<sup>+</sup> = Ca<sup>2+</sup> + HCO<sub>3</sub><sup>-</sup></td><td>1.7130</td></tr></tbody></table>

<table><caption>**Table 3.** Mineral, precipitation and dissolution parameters.</caption><thead><tr><th>Mineral</th><th>log<sub>10</sub> (K<sup>s</sup>)<br>[mol m<sup>-2</sup> s<sup>-1</sup>]</th><th>A<sup>s</sup><br>[m<sup>2</sup> m<sup>-3</sup>]</th><th>log<sub>10</sub> (K)</th><th>Init.<br>[mol]</th></tr></thead><tbody><tr><td>Calcite</td><td>-8.80</td><td>8.8</td><td>1.7130</td><td>0.1</td></tr></tbody></table>

<table><caption>**Table 4.** Input parameters for each ion.</caption><thead><tr><th>Ion</th><th>Init. molality<br>[mol kg<sup>-1</sup>]</th><th>Charge z</th><th>Size (Å)</th></tr></thead><tbody><tr><td>H<sup>+</sup></td><td>9.787 × 10<sup>-7</sup></td><td>1</td><td>9</td></tr><tr><td>CO<sub>2(l)</sub></td><td>4.393 × 10<sup>-4</sup></td><td>0</td><td>3</td></tr><tr><td>Ca<sup>2+</sup></td><td>3.723 × 10<sup>-4</sup></td><td>2</td><td>6</td></tr><tr><td>Cl<sup>-</sup></td><td>5.73 × 10<sup>-4</sup></td><td>3</td><td>9</td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>2.102 × 10<sup>-4</sup></td><td>-1</td><td>4.5</td></tr><tr><td>CaHCO<sub>3</sub><sup>-</sup></td><td>1.126 × 10<sup>-6</sup></td><td>1</td><td>4</td></tr><tr><td>CaCl<sup>+</sup></td><td>8.649 × 10<sup>-7</sup></td><td>1</td><td>4.1</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>4.23 × 10<sup>-7</sup></td><td>-2</td><td>4.5</td></tr><tr><td>CaCO<sub>3(l)</sub></td><td>5.131 × 10<sup>-9</sup></td><td>0</td><td>3.4</td></tr><tr><td>OH<sup>-</sup></td><td>1.133 × 10<sup>-8</sup></td><td>-1</td><td>3.5</td></tr></tbody></table>

two variants of complex geochemical systems that include both equilibrium and kinetic reactions. We consider the test named “six-element model”, whose reactions are displayed in Table 6. It involves four equilibrium reactions (the first four reactions) and three kinetic reactions (the last three ones) of mineral dissolution/precipitation. For the first reaction, the solubility law for CO₂ is implemented according to [65]. Mineral data for the kinetic reactions are summarized in Table 7.

A three-dimensional domain that is 15 km in both the *x* and *y*-directions and 100 m in the *z*-direction is considered. A well perforated in a single grid block located 25 m from the top of the aquifer injects a pure CO₂ stream at constant rate during the first 20 years. After the 20 years injection period, a total of 18.6 × 10⁹ kg of CO₂ is injected.

<table><caption>**Table 5.** B-dot model parameters from the EQ3/6 database [66].</caption><thead><tr><th>T (°C)</th><th>25</th><th>60</th><th>100</th><th>150</th><th>200</th></tr></thead><tbody><tr><td>A</td><td>0.4939</td><td>0.5465</td><td>0.5995</td><td>0.6855</td><td>0.7994</td></tr><tr><td>B</td><td>0.3253</td><td>0.3346</td><td>0.3421</td><td>0.3525</td><td>0.3639</td></tr><tr><td>˙B</td><td>0.041</td><td>0.0438</td><td>0.046</td><td>0.047</td><td>0.047</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Evolution of the molality of the different components versus time.</figcaption></figure>

As initial conditions for the two-phase two-component H₂O-CO₂ flow we have used hydrostatic condition for liquid pressure P₁, initial liquid saturation S₁ = 1 and initial CO₂ molality in liquid phase equals 3.55 × 10⁻³ mol kg⁻¹. Initial conditions for the reactive transport problem and parameters for the B-dot model used as activity model are shown in Table 8.

Impermeable Neumann boundary conditions are enforced on the boundaries of the domain. Constitutive laws and physical parameters are given in Table 9.

Figure 4 displays concentrations of calcite, anorthite and kaolinite at 20 and 2000 years for a grid composed of 1.6 × 10⁵ elements (100 × 100 × 16). Initially, their concentration were respectively 238, 87 and 88 mol m⁻³. We can see that calcite is dissolved near the injection of CO₂ and precipitated far from the injection while anorthite and
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the number of moles of dissolved calcite.</figcaption></figure>

<table><caption>**Table 6.** Chemical reactions.</caption><thead><tr><th>Reactions</th><th>log<sub>10</sub> (K)</th></tr></thead><tbody><tr><td>CO<sub>2(g)</sub> = CO<sub>2(l)</sub></td><td>-</td></tr><tr><td>CO<sub>2(l)</sub> + H<sub>2</sub>O = H<sup>+</sup> + HCO<sub>3</sub><sup>-</sup></td><td>-10.23</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup> + H<sup>+</sup> = HCO<sub>3</sub><sup>-</sup></td><td>-6.32</td></tr><tr><td>OH<sup>-</sup> + H<sup>+</sup> = H<sub>2</sub>O</td><td>-13.26</td></tr><tr><td>Ano + 8H<sup>+</sup> = 4H<sub>2</sub>O + Ca<sup>2+</sup> + 2Al<sup>3+</sup> + 2SiO<sub>2(l)</sub></td><td>25.82</td></tr><tr><td>Cal + H<sup>+</sup> = Ca<sup>2+</sup> + HCO<sub>3</sub><sup>-</sup></td><td>1.6</td></tr><tr><td>Kao + 6H<sup>+</sup> = 5H<sub>2</sub>O + 2Al<sup>3+</sup> + 2SiO<sub>2(l)</sub></td><td>6.82</td></tr></tbody></table>

<table><caption>**Table 7.** Mineral, precipitation and dissolution parameters.</caption><thead><tr><th>Mineral</th><th>log<sub>10</sub> (K<sup>s</sup>)</th><th>A<sup>s</sup></th><th>Init. conc.</th></tr></thead><tbody><tr><td>Ano</td><td>-12.0</td><td>88</td><td>87</td></tr><tr><td>Cal</td><td>-8.80</td><td>88</td><td>238</td></tr><tr><td>Kao</td><td>-13.0</td><td>17 600</td><td>88</td></tr></tbody></table>

kaolinite are respectively dissolved and precipitated everywhere.

Figure 5 depicts the molality of aqueous CO₂, the pH and the gas saturation at 20 and 2000 years. The pH and the molality of aqueous CO₂ are strongly correlated since high concentration of aqueous CO₂ acidifies the liquid phase. After the injection, the gaseous CO₂ migrates upward and spreads laterally when reaching the top of the aquifer that is impermeable.

Figure 6 exhibits the mineral evolution and carbon distribution over time. Negative value corresponds to dissolution while positive value means precipitation. We can see that globally, calcite and kaolinite precipitate while

<table><caption>**Table 8.** Input parameters for each ion.</caption><thead><tr><th>Ion</th><th>Charge z</th><th>Size (Å)</th><th>Init. mol.</th></tr></thead><tbody><tr><td>CO<sub>2(l)</sub></td><td>0</td><td>3</td><td>3.55 × 10<sup>-3</sup></td></tr><tr><td>H<sup>+</sup></td><td>1</td><td>9</td><td>5.71 × 10<sup>-7</sup></td></tr><tr><td>Al<sup>3+</sup></td><td>3</td><td>9</td><td>3.13 × 10<sup>-12</sup></td></tr><tr><td>SiO<sub>2(l)</sub></td><td>0</td><td>3</td><td>4.73 × 10<sup>-4</sup></td></tr><tr><td>Ca<sup>2+</sup></td><td>2</td><td>6</td><td>2.52 × 10<sup>-2</sup></td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>-1</td><td>4.5</td><td>2.14 × 10<sup>-3</sup></td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>-2</td><td>4.5</td><td>4.23 × 10<sup>-7</sup></td></tr><tr><td>OH<sup>-</sup></td><td>-1</td><td>3.5</td><td>4.49 × 10<sup>-8</sup></td></tr><tr><td>Al(OH)<sub>2</sub><sup>+</sup></td><td>1</td><td>5</td><td>1.28 × 10<sup>-9</sup></td></tr></tbody></table>

<table><caption>**Table 9.** Physical parameters for the test case of CO₂ injection.</caption><thead><tr><th>Constitutive law</th><th>Parameters</th></tr></thead><tbody><tr><td>Capillary pressure law</td><td>P<sub>c</sub> = 0 Pa</td></tr><tr><td>Absolute permeability</td><td>K = 10<sup>-13</sup> [m<sup>2</sup>]</td></tr><tr><td>Relative permeability</td><td></td></tr><tr><td>k<sub>rl</sub> = (S<sub>l</sub><sup>*</sup>)<sup>4</sup></td><td>S<sub>l</sub><sup>*</sup> = S<sub>l</sub> - S<sub>lr</sub></td></tr><tr><td>k<sub>rg</sub> = 0.4(1 - S<sub>l</sub><sup>*</sup>)<sup>2</sup>(1 - (S<sub>l</sub><sup>*</sup>)<sup>2</sup>)</td><td>S<sub>lr</sub> = 0.2</td></tr><tr><td>Liquid diffusion tensor</td><td></td></tr><tr><td>D<sub>l</sub> = D<sub>m</sub></td><td>D<sub>m</sub> = 10<sup>-9</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Gas diffusion tensor</td><td>Model based on [67]</td></tr><tr><td>Porosity</td><td>φ = 0.18</td></tr><tr><td>Temperature</td><td>T = 50 °C</td></tr><tr><td>Liquid density</td><td>Model based on [68]</td></tr><tr><td>Liquid viscosity</td><td>μ<sub>l</sub> = 4.8 10<sup>-4</sup> Pa s<sup>-1</sup></td></tr><tr><td>Gas density</td><td>Model based on [69]</td></tr><tr><td>Gas viscosity</td><td>Model based on [70]</td></tr></tbody></table>

anorthite is dissolved. The distribution of CO₂ shows that most of CO₂ is present in gas phase but a non negligible part of carbon is stored in mineral form. It is emphasized in Table 10 that summarizes the carbon distribution at t = 2000 years. Numerical results are comparable with those in [16].

Parallel performance

Parallel computations up to 768 processors have been performed on several grids. The parallel efficiency of our strategy is illustrated by solving 10 time steps. The code ran on a Bull cluster named OCCIGEN and Intel "Haswell" 12-Core E5-2690 V3 processors. In parallel computing, two types of scalability are defined. The first is the strong scaling, which represents the relation between the computation time and the number of processors for a fixed total problem size. The second is the weak scaling, for which the load per processor is fixed.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Profiles of mineral concentrations.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Profiles of pH, aqueous CO₂ molality and gas saturation.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Mineral net molar changes and evolution of CO₂ in time.</figcaption></figure>

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

### Strong scaling

Figure 7a displays on a logarithmic scale, CPU time as a function of the number of processors for three size problems of $1.6 \times 10^5$, $1.28 \times 10^6$ and $5.76 \times 10^6$ elements corresponding respectively to $1.92 \times 10^6$, $1.536 \times 10^7$ and $6.912 \times 10^7$ degrees of freedom (dof). The dashed lines represent an ideal behaviour.

Strong efficiency is given by:

$$SE(N) = \frac{\text{CPU time on } p \text{ processors} \times p}{\text{CPU time on } N \text{ processors} \times N}, \quad (31)$$

where *p* denotes the number of processors used for the reference time (not always equal to one for heavy computations). For this calculation, we took *p* = 12, 24 and 48. It points out an optimal use of the parallel resources. Efficiency equal to one indicates that communications and synchronizations between processors are negligible.

Figure 7b represents the strong scaling versus the number of processors. A high efficiency (greater than 0.70) is observed up to 256 processors for the computations involving $1.536 \times 10^7$ and $6.912 \times 10^7$ dof. For the simulation with $1.92 \times 10^6$ dof, the efficiency is good up to 72 processors. The loss of efficiency is mainly due to the increase of the communications between processors in comparison with the load of each processor. In [71, 72], the authors evaluate the parallel performance of the simulators PFLOTRAN. In [71], the authors assert that “as a general rule of thumb a minimum of 10,000 dof per core is needed to obtain good scaling performance”. In [73], for the simulator ParCrunchFlow, strong scaling breaks down somewhere between 69 000 and 40 000 dof per processor. Here, a minimum of 30 000 dof per processor seems to be required to maintain a good scaling.

### Weak scaling

Figure 8a displays CPU time as a function of the number of processors, with around 10 000, 20 000 and 40 000 elements per processor.

<table><caption>Table 10. Final carbon distribution at 2000 years.</caption><thead><tr><th>In CO<sub>2(g)</sub></th><th>In CO<sub>2(l)</sub></th><th>In ions</th><th>In minerals</th></tr></thead><tbody><tr><td>69.2%</td><td>20.9%</td><td>0.9%</td><td>9.0%</td></tr></tbody></table>

Weak efficiency is given by:

$$WE(N) = \frac{\text{CPU time on } p \text{ processors}}{\text{CPU time on } N \text{ processors}}, \quad (32)$$

where *p* still denotes the number of processors used for the reference time. Here, *p* = 1 for the three scenarios. Efficiency equal to one indicates an optimal behaviour for the algorithm and the computer architecture. Indeed, CPU time remains constant, equal to the reference time, while the total size of the problem increases with the number of processors. Usually, this property is hardly verified and curves with plateaus can be observed. Values of the plateaus rise toward one with the load of each processor. This phenomenon is illustrated in Figure 8b.

### 4.2.3 Comparison between direct substitution and sequential iterative approaches

This subsection aims at comparing the DSA used in this work with the SIA developed in [57] for solving the reactive transport sub-problem for the example presented above. Both approaches adopt an adaptive time-stepping. In the DSA, the control of the time-step is based on the number of iterations required by the Newton method to achieve convergence while in the SIA, it is based on the number of iterations required in the iterative algorithm to reach the tolerance $\epsilon_{SIA}$. In the sequel, tolerances for the Newton method and iterative algorithm are respectively $\epsilon_{Newton} = 10^{-8}$ and $\epsilon_{SIA} = 10^{-8}$.

Figure 9 compares the evolution of the molalities of H⁺ and Ca²⁺ obtained with the DSA and the SIA close to the injection with a mesh composed of 40 000 cells during
<figure><img src="image_9.png" /><figcaption>Fig. 7. CPU time and strong parallel efficiency as a function of the number of processors.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. CPU time and weak parallel efficiency as a function of the number of processors.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Comparison of the molalities of H⁺ and Ca²⁺ obtained with the DSA and the SIA.</figcaption></figure>
<table><caption>**Table 11.** CPU time (s) for the DSA and the SIA.</caption><thead><tr><td>Number of cells</td><td>DSA</td><td>SIA</td></tr></thead><tbody><tr><td>10 000</td><td>982.2</td><td>1518.8</td></tr><tr><td>40 000</td><td>3847.5</td><td>6235.5</td></tr><tr><td>160 000</td><td>15901.9</td><td>27235.5</td></tr></tbody></table>

the first year of injection. We can observe that the results are in good accordance and that both methods provide comparable results.

Table 11 displays the CPU time required for the DSA and the SIA to achieve the first year of injection years on several meshes. We can see that for this example, the DSA is faster than the SIA whatever the size of the mesh. The SIA is more CPU consuming in comparison with the DSA because many iteration steps in the iterative procedure and smaller time steps are required.

## 5 Conclusion

A mathematical formulation and a finite volume approximation were presented for a system of coupled partial differential and algebraic or ordinary differential equations describing two-phase flow with reactive geochemical transport in the subsurface. The proposed reactive transport software module coupling to the two-phase flow was implemented in the framework of the parallel open-source platform DuMu<sup>X</sup>. Numerical results are in good agreement with those in [16, 18] and provide additional validation of the method. Moreover, this simulation study demonstrates the ability of the model and the method to accurately capture flow behaviour in 3D reservoirs. Due to the strong coupling of the flow and reactive transport equations, a standard approach is to use a fully implicit approach to ensure stability in the solution. Although this guarantees numerical stability in the solution, this does not guarantee nonlinear convergence. This complexity makes analysing the entire nonlinear problem difficult. A separation of the different physics could improve the understanding and result in a better design of nonlinear solvers for the reactive transport problem. The analyses on the sequential scheme were only performed for numerical tests. Mathematical analysis of the different coupling schemes is ongoing research. The study still needs to be improved by developing a fully implicit scheme and compare it to the sequential approaches. We note finally that we encountered difficulties to find reliable and well documented benchmarks. In many articles, some data are missing. We think that a well documented benchmark for two-phase flow with reactive transport in porous media would be very useful for the community. Further work on these important issues is in progress.

**Acknowledgments.** This work was partially supported by the Carnot Institute, ISIFOR project (Institute for the Sustainable Engineering of Fossil Resources) and CDAPP (Communauté d'Agglomération de Pau Pyrénées). Their supports are gratefully acknowledged. We also thank the DuMu<sup>X</sup> and DUNE

teams for their help during the development of our reactive transport module and the CINES (National Computing Center for Higher Education) to give us access to their computing resources facility. This work was granted access to the HPC resources of CINES under the allocations 2017-A0010610019 and 2018-A0040610019 made by GENCI. The authors gratefully thank the anonymous referees for their insightful comments and suggestions.

## References

1. Niemi A., Bear J., Bensabat J. (2017) *Geological storage of CO₂ in deep saline formations*, Springer.

2. Zhang F., Yeh G.T., Parker J.C. (2012) *Groundwater reactive transport models*, Bentham e-books

3. Steefel C.I., Appelo C.A.J., Arora B., Jacques D., Kalbacher T., Kolditz O., Lagneau V., Lichtner P.C., Mayer K.U., Meeussen J.C.L., Molins S., Moulton D., Shao H., Šimunek J., Spycher N., Yabusaki S.B., Yeh G.T. (2015) Reactive transport codes for subsurface environmental simulation, *Comput. Geosci.* **19**, 445–478.

4. Jiang X. (2011) A review of physical modelling and numerical simulation of long-term geological storage of CO₂, *Appl. Energy* **88**, 3557–3566.

5. Intergovernmental Panel on Climate Change (IPCC). (2005) IPCC special report on carbon dioxide capture and storage, in: Metz B., Davidson O., de Coninck H.C., Loos M., Loos M., Meyer L.A. (eds.), *IPCC special report on carbon dioxide capture and storage*, Cambridge University Press. Prepared by Working Group III of the Intergovernmental Panel on Climate Change.

6. Al-Khoury R., Bundschuh J. (2014) *Computational models for CO₂ geo-sequestration and compressed air energy storage*, Sustainable Energy Developments, CRC Press.

7. Bear J., Carrera J. (2017) *Mathematical modeling of CO₂ storage in a geological formation*, Springer.

8. Haeberlein F. (2011) Time space domain decomposition methods for reactive transport – application to CO₂ geological storage, *PhD Thesis*, Université Paris-Nord – Paris XIII.

9. Lagneau V., Pipart A., Catalette H. (2005) Reactive transport modelling of CO₂ sequestration in deep saline aquifers, *Oil Gas Sci. Technol.* – *Rev. IFP Energies nouvelles* **60**, 231–247.

10. Haeberlein F., Michel A., Halpern L. (2009) A test case for multi-species reactive-transport in heterogeneous porous media applied to CO₂ geological storage. https://www.ljll.math.upmc.fr/mcparis09/Files/haeberlein_poster.pdf.

11. Ahmad N., Wörman A., Bottacin-Busolin A., Sanchez-Vila X. (2015) Reactive transport modeling of leaking CO₂-saturated brine along a fractured pathway, *Int. J. Greenh. Gas Con.* **42**, 672–689.

12. Pool M., Carrera J., Vilarrasa V., Silva O., Ayora C. (2013) Dynamics and design of systems for geological storage of dissolved CO₂, *Adv. Water Resour.* **62**, 533–542.

13. Ahmad N., Wörman A., Sanchez-Vila X., Jarsjö J., Bottacin-Busolin A., Hellevang H. (2016) Injection of CO₂ saturated brine in geological reservoir: A way to enhanced storage safety, *Int. J. Greenh. Gas Con.* **54**, 129–144.

14. Nicot J.P., Hosseini S.A., Solano S.V. (2011) Are single-phase flow numerical models sufficient to estimate pressure distribution in CO₂ sequestration projects?, *Energy Procedia* **4**, 3919–3926.
15 Audigane P., Gaus I., Czernichowski-Lauriol I., Pruess K., Xu T. (2007) Two-dimensional reactive transport modelling of CO₂ injection in a saline aquifer at the Sleipner site, North Sea, *Am. J. Sci.* **307**, 974–1008.

16 Fan Y., Durlofsky L.J., Tchelepi H.A. (2012) A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO₂ storage simulations, *Adv. Water Resour.* **42**, 47–61.

17 Leal A.M.M., Blunt M.J., LaForce T.C. (2013) A robust and efficient numerical method for multiphase equilibrium calculations: Application to CO₂-brine-rock systems at high temperatures, pressures and salinities, *Adv. Water Resour.* **62**, 409–430.

18 Nghiem L., Sammon P., Grabenstetter J., Ohkuma H. (2004) Modeling CO₂ storage in aquifers with a fully-coupled geochemical eos compositional simulator, *SPE - DOE Improved Oil Recovery Symposium Proceedings*.

19 Nghiem L., Shrivastava V., Kohse B. (2011) Modeling aqueous phase behavior and chemical reactions in compositional simulation, *Society of Petroleum Engineers - SPE Reservoir Simulation Symposium* **1**, 454–468.

20 Saaltink M., Vilarrasa V., De Gaspari F., Silva O., Carrera J., Rötting T.S. (2013) A method for incorporating equilibrium chemical reactions into multiphase flow models for CO₂ storage, *Adv. Water Resour.* **62**, 431–441.

21 Thibeau S., Nghiem L.X., Ohkuma H. (2007) A modeling study of the role of selected minerals in enhancing CO₂ mineralization during CO₂ aquifer storage, *Proceedings - SPE Annual Technical Conference and Exhibition* **2**, 906–922.

22 Huet B.M., Prevost J.H., Scherer G.W. (2010) Quantitative reactive transport modeling of portland cement in CO₂-saturated water, *Int. J. Greenh. Gas Con.* **4**, 561–574.

23 Jacquemet N., Pironon J., Lagneau V., Saint-Marc J. (2012) Armouring of well cement in H₂S-CO₂ saturated brine by calcite coating - experiments and numerical modelling, *Appl. Geochem.* **27**, 782–795.

24 Berner U., Kulik D.A., Kosakowski G. (2013) Geochemical impact of a low-pH cement liner on the near field of a repository for spent fuel and high-level radioactive waste, *Phys. Chem. Earth* **64**, 46–56.

25 Mon A., Samper J., Montenegro L., Naves A., Fernández J. (2017) Long-term non-isothermal reactive transport model of compacted bentonite, concrete and corrosion products in a hlw repository in clay, *J. Contam. Hydrol.* **197**, 1–16.

26 Sedighi M., Thomas H.R., Al Masum S., Vardon P.J., Nicholson D., Chen Q. (2015) Geochemical modelling of hydrogen gas migration in an unsaturated bentonite buffer, *Geol. Soc. Spec. Publ.* **415**, 189–201.

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

28 Xu T., Senger R., Finsterle S. (2011) Bentonite alteration due to thermal-hydro-chemical processes during the early thermal period in a nuclear waste repository, *Nucl. Technol.* **174**, 438–451.

29 Shao H., Dmytrieva S.V., Kolditz O., Kulik D.A., Pfingsten W., Kosakowski G. (2009) Modeling reactive transport in non-ideal aqueous-solid solution system, *Appl. Geochem.* **24**, 1287–1300.

30 Spycher N.F., Sonnenthal E.L., Apps J.A. (2003) Fluid flow and reactive transport around potential nuclear waste emplacement tunnels at yucca mountain, Nevada, *J. Contam. Hydrol.* **62**, 653–673.

31 Viswanathan H.S., Robinson B.A., Valocchi A.J., Triay I.R. (1998) A reactive transport model of neptunium migration from the potential repository at yucca mountain, *J. Hydrol.* **209**, 251–280.

32 De Windt L., Pellegrini D., Van Der Lee J. (2004) Coupled modeling of cement/claystone interactions and radionuclide migration, *J. Contam. Hydrol.* **68**, 165–182.

33 Lichtner P.C., Yabusaki S., Pruess K., Steefel C.I. (2004) Role of competitive cation exchange on chromatographic displacement of cesium in the vadose zone beneath the hanford s/sx tank farm, *Vadose Zone Journal* **3**, 203–219.

34 Steefel C.I., Carroll S., Zhao P., Roberts S. (2003) Cesium migration in hanford sediment: A multi-site cation exchange model based on laboratory transport experiments, *J. Contam. Hydrol.* **67**, 219–246.

35 Appelo C.A.J., Postma D. (2005) *Geochemistry, Groundwater and Pollution*, 2nd edn., Taylor & Francis.

36 Bear J., Cheng A.H.-D. (2010) *Modeling groundwater flow and contaminant transport*, Springer.

37 Zheng C., Bennett G.D. (2002) *Applied contaminant transport modeling*, John Wiley and Sons, New York.

38 Lichtner P.C. (1985) Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems, *Geochim. Cosmochim. Acta* **49**, 779–800.

39 Molins S., Carrera J., Ayora C., Saaltink M.W. (2004) A formulation for decoupling components in reactive transport problems, *Water Resour. Res.* **40**, 1–13.

40 Lasaga A.C., Soler J.M., Ganor J., Burch T.E., Nagy K.L. (1994) Chemical weathering rate laws and global geochemical cycles, *Geochim. Cosmochim. Acta* **58**, 2361–2386.

41 Steefel C.I., MacQuarrie K.T.B. (1996) Approaches to modeling of reactive transport in porous media, *Rev. Mineral.* **34**, 82–129.

42 Yeh G.T., Tripathi V.S. (1991) A model for simulating transport of reactive multispecies components: Model development and demonstration, *Water Resour. Res.* **27**, 3075–3094.

43 Barry D.A., Miller C.T., Culligan-Hensley P.J. (1996) Temporal discretisation errors in non-iterative split-operator approaches to solving chemical reaction/groundwater transport models, *J. Contam. Hydrol.* **22**, 1–17.

44 Valocchi A.J., Malmstead M. (1992) Accuracy of operator splitting for advection-dispersion-reaction problems, *Water Resour. Res.* **28**, 1471–1476.

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

46 Carrayrou J. (2010) Looking for some reference solutions for the reactive transport benchmark of MoMaS with SPECY, *Comput. Geosci.* **14**, 393–403.

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

48 Amir L., Kern M. (2010) A global method for coupling transport with chemistry in heterogeneous porous media, *Comput. Geosci.* **14**, 465–481.

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

50 Hoffmann J., Kräutle S., Knabner P. (2010) A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the MoMaS benchmark problem, *Comput. Geosci.* **14**, 421–433.
51 Mayer K.U., MacQuarrie K.T.B. (2010) Solution of the MoMaS reactive transport benchmark with MIN3P-model formulation and simulation results, *Comput. Geosci.* **14**, 405–419.

52 Erhel J., Sabit S. (2017) Analysis of a global reactive transport model and results for the MoMaS benchmark, *Math. Comput. Simulat.* **137**, 286–298.

53 Kräutle S., Knabner P. (2005) A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media, *Water Resour. Res.* **41**, 1–17.

54 Kräutle S., Knabner P. (2007) A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions, *Water Resour. Res.* **43**.

55 Carrayrou J., Hoffmann J., Knabner P., Kräutle S., de Dieuleveult C., Erhel J., Van der Lee J., Lagneau V., Mayer K.U., MacQuarrie K.T.B. (2010) Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems – the MoMaS benchmark case, *Comput. Geosci.* **14**, 483–502.

56 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. Surveys* **50**, 21–39.

57 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. Simulat.* **137**, 71–89.

58 DuMu$^X$ (2018) *DUNE for multi-Phase, Component, Scale, Physics, ... flow and transport in porous media*, https://www.dumux.org, last accessed February 1, 2018.

59 Flemisch B., Darcis M., Erbertseder K., Faigle B., Lauser A., Mosthaf K., Muthing S., Nuske P., Tatomir A., Wolf M., Helmig R. (2011) DuMu$^X$: DUNE for multi-Phase, Multicomponent, Scale, Physics, ...} flow and transport in porous media, *Adv. Water Resour.* **34**, 1102–1112.

60 Helmig R. (1997) *Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems*, Springer.

61 GSL - GNU Scientific LibraryMultidimensional Root-Finding. https://www.gnu.org/software/gsl/, Last accessed February 1, 2018.

62 Ahusborde E., Amaziane B., El Ossmani M. (2017) Finite volume scheme for coupling two-phase flow with reactive transport in porous media, *Finite volumes for complex applications VIII-hyperbolic, elliptic and parabolic problems*,

*Springer Proceedings in Mathematics and Statistics* **200**, 407–415.

63 Kirkham D.H., Helgeson H.C., Flowers G.C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures: IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600 °C and 5 KB, *Am. J. Sci.* **281**, 1249–1516.

64 Bethke C., Farrell B., Yeakel S., Yeakel S. (2018) *The Geochemist's Workbench® Release 12 - GWB Essentials Guide*. https://www.gwb.com/pdf/GWB12/GWBessentials.pdf.

65 Spycher N., Pruess K. (2005) CO₂-H₂O mixtures in the geological sequestration of CO₂. II. Partitioning in chloride brines at 12–100 °C and up to 600 bar, *Geochim. et Cosmochim. Acta* **69**, 3309–3320.

66 Wolery T.J. (1992) *EQ3/6 software package for geochemical modeling of aqueous systems: Package overview and installation guide (version 8.0)*, Lawrence Livermore National Laboratory Report UCRL-MA-110662 PT I.

67 Xu B., Nagashima K., DeSimone J.M., Johnson C.S. (2003) Diffusion of water in liquid and supercritical carbon dioxide: an NMR study, *J. Phys. Chem. A* **107**, 1–3.

68 Adams J.J., Bachu S. (2002) Equations of state for basin geofluids: algorithm review and intercomparison for brines, *Geofluids* **2**, 257–271.

69 Span R., Wagner W. (1996) A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa, *J. Phys. Chem. Ref. Data* **25**, 1–88.

70 Fenghour A., Wakeham W.A., Vesovic V. (1998) The viscosity of carbon dioxide, *J. Phys. Chem. Ref. Data* **27**, 31–44.

71 Hammond G.E., Lichtner P.C., Lu C., Mills R.T. (2012) PFLOTRAN: Reactive flow & transport code for use on laptops to leadership-class supercomputers, *Groundwater Reactive Transport Models* 141–159.

72 Hammond G.E., Lichtner P.C., Mills R.T. (2014) Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN, *Water Resour. Res.* **50**, 208–228.

73 Beisman J.J., Maxwell R.M., Navarre-Sitchler A.K., Steefel C.I., Molins S. (2015) ParCrunchFlow: an efficient, parallel reactive transport simulation tool for physically and chemically heterogeneous saturated subsurface environments, *Comput. Geosci.* **19**, 403–422.