E Ahusborde

B Amaziane

M El Ossmani

M Id Moulay

December

Numerical modeling and simulation of fully coupled processes of reactive multiphase flow in porous media

In this paper, we consider a finite volume approach for modelling multiphase flow coupled to geochemistry in porous media. Reactive multiphase flows are modelled by a highly nonlinear system of degenerate partial differential equations coupled with algebraic and ordinary differential equations. We propose a fully implicit scheme using a direct substitution approach (DSA) implemented in the framework of the parallel open-source platform DuMu X . We focus on the particular case where porosity changes due to mineral dissolution/precipitation are taken into account. This alteration of the porosity can have significant effects on the permeability and the tortuosity. The accuracy and effectiveness of the implementation of permeability/porosity and tortuosity/porosity relationships related to mineral dissolution/precipitation for single-phase and two-phase flows are demonstrated through numerical simulations.

Introduction

Reactive transport modelling is involved in many applications related to subsurface energy and environmental issues. We can mention, no exhaustively, the geological sequestration of CO 2 in saline aquifers, the management of nuclear waste, enhancement of oil recovery, groundwater remediation or deep geothermal energy. A detailed description of these applications and the numerical codes dedicated to reactive transport modelling can be found for instance in [START_REF] Zhang | Groundwater reactive transport models[END_REF], [2] or [START_REF]Reactive Transport Modeling: Applications in Subsurface Energy and Environmental Problems[END_REF].

