ogst
Home All issues Volume 73 (2018) Oil & Gas Science and Technology - Rev. IFP Energies nouvelles, 73 (2018) 73 Full HTML
Numerical methods and HPC
Open Access
Issue
Oil & Gas Science and Technology - Rev. IFP Energies nouvelles
Volume 73, 2018
Numerical methods and HPC
Article Number 73
Number of page(s) 16
DOI https://doi.org/10.2516/ogst/2018033
Published online 11 December 2018

© E. Ahusborde et al., published by IFP Energies nouvelles, 2018

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 CO2 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 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 CO2: 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 SHPCO2 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 CO2 is converted into a source term of liquid CO2 and then the authors study the transport of the dissolved CO2 and the precipitation/dissolution process of minerals. In [11], the authors employ a single-phase reactive flow to model the leaking of CO2-saturated brine in a fractured pathway once supercritical CO2 is totally dissolved. CO2 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 CO2 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 CO2 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 theses references being the complexity of the geochemical system.

During its storage, liquid or supercritical CO2 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 CO2, 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 do 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 CO2 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 Ip and Is such that I = Ip ∪ Is. We note by Ise the set of components of Is that are in equilibrium reactions and by Isk those in kinetic reactions (Is = Ise ∪ Isk).

In order to use a sequential approach to decouple two-phase flow and reactive transport, we introduce two subsets I2φ$ {I}_{2\phi }$ and Irt$ {I}_{{rt}}$. I2φ$ {I}_{2\phi }$ contains all the species present in both phases while Irt$ {I}_{{rt}}$ consists of the remaining ones such that I=I2φ∪Irt$ I={I}_{2\phi }\cup {I}_{{rt}}$. So Nc=card{I}$ {N}_{\mathrm{c}}=\mathrm{card}\{I\}$ is the total number of chemical components involved in Nr=card{Is}$ {N}_{\mathrm{r}}=\mathrm{card}\{{I}_{\mathrm{s}}\}$ chemical reactions (Nr=Ne+Nk$ {N}_{\mathrm{r}}={N}_{\mathrm{e}}+{N}_{\mathrm{k}}$), with Ne=card{Ise}$ {N}_{\mathrm{e}}=\mathrm{card}\{{I}_{\mathrm{se}}\}$ is the number of reactions at equilibrium and Nk=card{Isk}$ {N}_{\mathrm{k}}=\mathrm{card}\{{I}_{\mathrm{sk}}\}$ is the number of kinetic reactions. The chemical system writes:∑j=1NcνijAj=0, i=1,…,Nr,$$ \sum_{j=1}^{{N}_{\mathrm{c}}} {\nu }_{{ij}}{A}_j=0,\enspace i=1,\dots,{N}_r, $$where νij$ {\nu }_{{ij}}$ is the stoichiometric coefficient of the species Aj$ {A}_j$ in the reaction i$ i$.

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

Table 1.

Chemical reactions.

Ano, Cal and Kao refer respectively to Anorthite, Calcite and Kaolinite. In this case, Nc = 14 and Nr=8$ {N}_{\mathrm{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:Is=[H2O(g) CO2(g) HCO3- CO32- OH- Ano Cal Kao],Ip=[H2O(l) CO2(l) H+ Ca2+ Al3+ SiO2(l)],I2φ=[H2O(l) CO2(l) H2O(g) CO2(g)],Irt=[H+ Ca2+ Al3+ SiO2(l) HCO3- CO32- OH- Ano Cal Kao],Ise=[H2O(g) CO2(g) HCO3- CO32- OH-],Isk=[Ano Cal Kao].$$ \begin{array}{ll}& {I}_{\mathrm{s}}=\left[{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{g}\right)}\enspace \mathrm{C}{\mathrm{O}}_{2\left(\mathrm{g}\right)}\enspace \mathrm{HC}{\mathrm{O}}_3^{-}\enspace \mathrm{C}{\mathrm{O}}_3^{2-}\enspace \mathrm{O}{\mathrm{H}}^{-}\enspace \mathrm{Ano}\enspace \mathrm{Cal}\enspace \mathrm{Kao}\right],\\ & {I}_{\mathrm{p}}=\left[{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{l}\right)}\enspace \mathrm{C}{\mathrm{O}}_{2\left(\mathrm{l}\right)}\enspace {\mathrm{H}}^{+}\enspace \mathrm{C}{\mathrm{a}}^{2+}\enspace \mathrm{A}{\mathrm{l}}^{3+}\enspace \mathrm{Si}{\mathrm{O}}_{2\left(\mathrm{l}\right)}\right],\\ & {I}_{2\mathrm{\phi }}=\left[{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{l}\right)}\enspace \mathrm{C}{\mathrm{O}}_{2\left(\mathrm{l}\right)}\enspace {\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{g}\right)}\enspace \mathrm{C}{\mathrm{O}}_{2\left(\mathrm{g}\right)}\right],\\ & {I}_{\mathrm{rt}}=\left[{\mathrm{H}}^{+}\enspace \mathrm{C}{\mathrm{a}}^{2+}\enspace \mathrm{A}{\mathrm{l}}^{3+}\enspace \mathrm{Si}{\mathrm{O}}_{2\left(\mathrm{l}\right)}\enspace \mathrm{HC}{\mathrm{O}}_3^{-}\enspace \mathrm{C}{\mathrm{O}}_3^{2-}\enspace \mathrm{O}{\mathrm{H}}^{-}\enspace \mathrm{Ano}\enspace \mathrm{Cal}\enspace \mathrm{Kao}\right],\\ & {I}_{\mathrm{se}}=\left[{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{g}\right)}\enspace \mathrm{C}{\mathrm{O}}_{2\left(\mathrm{g}\right)}\enspace \mathrm{HC}{\mathrm{O}}_3^{-}\enspace \mathrm{C}{\mathrm{O}}_3^{2-}\enspace \mathrm{O}{\mathrm{H}}^{-}\right],\\ & {I}_{\mathrm{sk}}=\left[\mathrm{Ano}\enspace \mathrm{Cal}\enspace \mathrm{Kao}\right].\end{array} $$

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αj=Kj∏i∈Ip(aαi)νji, j∈Ise,$$ {a}_{\alpha }^j={K}_j\prod_{i\in {I}_{\mathrm{p}}} ({a}_{\alpha }^i{)}^{{\nu }_{{ji}}},\enspace j\in {I}_{\mathrm{se}}, $$(1)where aαj$ {a}_{\alpha }^j$ is the activity of species j$ j$ in its phase α$ \alpha $, Kj$ {K}_j$ is the equilibrium constant of reaction j$ j$. The activity of water and solid species are set equal to 1$ 1$. The activity of aqueous species is often written in terms of molality: alj=γljmljm0$ {a}_{\mathrm{l}}^j={\gamma }_{\mathrm{l}}^j\frac{{m}_{\mathrm{l}}^j}{{m}_0}$, where γlj$ {\gamma }_{\mathrm{l}}^j$ is the activity coefficient for species j$ j$ in the aqueous phase, mlj$ {m}_{\mathrm{l}}^j$ is the molality of species j$ j$ [mol kg−1] in the aqueous phase, and m0$ {m}_0$ is standard molality often taken as 1 mol kg−1. For each aqueous species j$ j$, the molality mlj$ {m}_{\mathrm{l}}^j$ and the mole fraction of species j$ j$ in aqueous phase xlj$ {x}_{\mathrm{l}}^j$ are related by mlj=xljMH2OxlH2O$ {m}_{\mathrm{l}}^j=\frac{{x}_{\mathrm{l}}^j}{{M}_{{\mathrm{H}}_2\mathrm{O}}{x}_{\mathrm{l}}^{{\mathrm{H}}_2\mathrm{O}}}$, where MH2O$ {M}_{{\mathrm{H}}_2\mathrm{O}}$ is the molar mass of water and xlH2O$ {x}_{\mathrm{l}}^{{\mathrm{H}}_2\mathrm{O}}$ is the molar fraction of water in the aqueous phase. The activity of gaseous species is often represented by fugacity: agj=φgjpgjp0$ {a}_{\mathrm{g}}^j={\phi }_{\mathrm{g}}^j\frac{{p}_{\mathrm{g}}^j}{{p}_0}$, where φgj$ {\phi }_{\mathrm{g}}^j$ is the fugacity coefficient for species j$ j$ in the gas phase, pgj$ {p}_{\mathrm{g}}^j$ the partial pressure of gas j$ j$ and p0$ {p}_0$ the standard pressure.

