Regular Article
Improvement of numerical approximation of coupled multiphase multicomponent flow with reactive geochemical transport in porous media
^{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
^{2}
University Moulay Ismaïl, EMMACSENSAM, Marjane II, 50000 Meknès, Morocco
^{*} Corresponding author: brahim.amaziane@univpau.fr
Received:
20
February
2018
Accepted:
6
July
2018
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 nonisothermal multiphase multicomponent flow and reactive transport simulation in the framework of the parallel opensource platform DuMu^{X}. 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 cellcentred finite volume scheme and then an efficient sequential coupling has been implemented in DuMu^{X}. 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 largescale 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 longterm fate of injected CO_{2} 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).
© E. Ahusborde et al., published by IFP Energies nouvelles, 2018
This 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 CO_{2} 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 CO_{2}: 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 twophase flow (supercritical/gasliquid). Under certain simplifying conditions, singlephase flow is considered for instance in [8] and [9]. In [8], in the framework of SHPCO_{2} Project [10], the gas phase is assumed to be immobile and therefore gaseous carbon dioxide is considered as a fixed species neglecting the twophase flow effects. In [9], an initial amount of supercritical CO_{2} is converted into a source term of liquid CO_{2} and then the authors study the transport of the dissolved CO_{2} and the precipitation/dissolution process of minerals. In [11], the authors employ a singlephase reactive flow to model the leaking of CO_{2}saturated brine in a fractured pathway once supercritical CO_{2} is totally dissolved. CO_{2} is generally injected in its supercritical form. This injection may induce important pressure buildup that can damage the reservoir or induce fracturing and seismic events. Moreover, the supercritical CO_{2} 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_{2} to circumvent the abovementioned 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 onephase flow models can predict pressure buildup far from the injection as well as complex twophase flow models, most of the studies deal with twophase reactive flows. The following is a nonexhaustive list of references: [15–21], the main difference between theses references being the complexity of the geochemical system.
During its storage, liquid or supercritical CO_{2} 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 longterm safety of the disposal of nuclear waste is an important issue in all countries with a significant nuclear program. Repositories for the disposal of highlevel and longlived radioactive waste generally rely on a multibarrier system to isolate the waste from the biosphere. The multibarrier 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_{2}, 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 multibarrier 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, DarcyMuskat’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 NonIterative (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 operatorsplitting 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 onedimensional investigations” due to their complexity and their high computational requirements. Thanks to the advance of highperformance 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 twophase 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 interphase 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 subproblems. The first subproblem computes an implicit twophase 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 subproblem 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 subproblem. To improve the robustness of the scheme and the accuracy loss due to the timesplitting 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 twophase multicomponent flow with reactive transport. Our methodology consisting in decoupling the original problem is explained in Section 3 where the finitevolume discretization of each subproblem is also detailed. In Section 4, a description of the implementation of our strategy in the DuMu^{X} 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 longterm fate of injected CO_{2} for geological sequestration including parallel scaling results. For this test, a comparison between the DSA and the SIA for solving the reactive transport subproblem 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 twophase 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} ∪ 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} ∪ I_{sk}).
In order to use a sequential approach to decouple twophase flow and reactive transport, we introduce two subsets and . contains all the species present in both phases while consists of the remaining ones such that . So is the total number of chemical components involved in chemical reactions (), with is the number of reactions at equilibrium and is the number of kinetic reactions. The chemical system writes:where is the stoichiometric coefficient of the species in the reaction .
As an illustration example, we present in Table 1 a chemical system used in [16] in a scenario of CO_{2} injection in a deep saline aquifer.
Chemical reactions.
Ano, Cal and Kao refer respectively to Anorthite, Calcite and Kaolinite. In this case, N_{c} = 14 and . The first five reactions are at equilibrium while the last three ones are kinetic. Below we give the various corresponding sets abovementioned:
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:(1)where is the activity of species in its phase , is the equilibrium constant of reaction . The activity of water and solid species are set equal to . The activity of aqueous species is often written in terms of molality: , where is the activity coefficient for species in the aqueous phase, is the molality of species [mol kg^{−1}] in the aqueous phase, and is standard molality often taken as 1 mol kg^{−1}. For each aqueous species , the molality and the mole fraction of species in aqueous phase are related by , where is the molar mass of water and is the molar fraction of water in the aqueous phase. The activity of gaseous species is often represented by fugacity: , where is the fugacity coefficient for species in the gas phase, the partial pressure of gas and the standard pressure.
Each kinetic reaction leads to an ordinary differential equation:(2)where denotes the concentration of mineral and is the kinetic reaction rate depending on the activities of the species present in the reaction. Expression for can be found for instance in [40].
2.2 Mathematical model for twophase multicomponent flow with reactive transport
In the sequel, the index (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:(3)(4)where [–] denotes the volumetric content of phase (, ϕ [–] being the porosity of the medium and [–] the saturation of phase if and ), [mol m^{−3}] is the molar density of phase , is the molar fraction of species in phase , [m^{2} s^{−1}] denotes the diffusivity of the phase α:
For the sake of simplicity, the diffusion coefficient is assumed to be independent of the chemical species i and depends only on the phase . [m^{2} s^{−1}] is the molecular diffusion, d_{L} [m] and d_{T} [m] are the magnitudes of longitudinal and transverse dispersion respectively and [m s^{−1}] is the Darcy velocity of phase , expressed as follows:(5)where [–] denotes the relative permeability of phase , [Pa s] is the dynamic viscosity of phase , [m^{2}] is the absolute permeability tensor, P_{α} [Pa] is the pressure of phase α, [kg m^{−3}] is the mass density of phase α and [m s^{−2}] is the gravitational acceleration.
The phase pressures are connected by the capillary pressure law:(6)
Finally, [mol m^{−3} s^{−1}] is the rate of reaction (it can be equilibrium or kinetic).
We choose to reformulate equations (3) and (4) in term of molar concentration . By neglecting the gradient of molar density resulting from the diffusive flux, equations (3) and (4) write:(7)(8)
For sake of simplicity, we introduce the advectiondiffusion operator:(9)
We can remark that, since and , so that .
In order to eliminate the reaction rates in equations (7), we make linear combinations between equations (8) with each equation (7). This introduces new conservation laws that write:(10)
To retrieve the same number of equations as there are unknowns, the equations (8) are replaced by mass actions laws defined by (1) corresponding to the equilibrium reactions and 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 subproblems. The first subproblem is dedicated to a twophase 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 subproblem 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 twophase multicomponent flow with reactive transport in porous media that will be used in the sequel.
The system of equations for the twophase compositional flow is given by:(11)(12)
Equation (11) corresponds to equation (10) for the species present in liquid and gas phases where and represent respectively the species in the liquid and the gas phase that are in chemical equilibrium and the phase denotes the phase containing the species . Equation (12) is the solubility law where and 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:(13)(14)(15)where and are respectively the kineticrate constant [mol m^{−2} s^{−1}] and the reactive surface area [m^{2} m^{−3}] of mineral .
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 are provided after the computation of the subsystem (11)–(12). In this latter, subsystem (13)–(15) provides explicitly the contribution of species as well as a possible update of porosity thanks to the computation of the mineral concentrations.
3.2 Finite volume discretization of twophase 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 cellcentred FV method. It consists in integrating the equations (11)–(15) on a control volume (see Fig. 1) and evaluating the fluxes at the interface between two adjacent elements and .
Spatial discretization of twophase compositional flow:
Spatial discretization of reactive transport:
Cellcentred FV scheme of twophase compositional flow:
Cellcentred FV scheme for reactive transport:
Fig. 1. Discretization by the cellcentered finite volume method. 
Now to define the finite volume scheme it is enough to approach the convective and the diffusive fluxes on the interfaces . For this purpose, a fully upwinding scheme is used to calculate the numerical flux for the convective term. More precisely, the quantities (, , , ) are evaluated implicitly and upstream at the interface between two adjacent elements as:(28)
The gradient operators on the interfaces γ_{kl} are calculated by a finite element method with piecewise linear elements. An harmonic average of the values between two adjacent elements is used to calculated the absolute permeability and the diffusion coefficients at the interface γ_{kl}. is computed as the arithmetic average of two elements and .
4 Numerical simulations
4.1 Implementation
The physical models described above are built in DuMu^{X} [58, 59] framework, which handles general input/output, memory management, grid generation, parallelism, etc. The code is an objectoriented software written in C++ and equipped with efficient solvers and massively parallel computation capability.
The first subproblem consists in solving a convectiondiffusion system modeling a twophase compositional flow. To tackle this system, we have used a module implemented in DuMu^{X} called (twophase, ncomponent), where n denotes the number of components present in both phases. The approach is fully implicit. The spatial discretization is performed by the cellcentred 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 timestep is based on the number of iterations required by the Newton method to achieve convergence for the last time iteration. The timestep 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 subproblem consists in solving a singlephase reactive transport problem. In [57], we have developed in the DuMu^{X} framework a onephase multicomponent transport module named 1pmcreact. In this context, we have implemented a new module (onephase, mcomponent) 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 twophase compositional flow. Here, we propose a GIA for this reactive transport subproblem. 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_{2} and gas migration in a nuclear waste repository. Nonisothermal 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 Bdot 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_{2} 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.
Chemical reactions.
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:(29)where γ^{i} is the activity coefficient for species i in the aqueous phase. Several models of activity coefficient can be used. Here, the Bdot model is considered [63]:(30)where A, B and are constants depending on temperature, is the ion electrical charge for species i and is the ion size of species . Finally, is the ionic strength of the aqueous phase and expresses as .
Input parameters for each ion.
Tables 4 and 5 exhibit the parameters for the Bdot model. Initially, the system contains 10 g of calcite and 100 kg of H_{2}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.
Fig. 2. Evolution of the molality of the different components versus time. 
Fig. 3. Evolution of the number of moles of dissolved calcite. 
4.2.2 Application to geological CO_{2} sequestration
We now apply our simulator to a 3D coupled twophase 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 “sixelement 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_{2} is implemented according to [65]. Mineral data for the kinetic reactions are summarized in Table 7.
Chemical reactions.
Mineral, precipitation and dissolution parameters.
A threedimensional domain that is 15 km in both the x and ydirections and 100 m in the zdirection is considered. A well perforated in a single grid block located 25 m from the top of the aquifer injects a pure CO_{2} stream at constant rate during the first 20 years. After the 20 years injection period, a total of 18.6 × 10^{9} kg of CO_{2} is injected.
As initial conditions for the twophase twocomponent H_{2}OCO_{2} flow we have used hydrostatic condition for liquid pressure , initial liquid saturation and initial CO_{2} molality in liquid phase equals 3.55 × 10^{−3} mol kg^{−1}. Initial conditions for the reactive transport problem and parameters for the Bdot model used as activity model are shown in 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.
Physical parameters for the test case of CO_{2} injection.
Figure 4 displays concentrations of calcite, anorthite and kaolinite at 20 and 2000 years for a grid composed of 1.6 × 10^{5} 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 CO_{2} and precipitated far from the injection while anorthite and kaolinite are respectively dissolved and precipitated everywhere.
Fig. 4. Profiles of mineral concentrations. 
Figure 5 depicts the molality of aqueous CO_{2}, the pH and the gas saturation at 20 and 2000 years. The pH and the molality of aqueous CO_{2} are strongly correlated since high concentration of aqueous CO_{2} acidifies the liquid phase. After the injection, the gaseous CO_{2} migrates upward and spreads laterally when reaching the top of the aquifer that is impermeable.
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 CO_{2} shows that most of CO_{2} 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].
Fig. 6. Mineral net molar changes and evolution of CO_{2} in time. 
Fig. 7. CPU time and strong parallel efficiency as a function of the number of processors. 
Fig. 8. CPU time and weak parallel efficiency as a function of the number of processors. 
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” 12Core E52690 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 × 10^{5}, 1.28 × 10^{6} and 5.76 × 10^{6} elements corresponding respectively to 1.92 × 10^{6}, 1.536 × 10^{7} and 6.912 × 10^{7} degrees of freedom (dof). The dashed lines represent an ideal behaviour.
Strong efficiency is given by:(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, and . 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 × 10^{7} and 6.912 × 10^{7} dof. For the simulation with 1.92 × 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.
Weak efficiency is given by:(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 subproblem for the example presented above. Both approaches adopt an adaptive timestepping. In the DSA, the control of the timestep 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 . 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 Ca^{2+} 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.
Fig. 9. Comparison of the molalities of H^{+} and Ca^{2+} 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.
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 twophase flow with reactive geochemical transport in the subsurface. The proposed reactive transport software module coupling to the twophase flow was implemented in the framework of the parallel opensource platform DuMu^{X}. 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 twophase 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^{X} 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 2017A0010610019 and 2018A0040610019 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 CO_{2} in deep saline formations, Springer. [CrossRef] [Google Scholar]
 Zhang F., Yeh G.T., Parker J.C. (2012) Groundwater reactive transport models, Bentham ebooks [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 longterm geological storage of CO_{2}, 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]
 AlKhoury R., Bundschuh J. (2014) Computational models for CO_{2} geosequestration and compressed air energy storage, Sustainable Energy Developments, CRC Press. [CrossRef] [Google Scholar]
 Bear J., Carrera J. (2017) Mathematical modeling of CO_{2} storage in a geological formation, Springer. [Google Scholar]
 Haeberlein F. (2011) Time space domain decomposition methods for reactive transport – application to CO_{2} geological storage, PhD Thesis, Université ParisNord – Paris XIII. [Google Scholar]
 Lagneau V., Pipart A., Catalette H. (2005) Reactive transport modelling of CO_{2} 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 multispecies reactivetransport in heterogeneous porous media applied to CO_{2} geological storage. https://www.ljll.math.upmc.fr/mcparis09/Files/haeberlein_poster.pdf. [Google Scholar]
 Ahmad N., Wörman A., BottacinBusolin A., SanchezVila X. (2015) Reactive transport modeling of leaking CO_{2}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 CO_{2}, Adv. Water Resour. 62, 533–542. [Google Scholar]
 Ahmad N., Wörman A., SanchezVila X., Jarsjö J., BottacinBusolin A., Hellevang H. (2016) Injection of CO_{2} 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 singlephase flow numerical models sufficient to estimate pressure distribution in CO_{2} sequestration projects?, Energy Procedia 4, 3919–3926. [Google Scholar]
 Audigane P., Gaus I., CzernichowskiLauriol I., Pruess K., Xu T. (2007) Twodimensional reactive transport modelling of CO_{2} 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 fullycoupled flowreactivetransport formulation based on element conservation, with application to CO_{2} 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 CO_{2}brinerock 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 CO_{2} storage in aquifers with a fullycoupled 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 CO_{2} 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 CO_{2} mineralization during CO_{2} 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 CO_{2}saturated water, Int. J. Greenh. Gas Con. 4, 561–574. [CrossRef] [Google Scholar]
 Jacquemet N., Pironon J., Lagneau V., SaintMarc J. (2012) Armouring of well cement in H_{2}SCO_{2} 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 lowpH cement liner on the near field of a repository for spent fuel and highlevel radioactive waste, Phys. Chem. Earth 64, 46–56. [CrossRef] [Google Scholar]
 Mon A., Samper J., Montenegro L., Naves A., Fernández J. (2017) Longterm nonisothermal 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) Corrosioninduced 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 thermalhydrochemical 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 nonideal aqueoussolid 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 multisite 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. [CrossRef] [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., CulliganHensley P.J. (1996) Temporal discretisation errors in noniterative splitoperator 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 advectiondispersionreaction 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 globalimplicit 2D 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 MIN3Pmodel 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 transportreaction problems in porous media, Water Resour. Res. 41, 1–17. [Google Scholar]
 Kräutle S., Knabner P. (2007) A reduction scheme for coupled multicomponent transportreaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions, Water Resour. Res. 43. [CrossRef] [PubMed] [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 twophase multicomponent flow with reactive transport in porous media: application to geological sequestration of CO_{2}, ESAIM: Proc. Surveys 50, 21–39. [CrossRef] [Google Scholar]
 Ahusborde E., El Ossmani M. (2017) A sequential approach for numerical simulation of twophase multicomponent flow with reactive transport in porous media, Math. Comput. Simulat. 137, 71–89. [CrossRef] [Google Scholar]
 DuMu^{X} (2018) DUNE for multiPhase, 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) DuMu^{X}: 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 RootFinding. 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 twophase flow with reactive transport in porous media, Finite volumes for complex applications VIIIhyperbolic, 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) CO_{2}H_{2}O mixtures in the geological sequestration of CO_{2}. 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 UCRLMA110662 PT I. [CrossRef] [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 triplepoint 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 leadershipclass 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., NavarreSitchler 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
All Figures
Fig. 1. Discretization by the cellcentered finite volume method. 

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

In the text 
Fig. 3. Evolution of the number of moles of dissolved calcite. 

In the text 
Fig. 4. Profiles of mineral concentrations. 

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

In the text 
Fig. 6. Mineral net molar changes and evolution of CO_{2} in time. 

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

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

In the text 
Fig. 9. Comparison of the molalities of H^{+} and Ca^{2+} obtained with the DSA and the SIA. 

In the text 