Equations governing such phenomena consists of a set of nonlinear degenerate system of advectiondiffusion partial differential equations (related to the flow) coupled to algebraic relations and ordinary differential equations (related to the chemistry). In the literature, several strategies can be considered to deal with this set of coupled equations. In the pioneering contribution [START_REF] Yeh | A model for simulating transport of reactive multispecies components: Model development and demonstration[END_REF], the authors present several implicit and sequential approaches for solving reactive transport problems. Global implicit approaches (GIA) tackle the full system of equations while sequential approaches decouple the flow and the chemistry. Due to the strong coupling of the flow and reactive transport equations, a standard approach is to use a GIA to ensure stability in the solution. Although this guarantees numerical stability of the solution, it does not guarantee a nonlinear convergence. This complexity makes difficult the analysis of the entire nonlinear problem. A separation of the different physics can improve the understanding and result in a better design of nonlinear solvers for the reactive transport problem. It is why sequential approaches are more widespread than GIA. Moreover, sequential approaches are easier for implementation since existing codes and specific methods can be used for each sub-problem (flow, transport, chemistry). In [START_REF] Ahusborde | Numerical simulation of two-phase multi-component flow with reactive transport in porous media: application to geological storage of CO 2[END_REF][START_REF] Ahusborde | A sequential approach for numerical simulation of two-phase multicomponent flow with reactive transport in porous media[END_REF][START_REF] Ahusborde | Improvement of numerical approximation of coupled two-phase multicomponent flow with reactive geochemical transport in porous media[END_REF], we developed and integrated in the DuMu X framework a sequential approach. DuMu X [START_REF] Dumu | DUNE for Multi-{Phase, Component, Scale, Physics[END_REF][START_REF] Flemisch | DuMu X : DUNE for Multi-{Phase, Component, Scale, Physics[END_REF] is a free and open-source parallel simulator for flow and transport processes in porous media, based on the Distributed and Unified Numerics Environment DUNE [START_REF] Dune | the Distributed and Unified Numerics Environment[END_REF]. Our strategy splits the global problem into two subproblems. The first sub-problem computes a two-phase compositional flow where only species present in both phases are treated implicitly. 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's velocity for each phase, saturation of each phase, temperature, density,...) are given by the first step. Nonetheless, sequential approaches introduce operator splitting errors [START_REF] Barry | Temporal discretisation errors in non-iterative split-operator approaches to solving chemical reaction/groundwater transport models[END_REF][START_REF] Valocchi | Accuracy of operator splitting for advection-dispersionreaction problems[END_REF] and restrictions on the time-step are mandatory to ensure mass conservation. In [START_REF] Yeh | A model for simulating transport of reactive multispecies components: Model development and demonstration[END_REF], the authors described 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. So, to improve the robustness of the scheme and the possible accuracy loss due to the time-splitting involved by sequential approaches, we switched to a GIA for a single-phase multicomponent flow with reactive transport in [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF]. Our strategy has been validated by numerous numerical examples including 2D and 3D simulations and parallel calculations.

In this work, we propose to extend GIA developed in [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF] to deal with reactive two-phase flows and consequently to drop out the sequential approach considered in [START_REF] Ahusborde | Numerical simulation of two-phase multi-component flow with reactive transport in porous media: application to geological storage of CO 2[END_REF][START_REF] Ahusborde | A sequential approach for numerical simulation of two-phase multicomponent flow with reactive transport in porous media[END_REF][START_REF] Ahusborde | Improvement of numerical approximation of coupled two-phase multicomponent flow with reactive geochemical transport in porous media[END_REF]. Moreover, we are interested in expanding the range of possible applications by considering some examples where porosity and permeability changes must be taken into account. Indeed, the simulation of permeability and tortuosity evolution due to porosity changes can be of crucial importance in the simulation of several processes. These porosity changes can occur due to the dissolution or precipitation of minerals. If the porosity is increased, new pathways can develop, facilitating solute transport while the decrease of porosity can lead to a total clogging, with a possible annihilation of any flow and/or solute transport (see for instance [START_REF] Cochepin | Approaches to modelling coupled flow and reaction in a 2D cementation experiment[END_REF], [START_REF] Xie | Implementation and evaluation of permeability-porosity and tortuosity-porosity relationships linked to mineral dissolution-precipitation[END_REF] or [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF] where permeability-porosity and tortuosity-porosity changes related to mineral dissolutionprecipitation are studied).

The rest of the paper is organized as follows. In Section 2, we describe the governing equations for two-phase multicomponent flow with reactive transport where permeability-porosity and tortuosityporosity changes involved by mineral dissolution-precipitation are taken into account. In Section 3, we present some numerical results. More precisely, a brief description of the numerical strategy is proposed and our method is validated by two test cases performed in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF] including single-phase and two-phase flow simulations. The two-phase simulation is compared in terms of CPU time with a non-reactive two-phase simulation.

Mathematical formulation of the problem

In this section, we present the geochemical and mathematical models for two-phase multicomponent flow with reactive transport in porous media where permeability-porosity changes are considered. For a general discussion on the physical principles we refer e.g. to [START_REF] Helmig | Multiphase Flow and Transport Processes in the Subsurface[END_REF][START_REF] Bethke | Geochemical and biogeochemical reaction modeling: Second edition[END_REF]. We recall here the basic facts and we introduce notation to be used throughout this paper.

Geochemical model

We adopt here the same notation as in [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF]. Precisely, I denotes the set of all the N c chemical components involved in the N r chemical reactions. In this contribution, only equilibrium reactions are considered. Following the Morel formalism, these components are split into primary and secondary components that are noted respectively I p and I s such that I = I p ∪ I s . The set of primary components I p is then divided into mobile primary components I pm and immobile primary components I pi such that I p = I pm ∪ I pi . In the same way, the set of secondary components I s is decomposed into mobile secondary components I sm , immobile secondary components I si and components involved in equilibrium dissolution/precipitation reactions I spe (I s = I sm ∪ I si ∪ I spe ). The chemical system can be written as:

N c ∑ j=1 ν i j A j = 0, i = 1, . . . , N r ,
where ν i j is the stoichiometric coefficient of the component A j in the reaction i.

Each equilibrium reaction gives rise to an algebraic relation called mass action law, relating the activities of the components involved in the reaction:

a j α = K j ∏ i∈I p (a i α ) ν ji , j ∈ I sm ∪ I si , (1)
where a j α is the activity of component j in its phase α, K j is the equilibrium constant of reaction j. The activity of water and solid species are set equal to 1. The activity of aqueous species is often written in terms of molality: a j l = γ j l m j l m 0 , where γ 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 =

x j l

M H2O x H2O l
, where 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:

a j g = ϕ j g P j g
P 0 , where ϕ 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 (1 atm). In the sequel, we consider ideal models for liquid activity (γ j l = 1) and gas fugacity (ϕ j g =1).

For each solid species involved in an equilibrium dissolution/precipitation reaction, a solubility product must be respected:

if K j ∏ i∈I p (a i α ) ν ji < 1 then c j s = 0, else K j ∏ i∈I p (a i α ) ν ji = 1, j ∈ I spe , (2)
where c j s denotes the molar concentration of solid species j [mol.m -3 ]. This complementarity problem is often reformulated as:

min c j s , 1 -K j ∏ i∈I p (a i α ) ν ji = 0, (3)
or using for instance the Fischer-Burmeister complementarity function [START_REF] Fischer | A special Newton-type optimization method[END_REF]:

(c j s ) 2 + 1 -K j ∏ i∈I p (a i α ) ν ji 2 -c j s -1 -K j ∏ i∈I p (a i α ) ν ji = 0. ( 4
)

Mathematical model for two-phase multicomponent flow with reactive transport

In the sequel, the index α ∈ {l, g, s} (l for liquid, g for gas and s for solid) refers to the phase, while the superscript i refers to the component. We define the phase-species correspondence by setting α i to the index of the phase that contains species i. Using the convention used in [START_REF] Farshidi | Chemical reaction modeling in a compositional reservoir-simulation framework[END_REF], the general mass conservation equations for fluid species (α ∈ {l, g}) and solid species (α = s) write:

∂ ∂t (φ S α i ρ mol, α i x i α i ) + ∇ • (ρ mol, α i x i α i q α i ) -∇ • (ρ mol, α i D α i ∇x i α i ) = ∑ j∈I s ν ji r j , i ∈ I pm , (5)
∂ c i s ∂t = ∑ j∈I s ν ji r j , i ∈ I pi , (6)
∂ ∂t (φ S α i ρ mol, α i x i α i ) + ∇ • (ρ mol, α i x i α i q α i ) -∇ • (ρ mol, α i D α i ∇x i α i ) = -r i , i ∈ I sm , (7)
∂ c i s ∂t = -r i , i ∈ I si ∪ I spe , (8)
where φ [-] is the porosity of the medium, S α i [-] denotes the saturation of fluid phase α i , ρ mol, α i [mol.m -3 ] is the molar density of fluid phase α i , x i α i is the molar fraction of fluid species i in the phase

α i , q α i [m.s -1 ] is the Darcy-Muskat velocity of fluid phase α i , D α i [m 2 .s -1
] denotes the diffusivity of the fluid phase α i and c i s is the molar concentration of solid species i.

The Darcy-Muskat velocity of the fluid phase α is expressed as follows:

q α = - k rα (S l ) µ α K(∇P α -ρ α g), (9)
where k rα [-] denotes the relative permeability of fluid phase α, µ α [Pa.s] is the dynamic viscosity of fluid phase α, K [m 2 ] is the absolute permeability tensor, P α [Pa] is the pressure of fluid phase α, ρ α [kg.m -3 ] is the mass density of the fluid phase α and g [m.s -2 ] is the gravitational acceleration.

The phase pressures are connected by the capillary pressure law:

P c (S l ) = P g -P l . (10)
Finally, r j [mol.m -3 .s -1 ] is the rate of equilibrium reaction j.

We choose to reformulate equations ( 5) and ( 7) in terms of molar concentration for fluid species c i α = ρ mol, α x i α (α ∈ {l, g}). By neglecting the gradient of molar density ∇ρ mol, α i resulting from the diffusive flux, equations ( 5) -( 8) can be written as follows:

∂ ∂t (φ S α i c i α i ) + ∇ • (c i α i q α i ) -∇ • (D α i ∇c i α i ) = ∑ j∈I s ν ji r j , i ∈ I pm , (11)
∂ c i s ∂t = ∑ j∈I s ν ji r j , i ∈ I pi , (12)
∂ ∂t (φ S α i c i α i ) + ∇ • (c i α i q α i ) -∇ • (D α i ∇c i α i ) = -r i , i ∈ I sm , (13)
∂ c i s ∂t = -r i , i ∈ I si ∪ I spe . (14)
For the sake of simplicity, we introduce the advection-diffusion operator L α i given by:

L α i (c i α i ) = ∇ • (c i α i q α i ) -∇ • (D α i ∇c i α i ). (15)
Our strategy consists in eliminating the equilibrium reaction rates r j in equations ( 11)- [START_REF] Valocchi | Accuracy of operator splitting for advection-dispersionreaction problems[END_REF]. For this purpose, we make linear combinations between equations ( 13)-( 14) with each equation ( 11)- [START_REF] Valocchi | Accuracy of operator splitting for advection-dispersionreaction problems[END_REF]. This introduces N p new conservation laws written as:

∂ ∂t φ S α i c i α i + ∑ j∈I sm ν ji φ S α j c j α j + ∑ j∈I si ∪I spe ν ji c j s + L α i (c i α i ) + ∑ j∈I sm ν ji L α j (c j α j ) = 0, i ∈ I pm , (16)
∂ dt c i s + ∑ j∈I si ∪I spe ν ji c j s = 0, i ∈ I pi . ( 17
)
To retrieve the same number of equations as there are unknowns, the N s equations ( 13)-( 14) are replaced by card{I sm ∪ I si } mass actions laws defined by ( 1) and card{I spe } complementarity problems defined by ( 3) or ( 4) corresponding to the dissolution/precipitation equilibrium reactions.

Formulation of the global implicit approach

We adopt a global implicit approach using a direct substitution method as in [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF], where only single phase multicomponent flow was studied. To deal with equilibrium dissolution/precipitation reactions, we reformulate the complementarity problem using expression (3). The full system of equations describing the problem writes:

∂ ∂t φ S α i c i α i + ∑ j∈I sm φ ν ji S α i c j α j + ∑ j∈I si ∪I spe c j s + L α i (c i α i ) + ∑ j∈I sm ν ji L α j (c j α j ) = 0, i ∈ I pm , (18)
∂ ∂t c i s + ∑ j∈I si ∪I spe ν ji c j s = 0, i ∈ I pi , (19)
a j α j = K j ∏ i∈I p (a i α i ) ν ji , j ∈ I sm ∪ I si , (20)
min c j s , 1 -K j ∏ i∈I p (a i α i ) ν ji = 0, j ∈ I spe . (21)
The direct substitution method that is considered here, consists in incorporating mass actions laws (20) in mass balance equations ( 18)- [START_REF] Fischer | A special Newton-type optimization method[END_REF].

Evolution of transport properties

Mineral dissolution/precipitation affects the porosity of the rock matrix. More precisely, the porosity φ depends on the concentrations of the minerals according to the relationship:

φ = 1 - N m ∑ j=1 ν j c j s , (22)
where N m is the number of reactive minerals, ν j and c j s represent respectively the molar volume [m 3 .mol -1 ] and the molar concentration [mol.m -3 ] of reactive solid species j. The diffusion coefficient D α is also subject to modification due to the dissolution/precipitation reactions. This dependence is modelled by empirical relationships. For single phase flow, Archie's law is considered [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]:

D l = D 0 l φ φ 0 n , (23)
where D 0 l and φ 0 are respectively the initial diffusion coefficient in liquid phase and initial porosity while n is Archie's empirical coefficient. For two-phase flow, Millington-Quirk's relationship [START_REF] Millington | Permeability of Porous Solids[END_REF] is used:

D α = D 0 α φ 4/3 S 10/3 α , (24)
where D 0 α is the initial diffusion coefficient in phase α. Finally, Kozeny-Carman's relationship is used to deal with the porosity dependence of the permeability:

K = K 0 1 -φ 0 1 -φ 2 φ φ 0 3 , (25)
where K 0 is the initial permeability.

Let us note that the relationship [START_REF] Millington | Permeability of Porous Solids[END_REF], as well as all the dependencies between the transport properties and porosity are treated implicitly.

Numerical results

In this section, we present a brief description of our numerical strategy. In order to validate our strategy taking into account the evolution of porosity and permeability of the porous medium, we consider two tests cases proposed in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF]. The first one is devoted to a single-phase flow while the second one considers a two-phase flow. Since we considered only one-dimensional cases, all computations were performed on a laptop with Intel Core i7-5600U Processor (2.6 GHz) and 8 GB RAM.