Each kinetic reaction leads to an ordinary differential equation:dcsidt=-ri, i∈Isk,$$ \frac{\mathrm{d}{c}_{\mathrm{s}}^i}{\mathrm{d}t}=-{r}_i,\enspace i\in {I}_{\mathrm{sk}}, $$(2)where csi$ {c}_s^i$ denotes the concentration of mineral i$ i$ and ri$ {r}_i$ is the kinetic reaction rate depending on the activities of the species present in the reaction. Expression for ri$ {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 α∈{l,g,s}$ \alpha \in \{\mathrm{l},\mathrm{g},\mathrm{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:∂∂t(θαρmol, αxαi)+∇⋅(ρmol, αxαiq⃗α)-∇⋅(ρmol, αDα∇xαi)=∑j∈Isνjirj, i∈Ip,$$ \begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}({\theta }_{\alpha }{\rho }_{\mathrm{mol},\enspace \alpha }{x}_{\alpha }^i)+\nabla \cdot ({\rho }_{\mathrm{mol},\enspace \alpha }{x}_{\alpha }^i{\vec{q}}_{\alpha })-\nabla \cdot ({\rho }_{\mathrm{mol},\enspace \alpha }{D}_{\alpha }\nabla {x}_{\alpha }^i)\\ =\sum_{j\in {I}_{\mathrm{s}}} {\nu }_{{ji}}{r}_j,\enspace i\in {I}_{\mathrm{p}},\end{array} $$(3)∂∂t(θαρmol, αxαi)+∇⋅(ρmol, αxαiq⃗α)-∇⋅(ρmol, αDα∇xαi)=-ri, i∈Is,$$ \begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}({\theta }_{\alpha }{\rho }_{\mathrm{mol},\enspace \alpha }{x}_{\alpha }^i)+\nabla \cdot ({\rho }_{\mathrm{mol},\enspace \alpha }{x}_{\alpha }^i{\vec{q}}_{\alpha })-\nabla \cdot ({\rho }_{\mathrm{mol},\enspace \alpha }{D}_{\alpha }\nabla {x}_{\alpha }^i)\\ =-{r}_i,\enspace i\in {I}_{\mathrm{s}},\end{array} $$(4)where θα$ {\theta }_{\alpha }$ [–] denotes the volumetric content of phase α$ \alpha $ (θα=ϕSα$ {\theta }_{\alpha }=\phi {\mathrm{S}}_{\alpha }$, ϕ [–] being the porosity of the medium and Sα$ {\mathrm{S}}_{\alpha }$ [–] the saturation of phase α$ \alpha $ if α∈{l,g}$ \alpha \in \{l,g\}$ and θs=1-ϕ$ {\theta }_{\mathrm{s}}=1-\mathrm{\phi }$), ρmol, α$ {\rho }_{\mathrm{mol},\enspace \alpha }$ [mol m−3] is the molar density of phase α$ \alpha $, xαi$ {x}_{\alpha }^i$ is the molar fraction of species i$ i$ in phase α$ \alpha $, Dα$ {D}_{\alpha }$ [m2 s−1] denotes the diffusivity of the phase α:Dα=θα(DmI+dL|q⃗α|I+(dL-dT)q⃗αq⃗αT|q⃗α|).$$ {D}_{\alpha }={\theta }_{\alpha }\left({D}_{\mathrm{m}}\mathbb{I}+{d}_{\mathrm{L}}|{\vec{q}}_{\alpha }|\mathbb{I}+({d}_{\mathrm{L}}-{d}_{\mathrm{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 $. Dm$ {D}_{\mathrm{m}}$ [m2 s−1] is the molecular diffusion, dL [m] and dT [m] are the magnitudes of longitudinal and transverse dispersion respectively and q⃗α$ {\vec{q}}_{\alpha }$ [m s−1] is the Darcy velocity of phase α$ \alpha $, expressed as follows:q⃗α=-krα(Sl)μαK(∇Pα-ραg⃗),$$ {\vec{q}}_{\alpha }=-\frac{{k}_{{r\alpha }}({S}_l)}{{\mu }_{\alpha }}\mathbb{K}(\nabla {P}_{\alpha }-{\rho }_{\alpha }\vec{g}), $$(5)where krα$ {k}_{{r\alpha }}$ [–] denotes the relative permeability of phase α$ \alpha $, μα$ {\mu }_{\alpha }$ [Pa s] is the dynamic viscosity of phase α$ \alpha $, K$ \mathbb{K}$ [m2] is the absolute permeability tensor, Pα [Pa] is the pressure of phase α, ρα$ {\rho }_{\alpha }$ [kg m−3] is the mass density of phase α and g⃗$ \vec{g}$ [m s−2] is the gravitational acceleration.

The phase pressures are connected by the capillary pressure law:Pc(Sl)=Pg-Pl.$$ {P}_{\mathrm{c}}({S}_{\mathrm{l}})={P}_{\mathrm{g}}-{P}_{\mathrm{l}}. $$(6)

Finally, rj$ {r}_j$ [mol m−3 s−1] is the rate of reaction j$ j$ (it can be equilibrium or kinetic).

We choose to reformulate equations (3) and (4) in term of molar concentration cαi=ρmol, αxαi$ {c}_{\alpha }^i={\rho }_{\mathrm{mol},\enspace \alpha }{x}_{\alpha }^i$. By neglecting the gradient of molar density ∇ρmol, α$ \nabla {\rho }_{\mathrm{mol},\enspace \alpha }$ resulting from the diffusive flux, equations (3) and (4) write:∂∂t(θαcαi)+∇⋅(cαiq⃗α)-∇⋅(Dα∇cαi)=∑j∈Isνjirj, i∈Ip,$$ \frac{\mathrm{\partial }}{\mathrm{\partial }t}({\theta }_{\alpha }{\mathrm{c}}_{\alpha }^i)+\nabla \cdot ({\mathrm{c}}_{\alpha }^i{\vec{q}}_{\alpha })-\nabla \cdot ({D}_{\alpha }\nabla {\mathrm{c}}_{\alpha }^i)=\sum_{j\in {I}_s} {\nu }_{{ji}}{r}_j,\enspace i\in {I}_p, $$(7)∂∂t(θαcαi)+∇⋅(cαiq⃗α)-∇⋅(Dα∇cαi)=-ri, i∈Is.$$ \frac{\mathrm{\partial }}{\mathrm{\partial }t}({\theta }_{\alpha }{\mathrm{c}}_{\alpha }^i)+\nabla \cdot ({\mathrm{c}}_{\alpha }^i{\vec{q}}_{\alpha })-\nabla \cdot ({D}_{\alpha }\nabla {\mathrm{c}}_{\alpha }^i)=-{r}_i,\enspace i\in {I}_s. $$(8)

For sake of simplicity, we introduce the advection-diffusion operator:Lα(c)=∇⋅(cq⃗α)-∇⋅(Dα∇c).$$ {L}_{\alpha }(\mathrm{c})=\nabla \cdot (c{\vec{q}}_{\alpha })-\nabla \cdot ({D}_{\alpha }\nabla c). $$(9)

We can remark that, since Ds=0$ {D}_{\mathrm{s}}=0$ and q⃗s=0⃗$ {\vec{q}}_{\mathrm{s}}=\vec{0}$, so that Ls(.)=0$ {L}_{\mathrm{s}}(.)=0$.

In order to eliminate the reaction rates rj$ {r}_j$ in equations (7), we make linear combinations between equations (8) with each equation (7). This introduces Np$ {N}_{\mathrm{p}}$ new conservation laws that write:∂∂t(θαcαi+∑j∈Isθανjicαj)+Lα(cαi)+∑j∈IsLα(νjicαj)=0, i∈Ip.$$ \begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\theta }_{\alpha }{c}_{\alpha }^i+\sum_{j\in {I}_{\mathrm{s}}} {\theta }_{\alpha }{\nu }_{{ji}}{c}_{\alpha }^j\right)+{L}_{\alpha }({c}_{\alpha }^i)+\sum_{j\in {I}_{\mathrm{s}}} {L}_{\alpha }({\nu }_{{ji}}{c}_{\alpha }^j)\\ =0,\enspace i\in {I}_{\mathrm{p}}.\end{array} $$(10)

To retrieve the same number of equations as there are unknowns, the Ns$ {N}_{\mathrm{s}}$ equations (8) are replaced by Ne$ {N}_{\mathrm{e}}$ mass actions laws defined by (1) corresponding to the equilibrium reactions and Nk$ {N}_{\mathrm{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:∂∂t(θlcli+θgcgi¯)+Ll(cli)+Lg(cgi¯)+∑j∈Is\I2φνji[∂∂t(θαcαj)+Lα(cαj)]=0, i∈Ip∩I2φ,$$ \begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\theta }_{\mathrm{l}}{c}_l^i+{\theta }_{\mathrm{g}}{c}_{\mathrm{g}}^{\bar{i}}\right)+{L}_{\mathrm{l}}\left({c}_{\mathrm{l}}^i\right)+{L}_g\left({c}_{\mathrm{g}}^{\bar{i}}\right)\\ +\sum_{j\in {I}_{\mathrm{s}}\backslash {I}_{2\phi }} {\nu }_{{ji}}\left[\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\theta }_{\alpha }{c}_{\alpha }^j\right)+{L}_{\alpha }\left({c}_{\alpha }^j\right)\right]=0,\enspace i\in {I}_{\mathrm{p}}\cap {I}_{2\phi },\end{array} $$(11)agi¯=Kiali, i¯∈Ise∩I2φ.$$ {a}_{\mathrm{g}}^{\bar{i}}={K}_i{a}_{\mathrm{l}}^i,\enspace \bar{i}\in {I}_{\mathrm{se}}\cap {I}_{2\phi }. $$(12)

Equation (11) corresponds to equation (10) for the species present in liquid and gas phases where i$ i$ and i¯$ \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$ j$. Equation (12) is the solubility law where ali$ {a}_{\mathrm{l}}^i$ and agi¯$ {a}_{\mathrm{g}}^{\bar{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:∂∂t(θαcαi+∑j∈Isθανjicαj)+Lα(cαi)+∑j∈IsLα(νjicαj)=0, i∈Ip∩Irt,$$ \begin{array}{c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\theta }_{\alpha }{c}_{\alpha }^i+\sum_{j\in {I}_{\mathrm{s}}} {\theta }_{\alpha }{\nu }_{{ji}}{c}_{\alpha }^j\right)+{L}_{\alpha }({c}_{\alpha }^i)+\sum_{j\in {I}_{\mathrm{s}}} {L}_{\alpha }({\nu }_{{ji}}{c}_{\alpha }^j)\\ =0,\enspace i\in {I}_{\mathrm{p}}\cap {I}_{{rt}},\end{array} $$(13)aαj=Kj∏i∈Ip(aαi)νji, j∈Ise\I2φ,$$ {a}_{\alpha }^j={K}_j\prod_{i\in {I}_{\mathrm{p}}} ({a}_{\alpha }^i{)}^{{\nu }_{{ji}}},\enspace j\in {I}_{\mathrm{se}}\backslash {I}_{2\phi }, $$(14)dcsidt=-KisAis(1-Ki∏j∈Ip(alj)νij), i∈Isk,$$ \frac{\mathrm{d}{c}_{\mathrm{s}}^i}{\mathrm{d}t}=-{K}_i^s{A}_i^s\left(1-{K}_i\prod_{j\in {I}_p} ({a}_l^j{)}^{{\nu }_{{ij}}}\right),\enspace i\in {I}_{{sk}}, $$(15)where Kis$ {K}_i^s$ and Ais$ {A}_i^s$ are respectively the kinetic-rate constant [mol m−2 s−1] and the reactive surface area [m2 m−3] of mineral i$ 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∈I2φ$ 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∈Is/I2φ$ j\in {I}_{\mathrm{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 Vk$ {V}_{\mathrm{k}}$ (see Fig. 1) and evaluating the fluxes at the interface γkl$ {\gamma }_{{kl}}$ between two adjacent elements Vk$ {V}_{\mathrm{k}}$ and Vl$ {V}_{\mathrm{l}}$.

  • Spatial discretization of two-phase compositional flow:

∫Vk∂∂t{θlcli+θgcgi¯} dV+∑l∈V(k)∫γkl{cliq⃗l+cgi¯q⃗g}⋅n⃗kl dγ-∑l∈V(k)∫γkl{Dl∇cli+Dg∇cgi¯}⋅n⃗kl dγ=-∑j∈Is\I2φνji(∫Vk∂∂t(θαcαj) dV+∑l∈V(k)∫γkl{cαjq⃗α-Dα∇cαj}⋅n⃗kl d γ),i∈Ip∩I2φ,$$ \begin{array}{c}{\int }_{{V}_k} \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left\{{\theta }_{\mathrm{l}}{c}_{\mathrm{l}}^i+{\theta }_{\mathrm{g}}{c}_{\mathrm{g}}^{\bar{i}}\right\}\enspace \mathrm{d}V+\sum_{l\in V(k)} {\int }_{{\gamma }_{{kl}}} \left\{{c}_{\mathrm{l}}^i{\vec{q}}_l+{c}_{\mathrm{g}}^{\bar{i}}{\vec{q}}_{\mathrm{g}}\right\}\cdot {\vec{n}}_{{kl}}\enspace \mathrm{d}\gamma \\ -\sum_{l\in V(k)} {\int }_{{\gamma }_{{kl}}} \left\{{D}_l\nabla {c}_l^i+{D}_g\nabla {c}_g^{\bar{i}}\right\}\cdot {\vec{n}}_{{kl}}\enspace \mathrm{d}\gamma \\ \begin{array}{c}=-\sum_{j\in {I}_s\backslash {I}_{2\phi }} {\nu }_{{ji}}\left({\int }_{{V}_k} \frac{\mathrm{\partial }}{\mathrm{\partial }t}({\theta }_{\alpha }{c}_{\alpha }^j)\enspace \mathrm{d}V\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}}\enspace d\enspace \gamma \right),i\in {I}_p\cap {I}_{2\phi },\end{array}\end{array} $$(16){ai¯}k=Ki{ai}k, i¯∈Ise∩I2φ,$$ {\left\{{a}^{\bar{i}}\right\}}_{\mathrm{k}}={K}_i{\left\{{a}^i\right\}}_{\mathrm{k}},\enspace \bar{i}\in {I}_{\mathrm{se}}\cap {I}_{2\phi }, $$(17)q⃗α=-K λα(Sl)(∇Pα-ρα g⃗),$$ {\vec{q}}_{\alpha }=-\mathbb{K}\enspace {\lambda }_{\alpha }({S}_l)(\nabla {P}_{\alpha }-{\rho }_{\alpha }\enspace \vec{g}), $$(18)where λα(Sl)=krα(Sl)μα$ {\lambda }_{\alpha }({S}_l)=\frac{{k}_{{r\alpha }}({S}_l)}{{\mu }_{\alpha }}$ is the phase mobility,

  • Spatial discretization of reactive transport:

∫Vk∂∂t{θαcαi+∑j∈Isνjiθαcαj} dV+∑l∈V(k)∫γkl{cαiq⃗α+∑j∈Isνjicαjq⃗α}⋅n⃗kl dγ-∑l∈V(k)∫γkl{Dα∇cαi+∑j∈IsνjiDα∇cαj}⋅n⃗kldγ=0, i∈Ip∩Irt,$$ \begin{array}{c}\begin{array}{c}{\int }_{{V}_k} \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left\{{\theta }_{\alpha }{c}_{\alpha }^i+\sum_{j\in {I}_s} {\nu }_{{ji}}{\theta }_{\alpha }{c}_{\alpha }^j\right\}\enspace \mathrm{d}V+\\ \sum_{l\in V(k)} {\int }_{{\gamma }_{{kl}}} \left\{{c}_{\alpha }^i{\vec{q}}_{\alpha }+\sum_{j\in {I}_s} {\nu }_{{ji}}{c}_{\alpha }^j{\vec{q}}_{\alpha }\right\}\cdot {\vec{n}}_{{kl}}\enspace \mathrm{d}\gamma \\ -\sum_{l\in V(k)} {\int }_{{\gamma }_{{kl}}} \left\{{D}_{\alpha }\nabla {c}_{\alpha }^i+\sum_{j\in {I}_{\mathrm{s}}} {\nu }_{{ji}}{D}_{\alpha }\nabla {c}_{\alpha }^j\right\}\cdot {\vec{n}}_{{kl}}\mathrm{d}\gamma \end{array}\\ =0,\enspace i\in {I}_{\mathrm{p}}\cap {I}_{{rt}},\end{array} $$(19){aαj}k=Kj∏i∈Ip{(aαi)νji}k,j∈Ise/I2φ,$$ {\left\{{a}_{\alpha }^j\right\}}_{\mathrm{k}}={K}_j\prod_{i\in {I}_{\mathrm{p}}} {\left\{({a}_{\alpha }^i{)}^{{\nu }_{{ji}}}\right\}}_k,j\in {I}_{\mathrm{se}}/{I}_{2\phi }, $$(20)d{csi}kdt+KisAis(1-Ki∏j∈Ip{(alj)νij}k)=0, i∈Isk,$$ \frac{\mathrm{d}{\left\{{c}_s^i\right\}}_k}{\mathrm{d}t}+{K}_i^s{A}_i^s(1-{K}_i\prod_{j\in {I}_p} {\left\{({a}_l^j{)}^{{\nu }_{{ij}}}\right\}}_k)=0,\enspace i\in {I}_{\mathrm{sk}}, $$(21)where n⃗kl$ {\vec{n}}_{{kl}}$ denotes the unit outer normal to γkl$ {\gamma }_{{kl}}$, V(k)$ V(k)$ is the set of adjacent elements of Vk$ {V}_k$ and {⋅}k$ {\left\{\cdot \right\}}_k$ represents the average values on Vk$ {V}_k$. By using the implicit Euler scheme for the time discretization and due to the fact that the primary unknowns (Sα$ {\mathrm{S}}_{\alpha }$, Pα$ {\mathrm{P}}_{\alpha }$, cαi$ {\mathrm{c}}_{\alpha }^i$) and the physical parameters are constant on each element Vk$ {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:

|Vk|Δtn({θlcli+θgcgi¯}kn+1-{θlcli+θgcgi¯}kn)+∑l∈V(k)|γkl|({cli}kln+1 {q⃗l}kln+1+{cgi¯}kln+1 {q⃗g}kln+1)⋅n⃗kl-∑l∈V(k)|γkl|({Dl}kln+1{∇cli}kln+1+{Dg}kln+1{∇cgi¯}kln+1)⋅ n⃗kl=-∑j∈Is/I2φνji|Vk|Δtn({θαcαj}kn-{θαcαj}kn-1)-∑j∈Is/I2φνji∑l∈V(k)|γkl|{cαj}kln {q⃗α}kln+1⋅ n⃗kl+∑j∈Is/I2φνji∑l∈V(k)|γkl|{Dα}kln+1{∇cαj}kln⋅ n⃗kl,i∈Ip∩I2φ,$$ \begin{array}{c}\frac{|{V}_k|}{\Delta {t}^n}\left({\left\{{\theta }_{\mathrm{l}}{c}_{\mathrm{l}}^i+{\theta }_{\mathrm{g}}{c}_{\mathrm{g}}^{\bar{i}}\right\}}_k^{n+1}-{\left\{{\theta }_{\mathrm{l}}{c}_{\mathrm{l}}^i+{\theta }_{\mathrm{g}}{c}_{\mathrm{g}}^{\bar{i}}\right\}}_k^n\right)\\ +\sum_{l\in V(k)}|{\gamma }_{{kl}}|\left({\left\{{c}_{\mathrm{l}}^i\right\}}_{{kl}}^{n+1}\enspace {\left\{{\vec{q}}_{\mathrm{l}}\right\}}_{{kl}}^{n+1}+{\left\{{c}_{\mathrm{g}}^{\bar{i}}\right\}}_{{kl}}^{n+1}\enspace {\left\{{\vec{q}}_{\mathrm{g}}\right\}}_{{kl}}^{n+1}\right)\cdot {\vec{n}}_{{kl}}\\ \begin{array}{c}-\sum_{l\in V(k)}|{\gamma }_{{kl}}|\left({\left\{{D}_{\mathrm{l}}\right\}}_{{kl}}^{n+1}{\left\{\nabla {c}_{\mathrm{l}}^i\right\}}_{{kl}}^{n+1}+{\left\{{D}_{\mathrm{g}}\right\}}_{{kl}}^{n+1}{\left\{\nabla {c}_{\mathrm{g}}^{\bar{i}}\right\}}_{{kl}}^{n+1}\right)\cdot \enspace {\vec{n}}_{{kl}}\\ =-\sum_{j\in {I}_{\mathrm{s}}/{I}_{2\phi }}{\nu }_{{ji}}\frac{|{V}_k|}{\Delta {t}^n}\left({\left\{{\theta }_{\alpha }{c}_{\alpha }^j\right\}}_k^n-{\left\{{\theta }_{\alpha }{c}_{\alpha }^j\right\}}_k^{n-1}\right)\\ \begin{array}{c}-\sum_{j\in {I}_{\mathrm{s}}/{I}_{2\phi }}{\nu }_{{ji}}\sum_{l\in V(k)}|{\gamma }_{{kl}}|{\left\{{c}_{\alpha }^j\right\}}_{{kl}}^n\enspace {\left\{{\vec{q}}_{\alpha }\right\}}_{{kl}}^{n+1}\cdot \enspace {\vec{n}}_{{kl}}\\ +\sum_{j\in {I}_{\mathrm{s}}/{I}_{2\phi }}{\nu }_{{ji}}\sum_{l\in V(k)}|{\gamma }_{{kl}}|{\left\{{D}_{\alpha }\right\}}_{{kl}}^{n+1}{\left\{\nabla {c}_{\alpha }^j\right\}}_{{kl}}^n\cdot \enspace {\vec{n}}_{{kl}},i\in {I}_p\cap {I}_{2\phi },\end{array}\end{array}\end{array} $$(22){agi¯}kn+1=Ki{ali}kn+1,i¯∈Ise∩I2φ,$$ {\left\{{a}_{\mathrm{g}}^{\bar{i}}\right\}}_k^{n+1}={K}_i{\left\{{a}_{\mathrm{l}}^i\right\}}_k^{n+1},\bar{i}\in {I}_{\mathrm{se}}\cap {I}_{2\phi }, $$(23){q⃗α}kln+1=-{K}kl {λα(Sl)}kln+1({∇Pα}kln+1-{ρα}kln+1 g⃗),$$ \{{\vec{q}}_{\alpha }{\}}_{{kl}}^{n+1}=-\{\mathbb{K}{\}}_{{kl}}\enspace \{{\lambda }_{\alpha }({S}_l){\}}_{{kl}}^{n+1}\left(\{\nabla {P}_{\alpha }{\}}_{{kl}}^{n+1}-\{{\rho }_{\alpha }{\}}_{{kl}}^{n+1}\enspace \vec{g}\right), $$(24)

  • Cell-centred FV scheme for reactive transport:

|Vk|Δtn({θαcαi+∑j∈Isνjiθαcαj}kn+1-{θαcαi+∑j∈Isνjiθαcαj}kn)+∑l∈V(k)|γkl|({cαi}kln+1{q⃗α}kln+1+∑j∈Isνji{cαj}kln+1{q⃗α}kln+1)⋅n⃗kl-∑l∈V(k)|γkl|({Dα}kln+1{∇cαi}kln+1+∑j∈Isνji{Dα}kln+1{∇cαj}kln+1)⋅n⃗kl=0,i∈Ip∩Irt,$$ \begin{array}{c}\frac{|{V}_{\mathrm{k}}|}{\Delta {t}^n}\left({\left\{{\theta }_{\alpha }{c}_{\alpha }^i+\sum_{j\in {I}_{\mathrm{s}}} {\nu }_{{ji}}{\theta }_{\alpha }{c}_{\alpha }^j\right\}}_k^{n+1}-{\left\{{\theta }_{\alpha }{c}_{\alpha }^i+\sum_{j\in {I}_{\mathrm{s}}} {\nu }_{{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}_{\mathrm{s}}} {\nu }_{{ji}}{\left\{{c}_{\alpha }^j\right\}}_{{kl}}^{n+1}{\left\{{\vec{q}}_{\alpha }\right\}}_{{kl}}^{n+1}\right)\cdot {\vec{n}}_{{kl}}\\ \begin{array}{c}-\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.\\ \left.+\sum_{j\in {I}_{\mathrm{s}}} {\nu }_{{ji}}{\left\{{D}_{\alpha }\right\}}_{{kl}}^{n+1}{\left\{\nabla {c}_{\alpha }^j\right\}}_{{kl}}^{n+1}\right)\cdot {\vec{n}}_{{kl}}=0,i\in {I}_p\cap {I}_{{rt}},\end{array}\end{array} $$(25){aαj}kn+1=Kj∏i∈Ip{(aαi)νji}kn+1,j∈Ise/I2φ,$$ {\left\{{a}_{\alpha }^j\right\}}_k^{n+1}={K}_j\prod_{i\in {I}_p} {\left\{({a}_{\alpha }^i{)}^{{\nu }_{{ji}}}\right\}}_k^{n+1},j\in {I}_{\mathrm{se}}/{I}_{2\phi }, $$(26){csi}kn+1={csi}kn-ΔtnKisAis(1-Ki∏j∈Ip{(alj)νij}kn+1),i∈Isk.$$ {\left\{{c}_s^i\right\}}_k^{n+1}={\left\{{c}_s^i\right\}}_k^n-\Delta {t}^n{K}_i^s{A}_i^s\left(1-{K}_i\prod_{j\in {I}_p} {\left\{({a}_l^j{)}^{{\nu }_{{ij}}}\right\}}_k^{n+1}\right),i\in {I}_{{sk}}. $$(27)

thumbnail Fig. 1.

Discretization by the cell-centered finite volume method.

Now to define the finite volume scheme it is enough to approach the convective and the diffusive fluxes on the interfaces γkl$ {\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α$ {S}_{\alpha }$, cαi$ {c}_{\alpha }^i$, Pα$ {P}_{\alpha }$, λα$ {\lambda }_{\alpha }$) are evaluated implicitly and upstream at the interface γkl$ {\gamma }_{{kl}}$ between two adjacent elements as:{⋅}kln+1={{⋅}kn+1 if {q⃗α}kln+1⋅n⃗kl>0,{⋅}ln+1 else.$$ \{\cdot {\}}_{{kl}}^{n+1}=\left\{\begin{array}{ll}& \{\cdot {\}}_k^{n+1}\enspace {if}\enspace \{{\vec{q}}_{\alpha }{\}}_{{kl}}^{n+1}\cdot {\vec{n}}_{{kl}}>0,\\ & \{\cdot {\}}_l^{n+1}\enspace \mathrm{else}.\end{array}\right. $$(28)

The gradient operators on the interfaces γkl are calculated by a P1/Q1$ {\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 calculated the absolute permeability {K}kl$ \{\mathbb{K}{\}}_{{kl}}$ and the diffusion coefficients {Dα}kln+1$ \{{D}_{\alpha }{\}}_{{kl}}^{n+1}$ at the interface γkl. {ρα}kln+1$ \{{\rho }_{\alpha }{\}}_{\mathrm{kl}}^{n+1}$ is computed as the arithmetic average of two elements Vk$ {V}_{\mathrm{k}}$ and Vl$ {V}_{\mathrm{l}}$.

4 Numerical simulations

4.1 Implementation

The physical models described above are built in DuMuX [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 DuMuX called 2pnc$ 2{pnc}$ (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 DuMuX framework a one-phase multicomponent transport module named 1pmc-react. In this context, we have implemented a new 1pmc$ 1{pmc}$ 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 CO2 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 CO2 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.

Table 2.

Chemical reactions.

Table 3.

Mineral, precipitation and dissolution parameters.

Initial molalities [mol kg−1] are exhibited in Table 4. The activity of each aqueous species is written in terms of molality as:ali=γimli,$$ {a}_l^i={\gamma }^i{m}_l^i, $$(29)where γi 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]:γi=-AziI1+ȦiBI+ḂI,$$ {\gamma }^i=-\frac{A{z}_i\sqrt{I}}{1+{\dot{A}}_iB\sqrt{I}}+\dot{B}I, $$(30)where A, B and Ḃ$ \dot{B}$ are constants depending on temperature, zi$ {z}_i$ is the ion electrical charge for species i and Ȧi$ {\dot{A}}_i$ is the ion size of species i$ i$. Finally, I$ I$ is the ionic strength of the aqueous phase and expresses as I=0.5∑imlizi2 $ I=0.5{\sum }_i{m}_l^i{z}_i^2\enspace $.

Table 4.

Input parameters for each ion.

Tables 4 and 5 exhibit the parameters for the B-dot model. Initially, the system contains 10 g of calcite and 100 kg of H2O.

Table 5.

B-dot model parameters from the EQ3/6 database [66].

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.

thumbnail Fig. 2.

Evolution of the molality of the different components versus time.

thumbnail Fig. 3.

Evolution of the number of moles of dissolved calcite.

4.2.2 Application to geological CO2 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 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 CO2 is implemented according to [65]. Mineral data for the kinetic reactions are summarized in Table 7.

Table 6.

Chemical reactions.

Table 7.

Mineral, precipitation and dissolution parameters.

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 CO2 stream at constant rate during the first 20 years. After the 20 years injection period, a total of 18.6 × 109 kg of CO2 is injected.

As initial conditions for the two-phase two-component H2O-CO2 flow we have used hydrostatic condition for liquid pressure Pl$ {P}_{\mathrm{l}}$, initial liquid saturation Sl=1$ {S}_{\mathrm{l}}=1$ and initial CO2 molality in liquid phase equals 3.55 × 10−3 mol kg−1. Initial conditions for the reactive transport problem and parameters for the B-dot model used as activity model are shown in Table 8.

Table 8.

Input parameters for each ion.

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

Table 9.

Physical parameters for the test case of CO2 injection.

Figure 4 displays concentrations of calcite, anorthite and kaolinite at 20 and 2000 years for a grid composed of 1.6 × 105 elements (100 × 100 × 16). Initially, their concentration were respectively 238, 87 and 88 mol m−3. We can see that calcite is dissolved near the injection of CO2 and precipitated far from the injection while anorthite and kaolinite are respectively dissolved and precipitated everywhere.

thumbnail Fig. 4.

Profiles of mineral concentrations.

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

thumbnail Fig. 5.

Profiles of pH, aqueous CO2 molality and gas saturation.

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 anorthite is dissolved. The distribution of CO2 shows that most of CO2 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].

thumbnail Fig. 6.

Mineral net molar changes and evolution of CO2 in time.

thumbnail Fig. 7.

CPU time and strong parallel efficiency as a function of the number of processors.

thumbnail Fig. 8.

CPU time and weak parallel efficiency as a function of the number of processors.

Table 10.

Final carbon distribution at 2000 years.

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.

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 × 105, 1.28 × 106 and 5.76 × 106 elements corresponding respectively to 1.92 × 106, 1.536 × 107 and 6.912 × 107 degrees of freedom (dof). The dashed lines represent an ideal behaviour.

Strong efficiency is given by:SE(N)= CPU time on p processors ×p CPU time on N processors ×N,$$ {SE}(N)=\frac{\enspace \mathrm{CPU}\enspace \mathrm{time}\enspace \mathrm{on}\enspace p\enspace \mathrm{processors}\enspace \times p}{\enspace \mathrm{CPU}\enspace \mathrm{time}\enspace \mathrm{on}\enspace N\enspace \mathrm{processors}\enspace \times N}, $$(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$ 24$ and 48$ 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 × 107 and 6.912 × 107 dof. For the simulation with 1.92 × 106 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.

Weak efficiency is given by:WE(N)= CPU time on p processors  CPU time on N processors ,$$ \mathrm{WE}(N)=\frac{\enspace \mathrm{CPU}\enspace \mathrm{time}\enspace \mathrm{on}\enspace p\enspace \mathrm{processors}\enspace }{\enspace \mathrm{CPU}\enspace \mathrm{time}\enspace \mathrm{on}\enspace N\enspace \mathrm{processors}\enspace }, $$(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 εSIA$ {\epsilon }_{\mathrm{SIA}}$. In the sequel, tolerances for the Newton method and iterative algorithm are respectively εNewton = 10−8 and εSIA = 10−8.

Figure 9 compares the evolution of the molalities of H+ and Ca2+ obtained with the DSA and the SIA close to the injection with a mesh composed of 40 000 cells during the first year of injection. We can observe that the results are in good accordance and that both methods provide comparable results.

thumbnail Fig. 9.

Comparison of the molalities of H+ and Ca2+ obtained with the DSA and the SIA.

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.

Table 11.

CPU time (s) for the DSA and the SIA.

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

  • Niemi A., Bear J., Bensabat J. (2017) Geological storage of CO2 in deep saline formations, Springer. [Google Scholar]
  • Zhang F., Yeh G.T., Parker J.C. (2012) Groundwater reactive transport models, Bentham e-books [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • Jiang X. (2011) A review of physical modelling and numerical simulation of long-term geological storage of CO2, Appl. Energy 88, 3557–3566. [Google Scholar]
  • 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. [Google Scholar]
  • Al-Khoury R., Bundschuh J. (2014) Computational models for CO2 geo-sequestration and compressed air energy storage, Sustainable Energy Developments, CRC Press. [Google Scholar]
  • Bear J., Carrera J. (2017) Mathematical modeling of CO2 storage in a geological formation, Springer. [Google Scholar]
  • Haeberlein F. (2011) Time space domain decomposition methods for reactive transport – application to CO2 geological storage, PhD Thesis, Université Paris-Nord – Paris XIII. [Google Scholar]
  • Lagneau V., Pipart A., Catalette H. (2005) Reactive transport modelling of CO2 sequestration in deep saline aquifers, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 60, 231–247. [CrossRef] [Google Scholar]
  • Haeberlein F., Michel A., Halpern L.(2009) A test case for multi-species reactive-transport in heterogeneous porous media applied to CO2 geological storage. https://www.ljll.math.upmc.fr/mcparis09/Files/haeberlein_poster.pdf. [Google Scholar]
  • Ahmad N., Wörman A., Bottacin-Busolin A., Sanchez-Vila X. (2015) Reactive transport modeling of leaking CO2-saturated brine along a fractured pathway, Int. J. Greenh. Gas Con. 42, 672–689. [CrossRef] [Google Scholar]
  • Pool M., Carrera J., Vilarrasa V., Silva O., Ayora C. (2013) Dynamics and design of systems for geological storage of dissolved CO2, Adv. Water Resour. 62, 533–542. [Google Scholar]
  • Ahmad N., Wörman A., Sanchez-Vila X., Jarsjö J., Bottacin-Busolin A., Hellevang H. (2016) Injection of CO2 saturated brine in geological reservoir: A way to enhanced storage safety, Int. J. Greenh. Gas Con. 54, 129–144. [CrossRef] [Google Scholar]
  • Nicot J.P., Hosseini S.A., Solano S.V. (2011) Are single-phase flow numerical models sufficient to estimate pressure distribution in CO2 sequestration projects?, Energy Procedia 4, 3919–3926. [Google Scholar]
  • Audigane P., Gaus I., Czernichowski-Lauriol I., Pruess K., Xu T. (2007) Two-dimensional reactive transport modelling of CO2 injection in a saline aquifer at the Sleipner site, North Sea, Am. J. Sci. 307, 974–1008. [Google Scholar]
  • Fan Y., Durlofsky L.J., Tchelepi H.A. (2012) A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO2 storage simulations, Adv. Water Resour. 42, 47–61. [Google Scholar]
  • Leal A.M.M., Blunt M.J., LaForce T.C. (2013) A robust and efficient numerical method for multiphase equilibrium calculations: Application to CO2-brine-rock systems at high temperatures, pressures and salinities, Adv. Water Resour. 62, 409–430. [Google Scholar]
  • Nghiem L., Sammon P., Grabenstetter J., Ohkuma H. (2004) Modeling CO2 storage in aquifers with a fully-coupled geochemical eos compositional simulator, SPE - DOE Improved Oil Recovery Symposium Proceedings. [Google Scholar]
  • 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. [Google Scholar]
  • 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 CO2 storage, Adv. Water Resour. 62, 431–441. [Google Scholar]
  • Thibeau S., Nghiem L.X., Ohkuma H. (2007) A modeling study of the role of selected minerals in enhancing CO2 mineralization during CO2 aquifer storage, Proceedings - SPE Annual Technical Conference and Exhibition 2, 906–922. [Google Scholar]
  • Huet B.M., Prevost J.H., Scherer G.W. (2010) Quantitative reactive transport modeling of portland cement in CO2-saturated water, Int. J. Greenh. Gas Con. 4, 561–574. [CrossRef] [Google Scholar]
  • Jacquemet N., Pironon J., Lagneau V., Saint-Marc J. (2012) Armouring of well cement in H2S-CO2 saturated brine by calcite coating - experiments and numerical modelling, Appl. Geochem. 27, 782–795. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [PubMed] [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [PubMed] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [PubMed] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [PubMed] [Google Scholar]
  • Appelo C.A.J., Postma D. (2005) Geochemistry, Groundwater and Pollution, 2nd edn., Taylor & Francis. [CrossRef] [Google Scholar]
  • Bear J., Cheng A.H.-D. (2010) Modeling groundwater flow and contaminant transport, Springer. [Google Scholar]
  • Zheng C., Bennett G.D. (2002) Applied contaminant transport modeling, John Wiley and Sons, New York. [Google Scholar]
  • Lichtner P.C. (1985) Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems, Geochim. Cosmochim. Acta 49, 779–800. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • Steefel C.I., MacQuarrie K.T.B. (1996) Approaches to modeling of reactive transport in porous media, Rev. Mineral. 34, 82–129. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • Valocchi A.J., Malmstead M. (1992) Accuracy of operator splitting for advection-dispersion-reaction problems, Water Resour. Res. 28, 1471–1476. [Google Scholar]
  • Carrayrou J., Kern M., Knabner P. (2010) Reactive transport benchmark of MoMaS, Comput. Geosci. 14, 385–392. [Google Scholar]
  • Carrayrou J. (2010) Looking for some reference solutions for the reactive transport benchmark of MoMaS with SPECY, Comput. Geosci. 14, 393–403. [Google Scholar]
  • Lagneau V., van der Lee J. (2010) HYTEC results of the MoMaS reactive transport benchmark, Comput. Geosci. 14, 435–449. [Google Scholar]
  • Amir L., Kern M. (2010) A global method for coupling transport with chemistry in heterogeneous porous media, Comput. Geosci. 14, 465–481. [Google Scholar]
  • de Dieuleveult C., Erhel J. (2010) A global approach to reactive transport: Application to the MoMaS benchmark, Comput. Geosci. 14, 451–464. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • Erhel J., Sabit S. (2017) Analysis of a global reactive transport model and results for the MoMaS benchmark, Math. Comput. Simulat. 137, 286–298. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • 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 CO2, ESAIM: Proc. Surveys 50, 21–39. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • DuMuX (2018) DUNE for multi-Phase, Component, Scale, Physics, … flow and transport in porous media, https://www.dumux.org, last accessed February 1, 2018. [Google Scholar]
  • Flemisch B., Darcis M., Erbertseder K., Faigle B., Lauser A., Mosthaf K., Muthing S., Nuske P., Tatomir A., Wolf M., Helmig R. (2011) DuMuX: DUNE for multi-{Phase, Component, Scale, Physics, …} flow and transport in porous media, Adv. Water Resour. 34, 1102–1112. [Google Scholar]
  • Helmig R. (1997) Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems, Springer. [Google Scholar]
  • GSL - GNU Scientific LibraryMultidimensional Root-Finding. https://www.gnu.org/software/gsl/, Last accessed February 1, 2018. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • Spycher N., Pruess K. (2005) CO2-H2O mixtures in the geological sequestration of CO2. II. Partitioning in chloride brines at 12–100 °C and up to 600 bar, Geochim. et Cosmochim. Acta 69, 3309–3320. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • 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. [Google Scholar]
  • Adams J.J., Bachu S. (2002) Equations of state for basin geofluids: algorithm review and intercomparison for brines, Geofluids 2, 257–271. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • Fenghour A., Wakeham W.A., Vesovic V. (1998) The viscosity of carbon dioxide, J. Phys. Chem. Ref. Data 27, 31–44. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [PubMed] [Google Scholar]
  • 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. [Google Scholar]

All Tables

Table 1.

Chemical reactions.

In the text
Table 2.

Chemical reactions.

In the text
Table 3.

Mineral, precipitation and dissolution parameters.

In the text
Table 4.

Input parameters for each ion.

In the text
Table 5.

B-dot model parameters from the EQ3/6 database [66].

In the text
Table 6.

Chemical reactions.

In the text
Table 7.

Mineral, precipitation and dissolution parameters.

In the text
Table 8.

Input parameters for each ion.

In the text
Table 9.

Physical parameters for the test case of CO2 injection.

In the text
Table 10.

Final carbon distribution at 2000 years.

In the text
Table 11.

CPU time (s) for the DSA and the SIA.

In the text

All Figures

thumbnail Fig. 1.

Discretization by the cell-centered finite volume method.
In the text
thumbnail Fig. 2.

Evolution of the molality of the different components versus time.
In the text
thumbnail Fig. 3.

Evolution of the number of moles of dissolved calcite.
In the text
thumbnail Fig. 4.

Profiles of mineral concentrations.
In the text
thumbnail Fig. 5.

Profiles of pH, aqueous CO2 molality and gas saturation.
In the text
thumbnail Fig. 6.

Mineral net molar changes and evolution of CO2 in time.
In the text
thumbnail Fig. 7.

CPU time and strong parallel efficiency as a function of the number of processors.
In the text
thumbnail Fig. 8.

CPU time and weak parallel efficiency as a function of the number of processors.
In the text
thumbnail Fig. 9.

Comparison of the molalities of H+ and Ca2+ obtained with the DSA and the SIA.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.