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

. In this paper, we consider a parallel ﬁnite volume algorithm for modeling complex processes in porous media that include multiphase ﬂow and geochemical interactions. Coupled ﬂow and reactive transport phenomena often occur in a wide range of subsurface systems such as hydrocarbon reservoir production, groundwater management, carbon dioxide sequestration, nuclear waste repository or geothermal energy production. This work aims to develop and implement a parallel code coupling approach for non-isothermal multiphase multicomponent ﬂow and reactive transport simulation in the framework of the parallel open-source 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 ﬁrstly a multiphase compositional ﬂow problem and then a Direct Substitution Approach (DSA) is used to solve the reactive transport problem. Both subsystems are discretized by a fully implicit cell-centred ﬁnite volume scheme and then an efﬁcient sequential coupling has been implemented in DuMu X . We focus on the stability and robustness of the coupling process and the numerical beneﬁts of the DSA approach. Parallelization is carried out using the DUNE parallel library package based on MPI providing high parallel efﬁciency and allowing simulations with several tens of millions of degrees of freedom to be carried out, ideal for large-scale ﬁeld applications involving multicomponent chemistry. As we deal with complex codes, we have tested and demonstrated the correctness of the implemented software by benchmarking, including the MoMaS reactive transport benchmark, and comparison to existing simulations in the literature. The accuracy and effectiveness of the approach is demonstrated through 2D and 3D numerical simulations. Parallel scalability is investigated for 3D simulations with different grid resolutions. Numerical results for long-term fate of injected CO 2 for geological storage are presented. The numerical results have demonstrated that this approach yields physically realistic ﬂow ﬁelds in highly heterogeneous media and showed that this approach performs signiﬁcantly better than the Sequential Iterative Approach (SIA).


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/gas-liquid). Under certain simplifying conditions, single-phase 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 two-phase 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 single-phase 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 build-up 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 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][16][17][18][19][20][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 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 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][25][26][27][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][30][31][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][36][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][49][50][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 twophase 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 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 long-term 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 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.

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

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 I 2u and I rt . I 2u contains all the species present in both phases while I rt consists of the remaining ones such that I ¼ I 2u [ I rt . So N c ¼ cardfIg is the total number of chemical components involved in N r ¼ cardfI s g chemical reactions (N r ¼ N e þ N k ), with N e ¼ cardfI se g is the number of reactions at equilibrium and N k ¼ cardfI sk g is the number of kinetic reactions. The chemical system writes: where m ij is the stoichiometric coefficient of the species A j in the reaction i.
As an illustration example, we present in Table 1 a chemical system used in [16] in a scenario of CO 2 injection in a deep saline aquifer.
Ano, Cal and Kao refer respectively to Anorthite, Calcite and Kaolinite. In this case, N c = 14 and N r ¼ 8. The first five reactions are at equilibrium while the last three ones are kinetic. Below we give the various corresponding sets above-mentioned: 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: where a j a is the activity of species j in its phase a, K j is the equilibrium constant of reaction j. The activity of water and solid species are set equal to 1. The activity of aqueous species is often written in terms of molality: where c j l is the activity coefficient for species j in the aqueous phase, m j l is the molality of species j [mol kg À1 ] in the aqueous phase, and m 0 is standard molality often taken as 1 mol kg À1 . For each aqueous species j, the molality m j l and the mole fraction of species j in aqueous phase x j l are related by m j l ¼ M H 2 O is the molar mass of water and x H 2 O l is the molar fraction of water in the aqueous phase. The activity of gaseous species is often represented by fugacity: where u j g is the fugacity coefficient for species j in the gas phase, p j g the partial pressure of gas j and p 0 the standard pressure.
Each kinetic reaction leads to an ordinary differential equation: where c i s denotes the concentration of mineral i and r i is the kinetic reaction rate depending on the activities of the species present in the reaction. Expression for r i can be found for instance in [40].

Mathematical model for two-phase multicomponent flow with reactive transport
In the sequel, the index a 2 fl; g; sg (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: where h a [-] denotes the volumetric content of phase a (h a ¼ /S a , / [-] being the porosity of the medium and S a [-] the saturation of phase a if a 2 fl; gg and is the molar density of phase a, x i a is the molar fraction of species i in phase a, D a [m 2 s À1 ] denotes the diffusivity of the phase a: For the sake of simplicity, the diffusion coefficient is assumed to be independent of the chemical species i and depends only on the phase a. D m [m 2 s À1 ] is the molecular diffusion, d L [m] and d T [m] are the magnitudes of longitudinal and transverse dispersion respectively andq a [m s À1 ] is the Darcy velocity of phase a, expressed as follows:q where k ra [-] denotes the relative permeability of phase a, l a [Pa s] is the dynamic viscosity of phase a, K [m 2 ] is the absolute permeability tensor, P a [Pa] is the pressure of phase a, q a [kg m À3 ] is the mass density of phase a and g [m s À2 ] is the gravitational acceleration. The phase pressures are connected by the capillary pressure law: Finally, r j [mol m À3 s À1 ] is the rate of reaction j (it can be equilibrium or kinetic).
We choose to reformulate equations (3) and (4) in term of molar concentration c i a ¼ q mol; a x i a . By neglecting the gradient of molar density rq mol; a resulting from the diffusive flux, equations (3) and (4) For sake of simplicity, we introduce the advection-diffusion operator: L a ðcÞ ¼ r Á ðcq a Þ À r Á ðD a rcÞ: ð9Þ Reactions We can remark that, since D s ¼ 0 andq s ¼0, so that L s ð:Þ ¼ 0.
In order to eliminate the reaction rates r j in equations (7), we make linear combinations between equations (8) with each equation (7). This introduces N p new conservation laws that write: To retrieve the same number of equations as there are unknowns, the N s equations (8) are replaced by N e mass actions laws defined by (1) corresponding to the equilibrium reactions and N k ordinary differential equations corresponding to the kinetic reactions given by (2).

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.

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: Equation (11) corresponds to equation (10) for the species present in liquid and gas phases where i and " i represent respectively the species in the liquid and the gas phase that are in chemical equilibrium and the phase a denotes the phase containing the species j. Equation (12) is the solubility law where a i l and a where K s i and A s i are respectively the kinetic-rate constant [mol m À2 s À1 ] and the reactive surface area [m 2 m À3 ] of mineral i.
Equation (13) corresponds to equation (10) for the primary species present only in one phase whereas equation (14) is the mass action law involved by equilibrium reactions. Finally, equation (15) corresponds to the ordinary differential equations involved by kinetic reactions.
Subsystem (11)- (12) and (13)- (15) are solved sequentially. In equation (13), the Darcy velocities, the phase saturations, the molar densities and the activities or species i 2 I 2u are provided after the computation of the subsystem (11)- (12). In this latter, subsystem (13)-(15) provides explicitly the contribution of species j 2 I s =I 2u as well as a possible update of porosity thanks to the computation of the mineral concentrations.

Finite volume discretization of two-phase compositional flow and reactive transport
The spatial discretization of the coupled system (11)-(15) employs a conservative Finite Volume (FV) method based on a fully upwinding scheme to treat the convective terms and a conforming finite element scheme with piecewise linear elements for the diffusive terms. The time discretization is done by an implicit Euler method.
Here, we choose a cell-centred FV method. It consists in integrating the equations (11)-(15) on a control volume V k (see Fig. 1) and evaluating the fluxes at the interface c kl between two adjacent elements V k and V l .
-Spatial discretization of two-phase compositional flow: where k a ðS l Þ ¼ kraðS l Þ l a is the phase mobility, -Spatial discretization of reactive transport: whereñ kl denotes the unit outer normal to c kl , V ðkÞ is the set of adjacent elements of V k and Á f g k represents the average values on V k . By using the implicit Euler scheme for the time discretization and due to the fact that the primary unknowns (S a , P a , c i a ) and the physical parameters are constant on each element V k , the cell-centred FV schemes corresponding to the discretization of the flow and reactive transport are given by: -Cell-centred FV scheme of two-phase compositional flow: ð24Þ -Cell-centred FV scheme for reactive transport: Now to define the finite volume scheme it is enough to approach the convective and the diffusive fluxes on the interfaces c kl . For this purpose, a fully upwinding scheme is used to calculate the numerical flux for the convective term. More precisely, the quantities (S a , c i a , P a , k a ) are evaluated implicitly and upstream at the interface c kl between two adjacent elements as: The gradient operators on the interfaces c kl are calculated by a P 1 =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 fKg kl and the diffusion coefficients fD a g nþ1 kl at the interface c kl . fq a g nþ1 kl is computed as the arithmetic average of two elements V k and V l . 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 object-oriented software written in C++ and equipped with efficient solvers and massively parallel computation capability.
The first sub-problem consists in solving a convectiondiffusion system modeling a two-phase compositional flow. To tackle this system, we have used a module implemented in DuMu X called 2pnc (two-phase, n-component), where n denotes the number of components present in both phases. The approach is fully implicit. The spatial discretization is performed by the cell-centred finite volume approach described in Section (3.2) by equations (22)- (24). The nonlinear system is solved by a Newton method and a preconditioned BiConjugate Gradient STABilized (BiCGSTAB) method is used to solve the linear system. The control of the time-step is based on the number of iterations required by the Newton method to achieve convergence for the last time iteration. The time-step is reduced, if the number of iterations exceeds a specified threshold, whereas it is increased if the method converges within less iterations. The contribution of the other species is treated explicitly and put as a source term in right hand side of the discretized mass conservation laws (22). We have also implemented the change of porosity linked to the concentrations of minerals.
The second sub-problem consists in solving a singlephase reactive transport problem. In [57], we have developed in the DuMu X framework a one-phase multicomponent transport module named 1pmc-react. In this context, we have implemented a new 1pmc module (one-phase, m-component) that is coupled with the calculations of the chemical problem computed with GSL [61] through a SIA. In this step, the properties of the flow (velocity, saturation, density, concentration of species present in both phases) are given by the two-phase compositional flow. Here, we propose a GIA for this reactive transport sub-problem. More precisely, we used a DSA [42]. The DSA consists in integrating directly the mass action laws (26) in the discretized conservation laws (25). Again, the spatial discretization is carried out using the finite volume method described in Section (3.2) by equations (25)- (27).
A parallel simulator has been developed based on these approximation techniques and is being used to study the sequestration of CO 2 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.

Numerical results
To validate our strategy, several tests have been performed. In [62], we studied the migration of hydrogen produced by the corrosion reaction in deep geological radioactive waste repository and obtained some results coinciding with those presented in [27]. In this contribution, we focus on two tests cases. The first one aims to validate the implementation of the B-dot model [63] for the activity, chemical equilibrium calculations, mineral dissolution and precipitation reaction and the chosen numerical scheme to deal with geochemistry. The second test [16] deals with a scenario of injection of CO 2 in a deep saline aquifer where parallel computations are performed and their performances are detailed.

Validation of geochemistry calculations
We consider a test case presented in [18] which aims to validate the implementation of the activity, chemical equilibrium calculations, mineral dissolution and precipitation reaction and the chosen numerical scheme to deal with geochemistry. There is no flow, no injection and no production. In [18], the authors compared their results with those obtained by the geochemistry software ''The Geochemist's Workbench Ò (GWB)'' [64]. GWB models only the geochemistry but not the transport of components in the porous medium. Consequently, as in [18], we used a one dimensional grid for an example considering only chemistry with no space dependence. The test deals with the dissolution of the calcite at 25°C with the reactions presented in Table 2. Reactions (1)-(6) are in equilibrium while (7) is a kinetic reaction whose parameters are displayed in Table 3.
Initial molalities [mol kg À1 ] are exhibited in Table 4. The activity of each aqueous species is written in terms of molality as: where c 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]: where A, B and _ B are constants depending on temperature, z i is the ion electrical charge for species i and _ A i is the ion size of species i. Finally, I is the ionic strength of the aqueous phase and expresses as I ¼ 0:5 P i m i l z 2 i . Tables 4 and 5 exhibit the parameters for the B-dot model. Initially, the system contains 10 g of calcite and 100 kg of H 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.

Application to geological CO 2 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 CO 2 is implemented according to [65]. Mineral data for the kinetic reactions are summarized in Table 7.
A three-dimensional domain that is 15 km in both the x and y-directions and 100 m in the z-direction is considered. A well perforated in a single grid block located 25 m from the top of the aquifer injects a pure CO 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 two-phase two-component H 2 O-CO 2 flow we have used hydrostatic condition for liquid pressure P l , initial liquid saturation S l ¼ 1 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 B-dot model used as activity model are shown in Table 8.
Impermeable Neumann boundary conditions are enforced on the boundaries of the domain. Constitutive laws and physical parameters are given in Table 9. Figure 4 displays concentrations of calcite, anorthite and kaolinite at 20 and 2000 years for a grid composed of 1.6 · 10 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   1.22 (5) Ca 2+ + Cl À = CaCl + 0.7 (6) HCO À 3 + Ca 2+ = CaCO 3(l) + H + À7.128 (7) Calcite + H + = Ca 2+ + HCO À 3 1.7130 kaolinite are respectively dissolved and precipitated everywhere. 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. 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].

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.

Constitutive law Parameters
Capillary pressure law P c = 0 Pa Absolute permeability     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: where p denotes the number of processors used for the reference time (not always equal to one for heavy computations). For this calculation, we took p = 12, 24 and 48. It points out an optimal use of the parallel resources. Efficiency equal to one indicates that communications and synchronizations between processors are negligible. Figure 7b represents the strong scaling versus the number of processors. A high efficiency (greater than 0.70) is observed up to 256 processors for the computations involving 1.536 · 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. 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:

Weak scaling
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.

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 e SIA . In the sequel, tolerances for the Newton method and iterative algorithm are respectively e Newton = 10 À8 and e 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. 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.

Conclusion
A mathematical formulation and a finite volume approximation were presented for a system of coupled partial differential and algebraic or ordinary differential equations describing two-phase flow with reactive geochemical transport in the subsurface. The proposed reactive transport software module coupling to the two-phase flow was implemented in the framework of the parallel open-source platform DuMu 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 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.