Numerical developments

All our developments have been implemented in DuMu X [START_REF] Dumu | DUNE for Multi-{Phase, Component, Scale, Physics[END_REF][START_REF] Flemisch | DuMu X : DUNE for Multi-{Phase, Component, Scale, Physics[END_REF]. In [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF], we have developed and implemented in the DuMu X framework a GIA to deal with reactive single-phase multicomponent flow.

Then in [START_REF] Ahusborde | Improvement of numerical approximation of coupled two-phase multicomponent flow with reactive geochemical transport in porous media[END_REF], we coupled this fully implicit module with a fully implicit compositional two-phase flow through a sequential algorithm. Here, we propose to drop out the sequential strategy and to extend the GIA considered in [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF] to the full system of equations ( 18)- [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]. Moreover, we take into account some porosity changes due to dissolution/precipitation process that induce some permeability and tortuosity alterations.

The spatial discretization is performed using a fully coupled fully implicit conservative finite volume method. A fully upwinding scheme is implemented to treat the convective terms and a conforming finite element scheme with piecewise linear elements is used to calculate the diffusive terms, also referred to as the control volume finite element method [START_REF] Chen | On the control volume finite element methods and their applications to multiphase flow[END_REF]. The time discretization is done by an implicit Euler method. A detailed description of the finite-volume discretization can be found in [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF] and [START_REF] Ahusborde | Improvement of numerical approximation of coupled two-phase multicomponent flow with reactive geochemical transport in porous media[END_REF]. The non-linear system is solved by a Newton method and a preconditioned Biconjugate Gradient Stabilized (BiCGSTAB) method is used to solve the linear system. Numerical differentiation techniques are used to approximate the derivatives in the calculation of the Jacobian matrix. 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. Let us mention that the validation of numerical simulation for reactive multiphase multicomponent flow is a tricky task. Some benchmarks exist for the single-phase case. In [START_REF] Ahusborde | A fully implicit finite volume scheme for single phase flow with reactive transport in porous media[END_REF], we have validated our approach on several examples including the reactive transport benchmark of MoMaS [START_REF] Carrayrou | Reactive transport benchmark of MoMaS[END_REF]. For the two-phase configuration, we note 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. Up to now, to validate our approach, we can only compare with results previously published and by performing a convergence study.

Single-phase flow

Description of the test case

This first example is adapted from the single-phase example considered in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF]. The one-dimensional geometry is depicted in Figure 1. The domain consists of three zones with different diffusivity.

C aq aq D aq B

A min A aq 2.5 cm 2.5 cm 5 cm

x Constant concentration Constant pressure

Constant pressure Constant concentration

Center Right Left

Figure 1: Illustration of the geometry for the single-phase example.

The chemical system and the characteristics of the three different minerals are given in Table 1. Initial physical parameters are presented in Table 2.

Reactions log K Molar volume ν [m 3 .mol -1 ] A min ----A aq + 2 H 2 O 1 0.068 × 10 -3 BC min ----B aq + C aq + H 2 O 2 0.188 × 10 -3 AD min ----D aq + A aq 3 0.218 × 10 -3 OH -----H 2 O -H + 14 -
Table 1: Chemical reactions and mineral data for the single-phase example.

Dirichlet boundary conditions are enforced for the liquid pressure (P l = 10 5 ) and for the concentrations (equal to the initial conditions) on the left and right borders. Table 3 exhibits the chemical initial conditions expressed in molality [mol.kg -1 ] for the aqueous components and in molar concentration [mol.m -3 ] for the minerals. The time of simulation is 10 years.

The goal of this test is to compute the evolution of porosity due to mineral dissolution/precipitation. This changes in porosity and permeability will induce a flow in an problem that is initially purely diffusive. Aqueous component B aq and C aq will react to form mineral BC min , consuming water. Mineral 3: Chemical initial conditions for the single-phase example.

Results and analysis

Firstly, a numerical convergence analysis has been performed. Several one-dimensional meshes have been considered for this test case. Figure 2 represents the profiles for the porosity and the concentration of the mineral AD min at the end of the simulation for several grid resolutions. As expected, the quantities converge toward a reference solution when the grid resolution increases.

Figure 3 exhibits the evolution of the porosity and the concentrations of the three minerals for a total of 60 nodes as in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF]. As expected, from the beginning of the simulation, mineral A min is dissolved at both ends of the domain where a constant pressure is enforced. The porosity increases, reducing the pressure and generating a flow rate whose evolutions are depicted in Figure 4. Meanwhile, on the left part of the domain, A aq is released by the dissolution of A min and it reacts with D aq that is initially present to precipitate AD min . The precipitation of BC min occurs later, after the diffusion of B aq and C aq that are represented in Figure 5. This precipitation leads to a pressure build-up in the centre of the domain. These results are in good agreement with those obtained in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF] where a detailed description and interpretation of the results are provided.

Finally, Figure 6 represents the time-step used during the computations and the number of iterations in Newton's algorithm for a mesh composed of 1000 elements. We have to specify that a maximum time-step equal to 10 6 s was enforced. As expected, the time-step is reduced when a high number of 0,02 0,03 0,04 0,05 0,06 0,07 0,08

x-position (m)

0 0,5 1 1,5 2 2,5
A min (mol/l) 0,02 0,03 0,04 0,05 0,06 0,07 0,08

x-position (m) 0 0,1 iterations is required to reach convergence in Newton's algorithm.

Two-phase flow 3.3.1 Description of the test case

The second test deals with a two-phase reactive flow. It is adapted from the two-phase example presented in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF]. The chemical system is depicted in Table 4. In comparison with the previous example, a gas phase consisting of N 2(g) and D g is present. N 2(aq) is in equilibrium with N 2(g) and the molar fractions of N 2(g) in gas phase and the activity of N 2(aq) in liquid phase are related by Henry's law:

x N 2(g) g P g P 0 = K N 2 a N 2(aq) l
.

The gaseous component D g will be dissolved in water to form D -and H + . The molar fractions of D g in gas phase and the activity of D -in liquid phase are related by a mass action law that writes:

x D g g P g P 0 = K D a D - l a H + l .
D -will react with A aq to form AD min . Figure 7 represents the geometry of the example. Initially, the rock zone (right part) contains the mineral A min while D g is only present in the atmospheric zone (left part). On the left boundary, a constant saturation and a constant partial pressure of gaseous component D g are imposed while no-flow boundary conditions are enforced on the right border. Chemical initial conditions are provided in Table 5 and the other physical parameters are given in Table 6. Finally, Table 7 depicts the parameters for van Genuchten's relationships used for capillary pressure and relative permeabilities laws :

P c (S l ) = 1 α V G S l -1 m -1 1-m , (26)
k rl (S l ) = S l 1 2 1 -(1 -S l 1 m ) m 2 , (27)
k rg (S l ) = (1 -S l ) 1 3 1 -S l 1 m 2m , (28)
with S l = S l -S r l 1 -S r l -S r g . The time of simulation is 1 year.

Reactions log K Molar volume

ν (m 3 .mol -1 ) A min ----A aq + 2 H 2 O 1 0.068 × 10 -3 AD min ----D -+ A aq + H + 10 0.218 × 10 -3 OH -----H 2 O -H + 14 - D g ----D -+ H + -H 2 O 5 - N 2(g) ----N 2(aq)
3.189 -Table 4: Chemical reactions and mineral data for the two-phase example.

x Table 6: Initial hydrodynamic parameters for the two-phase example.

Results and analysis

As for the single-phase example, a numerical convergence analysis has been performed using several grid resolutions. Figure 8 represents the profiles for the porosity and the concentration of the mineral AD min at the end of the simulation as a function of the number of elements composing the meshes. We still observe the convergence for both quantities toward a reference solution when refining the mesh. For a mesh composed of 100 nodes, Figure 9 mineral A min is dissolved while mineral AD min precipitates, reducing the porosity close to the interface. Similar results can be observed in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF]. Table 8 compares the CPU time for the different meshes with the CPU time consumed for a simple two-phase flow. This highlights the complexity of the test case since even for one dimensional simulations, the computation time is far from negligible in comparison with the simple two-phase case. This is due to the strong non-linearities generated by the geochemical reactions and to the significant difference in sizes of the two problems (9 unknowns per element for the reactive case versus 2 unknowns per mesh for the non-reactive case). As for the single-phase simulation, Figure 10 displays the time-step used during the computations and the number of iterations in Newton's algorithm that are still strongly correlated. In this case, a mesh composed of 500 elements and a maximum time-step equal to 10 4 s have been considered.

Conclusion

We proposed a fully coupled, fully implicit finite volume approach for solving a system of coupled partial differential and algebraic or ordinary differential equations describing multiphase flow with reactive 14 geochemical transport in the subsurface. All the developments have been implemented in the framework of the parallel open-source platform DuMu X . To validate our strategy, we focused on examples where changes of permeability and tortuosity induced by porosity alteration are taken into account. Indeed, in many applications, it can be crucial to consider the feedback between multiphase flow and reactive transport: mineral dissolution will enhance porosity while precipitation can lead to a complete clogging. Single-phase and two-phase numerical examples have been performed. In both cases, numerical convergence analyses were carried out and numerical results are in good agreement with those obtained in [START_REF] Seigneur | Recoupling flow and chemistry in variably saturated reactive transport modelling -An algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow[END_REF]. They provided validation of our implementation of permeability/porosity and tortuosity/porosity relationships related to mineral dissolution/precipitation. Next, we are planning to validate our fully implicit approach for reactive multiphase flows through more complicated cases including complex 3D geometries or/and complex chemical systems. An advanced comparison with the sequential approach developed in [START_REF] Ahusborde | Improvement of numerical approximation of coupled two-phase multicomponent flow with reactive geochemical transport in porous media[END_REF] is in progress.

Figure 2 :

 2 Figure 2: Porosity and AD min profiles at t = 10 years for different grid resolutions (note that only an area of interest between x=0.02 and x=0.08 is presented).

Figure 3 :

 3 Figure 3: Evolution of porosity and mineral profiles (note that only an area of interest between x=0.02 and x=0.08 is presented).

Figure 4 :

 4 Figure 4: Left: Generated flow rate versus time (The arrows indicate the direction of the flow). Right: Pressure versus time.

Figure 5 :

 5 Figure 5: Evolution of the aqueous concentrations.

Figure 6 :

 6 Figure 6: Time step and number of iterations in Newton's algorithm versus time for the single-phase example.

10 - 4

 104 depicts the evolution of the saturation that illustrates the drying of the material. Moreover it shows the evolution of the different minerals and the porosity: Pa-1

Figure 8 :

 8 Figure 8: Porosity and AD min profiles at t = 1year for different meshes (only an area of interest close to interface is presented).

Figure 9 :

 9 Figure 9: Evolution of saturation, porosity and mineral profiles for the two-phase simulation (note that except for saturation, only an area of interest close to the interface is presented).

Figure 10 :

 10 Figure 10: Time-step and number of iterations in Newton's algorithm versus time for the two-phase example.

Table 2 :

 2 Initial Hydrodynamic parameters for the single-phase example.will dissolve, releasing A aq that will react with D aq to precipitate AD min .

	A min

Table

Table 5 :

 5 Chemical initial conditions for the two-phase example.

	Constant concentration	Constant pressure/saturation

Table 7 :

 7 Parameters for the van Genuchten relationship.

		0,52												0,2
		0,5												0,18									100 elements
																							200 elements
	Porosity (%)	0,42 0,44 0,46 0,48									1000 elements 500 elements 100 elements 200 elements	AD min (mol/l)	0,08 0,1 0,12 0,16 0,14									500 elements 1000 elements
		0,4
														0,06
		0,38												0,04
		0,36												0,02
		0,34	0	0,01	0,02	0,03	0,04	0,05	0,06	0,07	0,08	0,09	0,1	0	0	0,02	0,04	0,06	0,08	0,1	0,12	0,14	0,16	0,18	0,2
							x-position (m)										x-position (m)

Table 8 :

 8 CPU time (s) as a function of the grid resolution for reactive and non-reactive two-phase flows.

	Number of elements CPU time (s)	CPU time (s)
		for reactive	for non-reactive
		two-phase flow two-phase flow
	100	165.2	19.3
	200	379.4	37.8
	500	1187.7	92.6
	1000	4923.1	175.4

Acknowledgments

This work was partially supported by the Carnot Institute, ISIFoR project (Institute for the sustainable engineering of fossil resources). This support is gratefully acknowledged. We also thank the DuMu X and DUNE teams for their help during the development of our reactive transport module and Nicolas Seigneur for his helpful indications and fruitful discussions for the numerical tests. Finally, the authors gratefully thank the anonymous referees for their insightful comments and suggestions, as these comments led us to an improvement of the work.