Mikhail Panfilov

Alain Bourgeat

Leonid Pankratov

Sergey Shmarev

Tatyana Chechkina

Gregory Chechkin

Je Remercie Également

Caroline De Dieuleveult

Ghassan Jomaa

Anvar Farkhutdinov

Trevor Cox

Marco Campestrini

Dat Tri

Alain Ngo

Gaunand

Irina Sin
email: irina.sin@mines-paristech.fr

2D simulation of natural gas reservoir by multiphase multicomponent reactive flow and transportdescription of a benchmarking exercise

Keywords: Compressible two-phase flow, Reactive transport, Sequential iterative coupling, Operator splitting, HYTEC, EOS

mique ainsi que pour toutes les discussions productives, sa disponibilité et participation enthousiaste au travail, et pour m'avoir ouvert les portes du Centre Thermodynamique des Procédés de l'École des Mines (CTP). J'ai beaucoup apprécié bien la collaboration avec Martha Hajiw et Tianyuan Wang sur les données de CPA-PR. Mes remerciements vont également à Philippe Rivière du Centre Efficacité Énergétique des Systèmes (CES) pour ses conseils et son soutien permanent depuis le premier jour à l'École des Mines. Je remercie Claude Tadonki d'avoir partagé son expérience sur GMRES. Je suis très reconnaissante aux murs et aux pavés de l'École des Mines et en particulier, à tout le Centre de Géosciences à Fontainebleau pour leur accueil cordial. J'ai beaucoup apprécié le travail sur le benchmark avec Laurent de Windt, son initiative et son expérience remarquable en géochimie. Je remercie Raphaël Mathieu d'avoir été mon guide dans le monde d'HYTEC, je le remercie également pour tout le support informatique, toujours immédiat. Mes remerciements vont bien entendu aussi à ma chère Rosa-Lynda, inestimable et indispensable, pour ses efforts administratifs et pour sa générosité d'âme. Je tiens à remercier Anthony Michel pour ses conseils prudents et explications enrichissantes dès le début de ma thèse. J'ai tiré beaucoup d'enseignements de nos échanges sur le couplage et le benchmarking avec lui et Thibault Faney.

Dans les codes de transport réactif, le nombre de composants dans le problème d'écoulement multiphasique est généralement minimisé puisque chaque composant supplémentaire augmente la dimension de système non-linéaire. Pour réduire la complexité numérique, la formulation de conservation des phases couplée au transport compositionnel linéaire et aux équations d'état cubiques s'est employée à résoudre l'écoulement multicomposant multiphasique compressible. Ainsi, plus le système chimique est riche en composants, plus le couplage d'écoulement peut être avantageux en comparaison avec les méthodes basées sur la formulation compositionnelle. La résolution d'écoulement peut être découplée comme un module à l'aide de la méthode proposée et intégrée dans le cadre de codes de transport réactif. Une approche alternative a ensuite été développée pour modéliser séparément le vii viii transport multiphasique en utilisant des termes de taux de réaction numériques. Le système de couplage d'écoulement multiphasique et de transport réactif hérite la propriété de stabilité inconditionnelle du schéma complètement implicite d'écoulement, ce qui fait la modélisation des régimes d'advection dominante efficace.

Un problème de benchmark avec une solution auto-similaire a été modélisé pour vérifier la méthode. Par la suite, une comparaison avec une formulation globale implicite a été effectuée sur un problème de deux composants et deux phases, les résultats ont démontré la performance numérique de la méthode. Pour vérifier et valider les méthodes de couplage de simulateurs multiphasiques réactifs, le benchmarking est indispensable. Nous avons établi un exercice de benchmark basé sur la modélisation d'un réservoir de gaz naturel (avec trois gaz). Le problème suppose une relation étroite entre le transport et la géochimie et vise à examiner la capacité numérique des méthodes de couplage et leurs caractéristiques dans le contexte de l'écoulement et du transport multiphasiques multicomposants, solubilité de gaz et interaction de gaz-eau-roche.

Les équations d'état peuvent être employées globalement ou localement pour des propriétés thermodynamiques différentes au sein du couplage. La calibration de leurs paramètres varie en fonction des besoins. Un aperçu des équations d'état a été donné sous la forme d'une classification par type de modèle et son application à des problèmes différents. Une analyse comparative sur propriétés volumétriques de corps purs et de systèmes binaires a été effectuée.

Les développements numériques de ce travail et leur implémentation dans le code de transport réactif HYTEC ont été utilisés pour modéliser captage et stockage du carbone : en particulier, l'impact de mélange convectif et de présence des impuretés dans le gaz injecté. La méthode de couplage proposée a permis d'étudier des scenarios d'injection de gaz complexes, incluant sans s'y limiter des mélanges de quatre gaz, et de révéler de nouveaux comportements de la dynamique des fluides, gazeux et liquides -interaction entre le flux convectif densitaire et la séparation chromatographique, leurs effets sur la vitesse du front de courant de gaz et sur la dissolution globale de gaz.

Le problème de Hele-Shaw de mélange convectif dans la phase aqueuse en présence du gaz et de la frange capillaire a été modélisé. Les résultats ont indiqué qu'en plus des mouvements convectifs dans la phase aqueuse, les instabilités, qui continuaient à grandir sur la surface de contact irrégulière entre deux phases, ont engendré des perturbations cor-ix respondantes dans la phase gazeuse, favorisant le mélange de gaz. Par la suite, le mélange convectif dans les deux phases a aussi été observé à l'échelle du réservoir.

La modélisation 3D de stockage du carbone montre une dissolution convective relativement faible à l'échelle du réservoir, puisque maillage à haute résolution est critique pour caractériser précisément le flux massique. Nous avons donc examiné la sous-estimation convective numérique sur la géométrie radiale en utilisant l'approximation par morceaux en temps d'une solution analytique, non-convective, couplée à une expression analytique de transfert convectif. La méthode donne une prédiction du rayon de courant en haut de l'aquifère en tenant compte de la dissolution convective. L'étude a confirmé que la modélisation 2D radiale pouvait être employée pour estimer le transfert convectif, en réduisant la dimension de maillage et en améliorant le mélange convectif en même temps.

x

Long abstract

Human activity in the subsurface has been rapidly expanding and diversifying (waste disposal, new mining technologies, high-frequency storage of energy), while the public and regulatory expectations keep growing. The assessment of each step of underground operations requires careful safety and environmental impact evaluations. Multiphysics modeling provides an effective way to understand and predict the behavior of such complex systems at different time and space scales. It is of great importance to develop an efficient method for integrating multiphase flow into the reactive transport framework.

In this work, we proposed an operator-splitting based coupling approach that is devised to model multiphase chemical systems, conserving the versatile structure of an existing reactive transport code, and hence, it retains the general paradigm of the flow module independent of the geochemical processes. The method was implemented in the reactive transport simulator HYTEC, building new capabilities in terms of applications and further developments.

In the reactive transport codes, the number of components in the multiphase flow problem is usually minimized since any additional component increases the dimension of nonlinear system. To reduce the computational complexity, the phase conservation formulation coupled with the linear component transport and cubic equations of state was employed for solving compressible multiphase multicomponent flow. Thus, the more the chemical system is abundant in components, the more the proposed flow coupling can be advantageous in comparison with the methods based on the compositional formulation.

The flow solution can be split as a module, by means of the proposed method, and integrated into the reactive transport framework. Thereafter, an alternative way based on the sequential iterative approach was invented to separately model multiphase transport using numerical reaction rate terms. The coupled system of multiphase flow and reactive transxi xii port inherits the unconditional stability of the fully implicit scheme of flow that allows to efficiently model advective dominant regimes.

The method was verified by modeling a benchmark problem admitting a self-similar solution. Then, a comparison with the global implicit formulation was carried out, that demonstrated the computational efficiency of the method on the two-phase two-component problem. In order to verify and validate coupling approaches of reactive multiphase simulators, benchmarking is essential. Thus, we presented a benchmark exercise with tightly connected transport and geochemistry by modeling a natural gas reservoir (with three gases) that aims at specifically analyzing the numerical capacity of the coupling methods and their characteristics in the context of multiphase multicomponent flow and transport, gas solubility and gas-water-rock interaction.

Within the coupling, the equations of state can be deployed globally or locally for different thermodynamic properties. Their parameters and formulations are calibrated according to the modeled property. An overview of equations of state was given as a classification by type of model and its applicability to different problems. A comparative analysis on volumetric properties of pure components and binary systems was carried out.

The numerical developments of this work and their implementation in the reactive code HYTEC, were then used to model carbon capture and storage and to study the impact of convective mixing and gas impurities in the injected stream. The proposed coupling allowed to investigate complex gas injection scenarios, modeling a system including but not limited to four gases, and to reveal new behaviors of gas and liquid dynamics: high interaction between density driven flow and chromatographic partitioning and their effects on the velocity of the front of gas current and on the overall gas dissolution.

The Hele-Shaw problem of convective mixing in the aqueous phase in the presence of the gas and the capillary fringe was modeled. It was shown that in addition to the convective motion in the aqueous phase, the instabilities growing on the rough contact surface between two phases implied corresponding perturbations in the gas phase, leading to the gas mixing. Afterward, the convective mixing in both phases was also observed at the reservoir scale.

The 3D modeling of carbon storage provides relatively lower convective dissolution at the reservoir scale, since a high-resolution grid is critical to accurately characterize the mass flux. We examined the numerical convective underestimation on the radial geometry xiii applying an analytical formulation that was obtained by time piecewise approximation of non-convective analytical solution coupled with an analytical expression of convective transfer. The method yields a prediction of current radius at the top of the aquifer with convective dissolution. The study exhibited that the 2D radial modeling could be efficiently used to estimate the convective transfer, reducing the grid dimension and enhancing the convective mixing at the same time. for hydrocarbon fluids [START_REF] Twu | A new generalized alpha function for a cubic equation of state part 1. peng-robinson equation[END_REF]; (L), (M ), (N ) parameters, generalized for alkanes, aromatics, ketons, alcohols, refrigerants and gases at

Modeling multiphase flow and reactive transport

Coupling the multiphase flow and reactive transport has already been performed both by the reservoir and hydrogeochemical modeling communities. The compositional flow simulators, based on implicit methods, are usually extended to the chemical modeling by solving the fully coupled system, since the framework is initially devised and aimed at large nonlinear systems: GEM-GHG [START_REF] Nghiem | Modeling co2 storage in aquifers with a fully-coupled geochemical eos compositional simulator[END_REF], GPRS [START_REF] Fan | A fully-coupled flowreactive-transport formulation based on element conservation, with application to co 2 storage simulations[END_REF]. By contrast, MoReS [START_REF] Farajzadeh | Detailed modeling of the alkali/surfactant/polymer (asp) process by coupling a multipurpose reservoir simulator to the chemistry package phreeqc[END_REF]Wei, 2012) and UTCHEM [START_REF] Delshad | Utchem version 9.0 technical documentation[END_REF] use an operator-splitting approach. The element-based version of GPRS [START_REF] Fan | A fully-coupled flowreactive-transport formulation based on element conservation, with application to co 2 storage simulations[END_REF] and MoReS, coupled with PHREEQC [START_REF] Parkhurst | User's guide to PHREEQC (Version 2): A computer program for speciation, batch-reaction, onedimensional transport, and inverse geochemical calculations[END_REF], can solve the kinetic and equilibrium reactions in mineral phase. Only MoReS provides complex activity models and thermodynamic databases. Coupling existing codes is another alternative, it thus enforces the operator splitting routine; e.g., iCP [START_REF] Nardi | Interface comsol-phreeqc (icp), an e cient numerical framework for the solution of coupled multiphysics and geochemistry[END_REF], the interface of PHREEQC and COMSOL is an ongoing project1 .

Several reactive transport simulators have succeeded in integrating multiphase flow by an operator-splitting approach: PFLOTRAN [START_REF] Lu | High resolution numerical investigation on the e↵ect of convective instability on long term co2 storage in saline aquifers[END_REF], STOMP [START_REF] White | STOMP Subsurface Transport Over Multiple Phases, Version 4.0, User's Guide[END_REF], TOUGHREACT [START_REF] Xu | Coupled modeling of non-isothermal multiphase flow, solute transport and reactive chemistry in porous and fractured media: 1. model development and validation[END_REF]. In these codes, the twophase flow system is set for two components that are presented by water and air/CO 2 /CH 4 ;

the relevant thermodynamic/equation of state module should be chosen (A more detailed description of the coupling methods will be given in Sec. 3.1.3). In order to avoid large flow systems, the formulation can be reduced by considering only dominant components, TOUGHREACT. The applications of this method are limited.

A general method should be initiated to allow an arbitrary number of components, all the while preserving minimum dimension of the flow problem. This work proposes to employ a phase conservation formulation for solving compressible multiphase multicomponent flow by an operator-splitting approach. The phase formulation is coupled with the linear component conservation system, so that the flow module can be separated. Therefore, the dimension of multiphase system is independent of the number of components. The module can be integrated into the reactive transport framework as a conventional two-phase two-component or compositional flow module.

A new approach of multiphase modeling was implemented in the (SIA) reactive transport simulator, HYTEC (Lagneau and van der Lee, 2010a;[START_REF] Van Der Lee | Moduleoriented modeling of reactive transport with hytec[END_REF], et al., 2003;[START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the momas benchmark case[END_REF]Lagneau and van der Lee, 2010b;[START_REF] Trotignon | Design of a 2-d cementation experiment in porous medium using numerical simulation[END_REF] and it was applied in numerous applications, such as cement degradation (De Windt L., 2010), radioactive waste disposal [START_REF] Debure | Hlw glass dissolution in the presence of magnesium carbonate: Di↵usion cell experiment and coupled modeling of di↵usion and geochemical interactions[END_REF], geological storage of acid gases [START_REF] Corvisier | Impact of co-injected gases on co2 storage sites: Geochemical modeling of experimental results[END_REF], uranium in situ recovery processes [START_REF] Regnault | 3d reactive transport simulations of uranium in situ leaching : Forecast and process optimization[END_REF]. When modeling miscible flow, the choice of equation of state models can be decisive, especially at high pressure and temperature. In Chapter 4, the equation of state models are classified. For some of them, a comparative analysis on volumetric properties is carried out. The models of fluid properties are also discussed.

Outline

Chapter 5 concerns the numerical and analytical modeling of carbon capture and storage and the impact of gas impurity in the injected stream. The study starts from a 1D problem of chromatographic partitioning. Next, it is followed by a convective problem with a thin gas layer. Thereafter, we provide the numerical and analytical results of 2D radial CO 2 modeling. The problem is then extended to three scenarios of CO 2 injection containing gas impurities, modeling a system with up to four gases.

Chapter 6 draws some conclusions and recommendations for future work.

Chapter 2

Multiphase Flow

Le After introducing basic definitions related to hydrodynamics (Sec. 2.1), the general formulations of multiphase flow (Sec. 2.2) are presented in order to provide a foundation for the modeled processes and their further extension. One of these systems is then chosen for modeling under the hypotheses listed in Sec.

Definitions

Averaging over a representative elementary volume (REV), one can represent the fluid and rock properties passing from one region to another using a continuum approach. In this work, the macroscopic level is considered; i.e., approximately from 10 -3 to 10 m. We first introduce notions of phase, porosity, saturation and component, then the Darcy-Muskat velocity is given with the description of fluid-matrix parameters of multiphase systems.

In case two or more fluids are immiscible (i.e., they are separated by physically distinct interfaces), each of the fluids is a phase and the system is multiphase. The porosity φ is the ratio of porous space within an REV relative to the total volume of the REV. The saturation S α of phase α is the ratio of the volume occupied by the phase α in an REV relative to the total volume of pore space in the REV. Supposing that the fluids fill all porous space in the REV, the following relationship holds for the phase saturations:

α S α = 1 (2.1.1)
When fluid is displaced, a residual saturation S αr is typically retained; i.e., the residual saturation of the wetting phase during drainage, the residual saturation of the non-wetting phase during imbibition. Thus, the effective saturation is defined as

S α = S α -S αr 1 -β S βr , (2.1.2)
where α, β are the phases and then α S α = 1.

The composition of phase α can be described by mass/mole fractions, concentrations.

The mass fraction X α,k of the components k is the ratio of the mass of the component k to the total mass of the phase α, then k X α,k = 1.

(2.1.3) By extension of Darcy's law, the velocity for multiphase flow Darcy-Muskat [START_REF] Muskat | Flow of gas-liquid mixtures through sands[END_REF] is expressed as follows 2.1.4) where K is the intrinsic permeability tensor, k rα is the relative permeability of phase α, p α is the pressure of phase α, ρ α is the density of phase α, µ α is the viscosity of phase α and g is the gravity vector.

u α = - k rα µ α K (∇p α -ρ α g) , (
Due to the interfacial tension, the pressure changes from one phase to another and the difference is expressed by capillary pressure p c : (2.1.5) where p n is the pressure of the non-wetting phase and p w is the pressure of the wetting phase.

p n = p w + p c ,
Since the porous space is filled by several fluids, the relative permeability of each phase k rα represents a fraction of intrinsic permeability for phase α. The phase mobilities λ α = krα µα are also used, the total mobility λ tot is defined by

λ tot = α λ α . (2.1.6)
The most common empirical functions of relative permeability k rα and capillary pressure p c are Brooks-Corey (1954) and van Genuchten (1980) models.

The Brooks-Corey capillary pressure function [START_REF] Brooks | Hydraulic properties of porous media[END_REF]) is given by

p c (S w ) = p b S -1 λ w , (2.1.7)
where p b is the entry pressure, indicating a minimum non-wetting pressure required to enter into porous media saturated by the wetting phase; λ reflects pore size distribution,

[0.2, 3].
The van Genuchten capillary pressure model [START_REF] Van Genuchten | A closed form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] yields (2.1.8) where m = 1 -1 n , n > 1. The Brooks-Corey relative permeability model [START_REF] Brooks | Hydraulic properties of porous media[END_REF] uses the function of capillary pressure (2.1.7) with the same λ:

p c (S w ) = 1 α (S -1 m w -1) 1-m ,
k rw (S w ) = S 2+3λ λ w , (2.1.9) k rn (S n ) = S 2 n (1 -(1 -S n ) 2+λ λ ). (2.1.10)
The van Genuchten relative permeability [START_REF] Van Genuchten | A closed form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] is expressed by

k rw (S w ) = S ε w (1 -(1 -S n n-1 w ) n-1 n ) 2 , (2.1.11) k rn (S n ) = S γ n (1 -(1 -S n ) n n-1 ) 2(n-1) n , (2.1.12)
where ε = 1/2, γ = 1/3 [START_REF] Helmig | Multiphase flow and transport processes in the subsurface -a contribution to the modeling of hydrosystems[END_REF], [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated media[END_REF].

m = 1 -1 n , n > 1,

Overview of multiphase flow formulations

We introduce briefly a survey of formulations for multiphase flow, that takes place according to the modeled problems and fluid properties. If the fluids are supposed to be immisible, the general phase problem is solved, otherwise the compositional form should be applied to permit phase transition such as vaporization and dissolution. Detailed derivations of multiphase flow formulations are presented in [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF]; [START_REF] Barenblatt | Theory of unsteady filtration of fluids and gases[END_REF]; Bastian (1999); [START_REF] Bear | Dynamics of fluids in porous media[END_REF]; [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF]; [START_REF] Hassanizadeh | General conservation equations for multi-phase systems: 2. mass, momenta, energy, and entropy equations[END_REF] and Peaceman (1977).

Immiscible multiphase flow

By considering the immiscible multiphase flow (no mass phase transfer), the mass conservation equation is written by analogy to single-phase flow but for each phase α, stating

∂(φS α ρ α ) ∂t + ∇ • (ρ α u α ) = q α , (2.2.1)
where q α is the source term. Including the sum of saturations (2.1.1) and N f -1 capillary pressure equations (2.1.5), it constitutes a closed system of 2N f equations for 2N f unknowns. There exists several representations of system (2.2.1), so let us propose a concise list starting from the strongly coupled formulation to the weaker ones.

• Coupled general phase conservation equations (PS)

Substituting the velocity (2.1.4) in (2.2.1) results in the general phase formulation

∂(φS α ρ α ) ∂t -∇ • (ρ α λ α K (∇p α -ρ α g)) = q α . (2.2.2)
For each phase α, a so-called natural variable, S α or p α , is chosen, leading to N f unknowns. It is clear that these N f nonlinear equations (2.2.2) are strongly coupled.

Moreover, the nonlinear functions k rα and p c depend on the saturation. Supposing that p w and S n are the primary variables for the two-phase problem, the non-wetting phase conservation equation contains a term λ n ∇p c . Choosing either the Brooks-Corey models (BC) for relative permeabilities (2.1.9), (2.1.10) and capillary pressure (2.1.7) or the van Genuchten models (V G) (2.1.11), (2.1.12) and (2.1.8) leads to the following singularities of λ n ∇p c (µ α is neglected):

BC V G lim Sn→0 k rn p c (S n ) = 0 0 lim Sn→1 k rn p c (S n ) = ∞ ∞ (2.2.3)
Assuming (p n , S w ) as the primary variables yields the appearance of λ w ∇p c in the wetting phase conservation equation and implies

BC V G lim Sw→0 k rw p c (S w ) = 0 0 lim Sw→1 k rw p c (S w ) = < ∞ ∞ (2.2.4)
Hence, any choice of one of these pairs entails restrictions at the extreme values of saturation in the presence of capillary forces. To avoid this problem, the pressurepressure formulation can be used by extending the capillary pressure. In this case, the phase disappearance problem arises that will be discussed later, p.

∂φ ∂t + α 1 ρ α φS α ∂ρ α ∂t + ∇ρ α • u α + ∇ • u = α q α . (2.2.5)
Once the unknown for the pressure equation is chosen, the total velocity u = α u α and the saturation equation(s) can be reformulated. The phase pressure and global pressure formulations are discussed below.

-Phase pressure and saturation (DPPS)

Supposing the wetting phase pressure as primary variable allows to rewrite the velocities as

u = -λ tot K(∇p w + f n ∇p c - α f α ρ α g), (2.2.6
) 2.2.8) where f α = λ α / β λ β is the fractional flow. Substituting the non-wetting phase velocity, for example, yields the saturation equation under the form of

u w = f w u + f w λ n K (∇p c + (ρ w -ρ n )g) , (2.2.7) u n = f n u -f n λ w K (∇p c -(ρ n -ρ w )g) , (
∂(φS n ρ n ) ∂t + ∇ • (ρ n u n ) = q n .
(2.2.9)

-Global pressure and saturation (GPS)

By defining the global pressure p such that ∇p = α f α ∇p α , the total velocity simply becomes

u = -λ tot K(∇p - α f α ρ α g).
(2.2.10)

Using the total velocity, the wetting phase saturation equation can be expressed by following .2.11) The main difference between the two formulations is the capillary pressure term λ n ∇p c that appears in the pressure equation of the DPPS formulation through the total velocity (2.2.6) as well as in the PS formulation (2.2.2). In the incompressible case, this nonlinearity is eliminated (2.2.10) by definition of the global pressure.

∂(φS w ρ w ) ∂t + ∇ • (ρ w (f w u + f w λ n K (∇p c + (ρ w -ρ n )g))) = q w . ( 2
When the compressibility is taken into account, the pressure equation (2.2.5) involves the phase velocities and then λ w λ n ∇p c , where λ w λ n vanishes at S α = 0, ∀α.

However, the global pressure p is not a physical variable, consequently the initial and Dirichlet boundary conditions cannot be directly imposed in the presence of capillary forces.

About numerical methods

Hence, we are left with the (near-)elliptic pressure equation, the parabolic/hyperbolic saturation equation(s). It allows to employ different numerical methods dedicated for these types of equations and solve them separately to reduce the degree of freedom, especially when the flow is multicomponent and miscible. Finally, it is worth pointing out that the pressure equation (2.2.5) must be rewritten in a conservative way in order to apply the Finite Volume Method (FVM).

Being proposed in [START_REF] Chavent | A new formulation of diphasic incompressible flows in porous media[END_REF], the idea of decoupling the pressure equation is extended to the multicomponent case and used in multiscale methods since it allows to evaluate the velocity field. For decoupled formulations 2.2.1, the sequential solution is CHAPTER 2. MULTIPHASE FLOW usually applied -IMPES (IMplicit Pressure Explicit Saturation) method, where the pressure equation is first solved implicitly, the explicit discretization of the saturation equation is then employed with fixed fluxes (similar to sequential modeling single-phase flow and multicomponent transport). The FIM (Fully Implicit Method) is usually employed for the PS formulation that requires to solve the whole system of nonlinear equations simultaneously. The AIM (Adaptive Implicit Method) allows to combine these methods by applying different levels of implicitness to each grid block dynamically in time [START_REF] Thomas | Reservoir simulation using an adaptive implicit method[END_REF][START_REF] Bibliography | Stability analysis and switching criteria for adaptive implicit methods based on the cfl condition[END_REF]. The comparative study of three methods was presented in [START_REF] Marcondes | A comparative study of implicit and explicit methods using unstructured voronoi meshes in petroleum reservoir simulation[END_REF] and highlighted the advantages of the AIM.

Compositional flow

When the phase transition is possible, the mass balance is formulated for each component k, k = 1, . . . , N s . So, the compositional flow system consists of N s partial differential equations

∂ ∂t α φS α ρ α X α,k + ∇ • α ρ α X α,k u α = q k (2.2.12)
for 2N f + N f N s natural variables: 2N f phase pressures and saturations, N f N s mass fractions. To close the system, we add N f -1 capillary pressure relations (2.1.5), the total saturation (2.1.1), N f mass fraction (2.1.3) sums and (N f -1)N s relations of phase equilibria. For instance, assuming N f = 2 and α, β = {g, l}, the natural variables (Coats, 1980) can be taken as the primary variables: p g , S l and N s -2 mass fractions X g,k . While the phase α disappears, the set of variables changes to p g and N s -1 mass fractions X β,k , β = α. The set of phases and their composition are usually defined by the phase stability test (Michelsen, 1982a) and flash calculation [START_REF] Michelsen | The isothermal flash problem. part ii. phase-split calculation[END_REF]. The formulations and methods for solving compositional flow will be discussed later in Sec. 3.1.1.

Intermediate case

The black-oil model, aimed at modeling the hydrocarbon system, represents an alternative between immiscible and miscible flow. There are three phases (water, oil, gas) and three components (water, oil, gas). The oil-gas interactions are allowed, while the water phase is unreactive with the rest of two phases. As an extension of the black-oil model, the 2.3. IMMISCIBLE TWO-PHASE FLOW MODEL limited compositional model was developed to eliminates the last restriction. To introduce solubility, vaporization, the phase transfer is expressed by a volume ratio (normalized volume change due to reaction), which consequently appears in the accumulation term.

The complete description of three-phase flow for hydrocarbon systems can be found in Peaceman (1977).

Immiscible two-phase flow model

The present development of multiphase flow module is seen as a preliminary work before coupling with the reactive transport framework HYTEC. The following hypotheses are made:

• the immiscible two-phase flow formulation, α = {w, n};

• some amount of water is always present in REV (CHESS numerical limitation: the chemical speciation library needs non-zero quantities of water to avoid infinite concentrations);

• isotropy of the intrinsic permeability, K = KI;

• isothermal system, no heat conservation equation;

• no pore compressibility ∂φ/∂p = 0, no geochemical deformation;

• incompressible flow.

Since the wetting phase does not disappear, the PS formulation can be applied, where p = p w and S = S n are the primary variables. The flow problem consists of two PDEs on

Ω×(0,T ), Ω ∈ R N , N = 1, 2, 3 ∂ ∂t φ(1 -S)ρ w -∇ • (ρ w λ w K (∇p -ρ w g)) = q w , (2.3.1) ∂ ∂t φSρ n -∇ • (ρ n λ n K (∇(p + p c ) -ρ n g)) = q n (2.3.2)
with the boundary conditions

p(x, t) = p d (x, t) on Γ D w , ρ w u w • n = ψ w (x, t) on Γ N w , S(x, t) = S d (x, t) on Γ D n , ρ n u n • n = ψ n (x, t) on Γ N n (2.3.3) CHAPTER 2. MULTIPHASE FLOW
and with the following initial conditions p(x, 0) = p 0 (x), S(x, 0) = S 0 (x), x ∈ Ω. (2.3.4) Under the taken assumptions, porosity, viscosity and density are constant. The equations (2.3.1) and(2.3.2) can be divided by ρ w and ρ n , respectively. To close the system, the relative permeability and capillary pressure models should be imposed.

Further developments on the coupling (Ch. 3) will extend this formulation to compressible, viscous and miscible flow. Modeling non-isothermal flow and anisotropy is outside the scope of this work.

Numerical method

We present the description of numerical methodology chosen for solving the two-phase flow problem, for the PS formulation in particular: first, the discretization in space and time;

second, the linearization and Newton's algorithm; next, solving linear system and finally, the Jacobian structure.

Discretization

Volume averaging

The FVM based on Voronoï meshes is employed for Eqs. (2.3.1) and (2.3.2). Let us introduce the following definitions:

Ω i is the element with the center x i such that Ω = i∈{1,N } Ω i , |Ω i | = Ω i dω. (2.4.1)
Γ ij = ∂Ω i ∩ ∂Ω j is the interface between the elements Ω i and Ω j ,

|Γ ij | = Γ ij dγ.
(2.4.2)

I int i = {j : Γ ij = ∅}
is the index set of element Ω j that has the common edge with the element Ω i . .4.3) where the flux (2.4.4) (2.4.5) where F α is the approximated normal flux F α • n. The Neumann conditions are taken into account by integrating the imposed fluxes, while the Dirichlet approximation is given by (2.4.6)

I N α,il = {l : Γ il = ∂Ω i ∩ Γ N α,l = ∅}
Ω i ∇ • F α dω = ∂Ω i F α • n dγ, ( 2
F α = ρ α u α ,
|Ω i | ∂ ∂t φ i (1 -S i )ρ w,i + j∈I int ij |Γ ij |F w,ij + j∈I N w,il |Γ il |(ψ w ) ij = |Ω i |(q w ) i ,
|Ω i | ∂ ∂t φ i S i ρ n,i + j∈I int ij |Γ ij |F n,ij + j∈I N n,il |Γ il |(ψ n ) ij = |Ω i |(q n ) i ,
p i = p d,i ∀i ∈ I D w ,
S i = S d,i ∀i ∈ I D n .
(2.4.7)

Next, the two point flux approximation (TPFA) is employed for the interior fluxes: (2.4.8) (2.4.9) where (•) ij is the weighting scheme for interface parameters between Ω i and Ω j , g = g • n.

F w = -ρ w k rw µ w ij K ij p j -p i x j -x i -ρ w,i g ,
F n = -ρ n k rn µ n ij K ij p j -p i x j -x i + p c,j -p c,i x j -x i -ρ n,i g ,
The intrinsic permeability at the interface K ij can be evaluated as harmonic or upstream weighting, while for the relative permeability, the upstream treatment is prerequisite

k rα,ij =      k rα,i if p α,j -p α,i x j -x i -ρ α,i g ≥ 0,
k rα,j else.

(2.4.10)

Implicit (Backward Euler) in time

The Fully Implicit (Backward Euler) method is used because of its unconditional stability.

Assuming t n ∈ (0, T ) such that

0 = t 0 < • • • < t n < t n+1 < • • • < t M = T, ∆t n+1 = t n+1 -t n , n=0,...,M ∆t n+1 = T,
and the mass accumulation A α,i = φ i S α,i ρ α,i , the discretization yields .4.11) The relative permeability is discretized by implicit upstream treatment (update at every Newton iteration) guaranteeing stability but increasing truncation errors at the same time, as discussed below.

|Ω i | A n+1 α,i -A n α,i ∆t n+1 + j∈I int ij |Γ ij |F n+1 α + j∈I N α,il |Γ il |(ψ α ) n+1 ij = |Ω i |(q α ) n+1 i . ( 2
Weighting of nonlinearities, stability [START_REF] Settari | Treatment of nonlinear terms in the numerical solution of partial di↵erential equations for multiphase flow in porous media[END_REF], [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF], Peaceman (1977) demonstrated the impact of relative permeability approximation at the interface, in space (from midpoint rule to upstream) and in time (from explicit to implicit) on the accuracy and stability. [START_REF] Blair | Solution of two-phase flow problems using implicit di↵erence equations[END_REF] presented the numerical results of the implicit method with the implicit and explicit mobilities and highlighted the efficiency of the fully implicit scheme. Despite its unconditional stability, it is not reasonable to use any large time step. Since such stability is associated with numerical dispersion, the effect of smearing the sharp solution appears (in hyperbolic flows), whose degree can be reduced by grid refining, high-order discretization and total-variation-diminishing (TVD) schemes. Allen (1984) showed the numerical capillary diffusion at the Buckley-Leverett problem with upwind flux discretization based on FEM (Finite Element Method). Later, Bastian (1999) compared the implicit and the Crank-Nicolson scheme with upstream weighting of mobilities and also the implicit method with upstream and central differences treatment. We conclude that the fully implicit method with the implicit upstream mobilities can be carefully employed, as in the aimed applications, the capillary forces are indeed present (parabolic case). Moreover, our ultimate goal is to extend the two-phase flow module to multicomponent flow, Ch. 3.

Solving linearized system

Solving the nonlinear system: Newton's method

The discretization Eq. (2.4.11) yields the system of nonlinear algebraic equations (2.4.12) where F is the residual function with dim(F ) = 2N , ζ is the vector comprising 2N primary variables: .4.13) The alternative ordering ζ = (p, S) T is also possible. The approximative solution of the differential function F can be found by Newton's method. Arising from the 1 st Taylor polynomial, the estimation of F at (k + 1) th approximation of root ζ k+1 evolves into 2.4.14) where J is the Jacobian of the vector function F : .4.15) Since ζ k+1 is the approximative root, the reformulation of Eq. (2.4.14) leads to solving the linearized system 2.4.16) where δζ k+1 = ζ k+1ζ k . To summarize, let us present the iterative procedure in Algorithm 1. Once the error tolerance ε N f , usually ε N f ∈ [10 -8 , 10 -6 ], and the maximum number of iterations k max are imposed, we deal with the linearized problem to find δζ k+1 , line 5. The iterative method continues until the norm of the residual function is small enough, line 4. If k > k f l max , the time step ∆t n+1 in Eq. (2.4.11) is reduced and Algorithm 1 is repeated.

F (ζ) = 0,
ζ = (p 1 , S 1 , . . . , p N , S N ) T . ( 2
F (ζ k+1 ) = F (ζ k ) + J (ζ k )(ζ k+1 -ζ k ), (
J ij = ∂F i ∂ζ j . ( 2
J (ζ k )δζ k+1 = -F (ζ k ), (
Algorithm 1 Newton's method for Immiscible Two-Phase Flow

1: ε N f = 1 × 10 -6 2: k f l max = 9 3: k = 0 4: while F (ζ k ) ∞ ≥ ε N f F (ζ 0 ) ∞ ∧ k ≤ k f l max do 5: find δζ k+1 : J (ζ k )δζ k+1 = -F (ζ k ) 6: ζ k+1 ← ζ k + δζ k+1 7:
k ← k + 1

8: end while Solving the linear system: GMRES method

The linear system (2.4.16) is solved at each Newton iteration. The matrix J has a nonsymmetric sparse structure and becomes larger with increasing N as expected. Hence, the choice of linear solver is of great importance in multiphase modeling inasmuch as it affects the feasibility and the performance of the whole module. The overview of iterative methods for solving the linear systems can be found in [START_REF] Barrett | Templates for the solution of linear systems: building blocks for iterative methods[END_REF][START_REF] Benzi | Preconditioning techniques for large linear systems: a survey[END_REF]. It is well known that applying a preconditioner improves the condition number of the matrix and therefore, the convergence rate of compatible linear solver. The preconditioning is regularly employed with the Krylov methods; e.g., GMRES [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]Schultz, 1986), BiCGSTAB (Van der Vorst, 1992), that are principal solvers for non-symmetric system. [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] gave an excellent description of the Krylov subspace methods and preconditioning.

NUMERICAL METHOD 21

The GMRES method with ILU(0) provided by Sparselib++ [START_REF] Pozo | Sparselib++ v. 1. 5 sparse matrix class library. reference guide[END_REF] is employed for the system (2.4.16). The tolerance of GMRES ε lin impacts to the convergence rate of Newton's algorithm and it was suggested that ε lin ∈ [10 -8 , 10 -5 ]. The restart option can be defined in range [10 2 , 10 4 ] that depends on the condition number of the Jacobian. The split preconditioning is used for ill-conditioned matrix: left diagonal (explicit) and then right ILU, (Saad, 2003, Ch. 9).

Remarks: The two-stage constrained pressure residual (CPR) preconditioner, initially proposed by [START_REF] Wallis | Constrained residual acceleration of conjugate residual methods[END_REF], significantly accelerates the calculations of linear systems arising from the fully implicit scheme [START_REF] Cao | Parallel scalable unstructured cpr-type linear solver for reservoir simulation[END_REF][START_REF] Voskov | Comparison of nonlinear formulations for two-phase multi-component eos based simulation[END_REF].

The algorithm is designed to take into account the types of differential equations and to employ the relevant methods. For instance, the two-phase flow system comprises a near elliptic pressure equation and a near hyperbolic saturation equation (composition equations). After the linearisation, the system is decomposed into subsystems. First, the pressure subsystem is solved by a linear solver, an Algebraic Multigrid solver [START_REF] Stüben | Algebraic multigrid (amg): experiences and comparisons[END_REF]; the residual of the entire system is then corrected to the obtained pressure solution.

Second, a suitable preconditioner, such as ILU(0), is applied to the corrected system.

Analytical Jacobian

The two-phase PS formulation entails the numerical or analytical 2N × 2N Jacobian calculation. The matrix structure is 2 × 2 block whose dimension corresponds to the number of equations and to the number of primary variables. The Dirichlet conditions (2.4.6) yield trivial blocks.

The analytical estimation is applied in HYTEC. Now, we detail the Jacobian structure.

Neglecting the time indices, we derive the diagonal block of incompressible PS problem with the primary variables (p w , S n ) for Ω i (2.4.17) and therefore the contribution of i th block to the non-diagonal block j (2.4.18). Since the Jacobian should be invertible matrix, J 2i,2i = 0 in the diagonal block (2.4.17). Thus, there is a constraint at the extreme value of saturation S n = 1, when (p w , S n ) are primary variables. The restriction S w = 1 appears in case of (p n , S w ) choice by analogy.

The models of relative permeability and capillary pressure are user-defined. The models of fluid properties will be discussed in Ch. 4. Once all property models are selected, CHAPTER 2. MULTIPHASE FLOW the Jacobian is automatically assembled. The matrix is constructed by gathering the corresponding functions and their analytical derivatives and it is written in the CRS (compressed row storage) format [START_REF] Pozo | Sparselib++ v. 1. 5 sparse matrix class library. reference guide[END_REF]. The nonlinear functions must then be of class C 1 ([0, 1)). Some extensions and redefinitions of fluid-matrix functions and consequently their derivatives should be made at least on [0, 1ε], that is detailed for capillary pressure and relative permeability functions in Appendices A and B, respectively.

J 2i,2i = ∆t j∈I int ij |Γ ij | ρ w k rw µ w ij K ij 1 x j -x i J 2i,2i+1 = -|Ω i |φ i ρ w,i -∆t j∈I int ij |Γ ij | ρ w µ w ij ∂k rw,ij ∂S i K ij p j -p i x j -x i -ρ w,i g J 2i+1,2i = ∆t j∈I int ij |Γ ij | ρ n k rn µ n ij K ij 1 x j -x i J 2i+1,2i+1 = |Ω i |φ i ρ n,i -∆t j∈I int ij |Γ ij | ρ n µ n ij ∂k rn,ij ∂S i K ij p j -p i x j -x i + p c,j -p c,i x j -x i -ρ n,i g +∆t j∈I int ij |Γ ij | ρ n k rn µ n ij K ij ∂p c,i ∂S i 1 x j -x i (2.4.17) J 2i,2j = -∆t|Γ ij | ρ w k rw µ w ij K ij 1 x j -x i J 2i,2j+1 = -∆t|Γ ij | ρ w µ w ij ∂k rw,ij ∂S j K ij p j -p i x j -x i -ρ w,i g J 2i+1,2j = -∆t|Γ ij | ρ n k rn µ n ij K ij 1 x j -x i J 2i+1,2j+1 = -∆t|Γ ij | ρ n µ n ij ∂k rn,ij ∂S j K ij p j -p i x j -x i + p c,j -p c,i x j -x i -ρ g,i g -∆t|Γ ij | ρ n k rn µ n ij K ij ∂p c,j ∂S j 1 x j -x i .
(2.4.18)

Verification

We present the numerical results of a 1D radial parabolic problem with self-similar solution and a 2D hyperbolic problem.

Self-similar solution to two-phase displacement problem

Modeling the problem obeying the self-similarity allows not only to analyze the physical phenomena, intermediate asymptotic behavior [START_REF] Barenblatt | Scaling, self-similarity, and intermediate asymptotics: dimensional analysis and intermediate asymptotics[END_REF] but also to verify the numerical code. The invariant solution ξ = x/a √ t was initially discovered for the heat equation and then for several (parabolic PDE) infiltration problems and immiscible displacement [START_REF] Polubarinova-Kochina | Theory of ground water movement[END_REF]Barenblatt et al., 1972, Ch. 4 and 6). The latter carried out the self-similarity study considering the capillary pressure effects, radial flow, double porosity; afterward, [START_REF] Chen | Some invariant solutions to two-phase fluid displacement problems including capillary effect[END_REF] , 1981).

Model description

During the workshop Intercomparison of numerical simulation codes for geologic disposal of CO 2 initiated by Lawrence Berkeley National Laboratory (Pruess et al., 2002), the axisymmetric problem of constant injection in a long, horizontal, reservoir was stated.

Γ N : ψ l = 0 ψ g = 1/(2π0.3) [kg/m 2 /s] Γ D : p l = 120 [bar] S g = 0 ∂Ω\(Γ D ∪ Γ N ): ψ l = 0 ψ g =

Results

The numerical solution of gas saturation S g is plotted as a function of the self-similar variable R 2 /t in Fig. 2.4(a). The results are presented for the fixed time (30, 10 2 , 10 3 , 10 4 d)

and for the fixed radial distance (25.8, 1100.3 m). The numerical solution approximately admits the invariant variable. The time step varies from 5 s to the imposed maximum, 50 d.

Even taking into account the modifications of the problem (no evaporation of the residual water), the results obtained by HYTEC are in a good agreement with those of TOUGH2-ECO2, Fig. 2.4(b). We detail the differences below. The saturation profile modeled by HYTEC has two regions: to the neglected effects of solubility, fluid and rock compressibility, viscidity, whose role will be shown in Sec. 3.2. Nevertheless, the similarity is preserved and the solution is within the range of the results given by other codes (Fig. 2.6), [START_REF] Pruess | Code intercomparison builds confidence in numerical simulation models for geologic disposal of CO 2[END_REF]. Moreover, the HYTEC's profile is close to that of GEM simulator (ARC in the Fig. 2.6), which also lacks drying-out capabilities. By neglecting capillarity, the two-phase system results in the hyperbolic equation which transforms to a Riemann problem under discontinuous initial conditions. The solution is therefore analyzed by deriving characteristic families in order to determine the wave type -rarefaction, shock. The latter forms the discontinuity of solution, that can introduce numerical instability and oscillations. As was mentioned in Sec. 2.4.1, the chosen numerical method smothers the steep solutions. However, it must converge to the physical meaningful solution. The Riemann problem was extensively studied in gas thermodynamics. The detailed physical explanation of phenomena can be found in Loitsyansky (1950, Ch. 4) and Landau and Lifshitz (1986, Ch. 9-11); the mathematical and numerical study was given in Samarskii and Popov (1992, Ch. 1 and 5) and [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. To solve the numerical difficulties of Riemann problem, the Godunov method was proposed that is actually the integrated finite difference or finite volume method. Therefore, applying the FVM allows to model the hyperbolic two-phase flow problems.

• R 2 /t < 2 × 10 -2 [m 2 /s] that indicates two-phase regime and • R 2 /t ≥ 2 × 10 -2 [m 2 /
In the context of two-phase flow, we refer to the survey by Bastian (1999, Sec. 2.4.2), since the velocity (that becomes a fractional flow function in two-phase flow) depends on the mobility ratio (M = λ w /λ n for waterflooding case), and therefore on the relative permeability models and viscosity ratio. The characteristics analysis yields that the rarefaction and shock wave appearance is subject to the mobility ratio and the inflection point of fractional flow.

Model description

The model of waterflooding process in a 2D homogeneous horizontal reservoir was proposed by Bastian (1999, Sec. 7.2.1). The capillary pressure and gravity are then neglected. By transforming the notion of M , the viscosity ratio µ is defined as the ratio of the displaced phase viscosity to the injected one: µ = µ n /µ w = 20 when the water is injected in more viscous fluid. Applying the characteristic method reveals that the solution possesses CHAPTER 2. MULTIPHASE FLOW different wave states, the rarefaction zone turns to the shock through the tangential point of the fractional flow, see Bastian (1999, Sec. 2.4.2). Hence, reducing the viscosity ratio, for ex. to 1/20, provides the gradient decrease of rarefaction zone, the solution is steeper, more shock-like than that in case of µ = 20. So, the problem of µ = 1/20, is stated with the same parameters as those of case with µ = 20.

Results

The contour lines of wetting phase saturation S w for µ = 20 are given on the left column of Fig. 2.9. Refining the grid displays the convergence of the method. As discussed above, the rarefaction zone is followed by the shock front. The PS formulation with primary variable (p n , S w ) was used in Bastian (1999) and the initial condition is S n0 = 1, while the chosen formulation PS with pair (p w , S n ) enforces to adapt it, S n0 = 0.999. In spite of this difference, the results of UG (Bastian, 1999) and HYTEC are similar, Fig. 2.10.

The case of µ = 1/20 is demonstrated on the right column of Fig. 2.9. Again, the convergence of the method can be illustrated using grid refining. The solution exhibits the gradient decrease -the front sharpening. The injected non-wetting fluid of higher viscosity advances by occupying more volume and the front progress is then slower, comparing with the injection of less viscous wetting fluid.

The method allows to model the rarefaction and shock waves approximately. The convergence and the importance of grid scale were demonstrated The results are close to those in the reference and provide the physical solution.

Discontinuous permeability

Model description

The model of discontinuous permeability, posed in Bastian (1999, Sec. 7.2.3), points out the numerical difficulties associated with the discontinuity of intrinsic permeability. Based on the hyperbolic problem with homogeneous field Table 2.2, the next discontinuity subzones are added, Fig. 2.7(b):

Ω 1 = [37.5, 135] × [90, 232.5] [m], Ω 2 = [30, 142.5] × [172.5, 240] [m],
where (2.5.1) and

K =      10 -16 if x ∈ Ω 1 , 10 -10 otherwise,
S w0 =      0.8 (0.2) if x ∈ Ω 2 (and if µ = 1/20), S w0 otherwise.
(2.5.2) Then, the intersection of them Ω 3 = Ω 1 ∩Ω 2 presents a low permeability region, containing some amount of injected fluid -80% saturated. The analogous model is stated for µ = 1/20. Therefore, the initial heterogeneity is provided by the discontinuous saturation and absolute and relative permeability distribution. . 2.11 shows the results obtained by UG and HYTEC. The contour lines of wetting phase saturation are similar, although there are some differences. HYTEC's front of wetting fluid from Ω 2 \Ω 1 progresses faster than that of UG, that can be due to the different discretization methods and definition of Ω 1,2 (Ω 1,2 are slightly shifted comparing with those of Bastian (1999, Sec. 7.2.3)).

Results

Fig

For both viscosity ratios, the convergence of method can be seen on Fig. 2.12, where that is expressed by the front sharpening and alignment of Ω 1 , as expected. In waterflooding case µ = 20, the apparent impact of grid refining is additionally highlighted by the stretched form of front. We have seen the nature of solution with different viscosity ratios on the previous model, Sec. 2.5.2. The solution behavior in this model is analogous

CHAPTER 3. COUPLING MMF AND RT

Reactive transport simulators aim at modeling complex biogeochemical systems. It is of great importance to be able to extend its abilities to multiphase flow due to its relevance to applications covering a wide range of subsurface activities. This chapter is dedicated to the coupling of the multiphase multicomponent flow (MMF) and reactive transport (RT).

The main numerical methods are described for each problem in Sec. 3.1. Sec. 3.2 presents a new coupling approach and the first numerical results. The verification of the chosen methods including a benchmark proposition can be found in Sec. 3.3.

Overview of multiphase multicomponent flow and reactive transport

The section gives conceptual descriptions of MMF (Sec. 3.1.1) and RT (Sec. 3.1.2) modeling and their mathematical/numerical issues. A list of existing coupling methods is provided in Sec. 3.1.3, before presenting the actual coupling, Sec. 3.2.

Multiphase multicomponent flow

Numerical methods

The formulation of miscible/compositional flow, presented in Sec. 2.2.2, Eq. (2.2.12), involves the mass conservation equation for each component:

∂ ∂t α φS α ρ α X α,k + ∇ • ( α ρ α X α,k u α ) = q k (3.1.1)
and results in N s nonlinear equations. The problem can be solved by FIM or by implicit pressure explicit component (IMPEC) for decoupled formulation, similar to IMPES.

The fully implicit solution involves dividing the compositional problem into the primary and the secondary sets. The equations that need the information of the adjacent cells correspond to the primary set (e.g., the mass balance equations); otherwise, they are associated with the secondary set; e.g., the vapor-liquid equilibria (VLE) equations (that will be discussed below), the saturation, mole/mass fraction constraints. The last set is decoupled and used to eliminate the secondary variables.

By applying the IMPES/IMPEC scheme, the degree of freedom of the primary set is reduced at least twice in comparison with that of FIM, as only the pressure equations constitute the primary system. At the same time, explicitness imposes time step restriction, arising from the stability analysis. Analogous to the single phase flow, the Courant-Friedrichs-Lewy (CFL) condition for IMPES must be defined:

∆t ≤ CFL min V tot φ u in , V tot φ u out (3.1.2)
where CFL ≤ 1, min is the minimum function over all elements Ω i of the modeled domain, V tot is the volume of the cell, u in/out = ∂Ω i α u α,in/out • n dγ is the inflow/outflow rate. Peaceman (1977) and [START_REF] Settari | Treatment of nonlinear terms in the numerical solution of partial di↵erential equations for multiphase flow in porous media[END_REF] gave the stability analysis and the corresponding CFL condition with regard to the different approximations of relative permeability. [START_REF] Coats | Impes stability: selection of stable timesteps[END_REF] expanded the CFL definition to three phase multicomponent flow. In the advective dominant regime, the restriction can be an important factor if large time steps are required. The AIM combines the advantages of both of them: the method significantly reduces the number of implicit variables, linear solver and total time per Newton iteration and increases the accuracy of solution [START_REF] Zhou | Automatic differentiation framework for compositional simulation on unstructured grids with multi-point discretization schemes[END_REF].

Let us conclude, the IMPEC formulation is advantageous in terms of the number of nonlinear equations in the primary set at the cost of the time step restriction, while the FIM proposes the unconditional stability. Here, we have to cite [START_REF] Settari | Treatment of nonlinear terms in the numerical solution of partial di↵erential equations for multiphase flow in porous media[END_REF]:

As the implicitness of the method increases, stability improves, but truncation errors also increase.

Thus, constraints to the saturation, mole/mass fraction and pressure deviations should be imposed. When applying the AIM, the appropriate CFL switching criteria are calculated to define the implicit level of each cell [START_REF] Zhou | Automatic differentiation framework for compositional simulation on unstructured grids with multi-point discretization schemes[END_REF]. The detailed description of IMPEC, FIM can be found in [START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF].

Phase equilibria

We briefly recall the phase equilibria calculation, relying on the thermodynamics basics [START_REF] Firoozabadi | Thermodynamics of hydrocarbon reservoirs[END_REF][START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF]Smith and van Ness, 1987). The compositional flow system (3.

f v i (T, P, ȳ) = f l i (T, P, x) , (3.1.3)
where x and ȳ are the set of mole fractions in the liquid and gas phases. We first present the ϕϕ approach for solving the phase equilibria and give a general description of equation of state (EOS). Then the γϕ approach and the numerical methods are discussed.

By substituting the fugacity coefficients ϕ α i in both sides of the phase equilibirium equation (3.1.3) for species i in solution (ϕϕ approach), one can derive: 1.4) where x i and y i are the mole fractions in the liquid and vapor phases, respectively.

y i P ϕ v i = x i P ϕ l i , (3.
Eq. ( 3.1.4) can be rewritten: 1.5) where the fugacity coefficients can be calculated by EOS.

ln(ϕ v i y i ) = ln f v i P = ln f l i P = ln(ϕ l i x i ), (3.
An EOS establishes a relation between pressure, molar volume and temperature (P vT )

and allows to evaluate fugacity, solubility, enthalpy and density. For any pure homogeneous fluid in the equilibrium state and for the region of single phase, this relation can be expressed by f (P, v, T ) = 0. The simplest one is the ideal gas law:

P v = RT. (3.1.6)
The compressibility factor, defined as Z = P v/RT , is equal to unity if the gas is ideal.

Otherwise, if the gas is real, the models considering the intermolecular forces are used. In 1873, J.D. van der Waals proposed the cubic EOS, including the repulsive and attractive terms: The cubic EOS can be reformulated as a cubic equation on Z, that has three roots for a given P , two of them can be complex. Depending on the number of phases, the real roots 3.1. OVERVIEW OF MMF AND RT 41 correspond to the Z in the liquid or vapor phases. Note that the EOS is separately solved for liquid and vapor phases by using the set x and ȳ, respectively. The cubic equation can be solved for v either analytically or numerically (by an iterative method). Once the roots are found, the calculation of fugacity coefficients ln ϕ α i of a pure component or of a component in a mixture is straightforward by the EOS. The latter is performed by means of the binary interaction parameters and the mixing rules [START_REF] Peng | A new two-constant equation of state, industrial & engineering chemistry fundamentals[END_REF][START_REF] Firoozabadi | Thermodynamics of hydrocarbon reservoirs[END_REF]Smith and van Ness, 1987;[START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF]. Another advantage of EOS is the unique way of P vT properties and VLE calculation even in the critical region. The fully implicit solution of the EOS compositional flow was proposed in the pioneer work of Coats (1980).

P = P R + P A ,
The EOS equations cannot embrace the wide range of chemical systems, the activity models are then employed to describe the (real) behavior of species in solution. Therefore, the RHS of VLE Eq. (3.1.4) is expressed by applying either the symmetric or the asymmetric activity coefficient γ (γϕ approach). The latter uses the reference state in a solvent. Denoting the symmetric activity coefficient of species

i at infinite dilution in a solvent k by γ ∞ i such that γ ∞ i = lim x i →0 γ i , i = k yields f l i = K h,∞ i γi x i , (3.1.7)
where γi = γ i /γ ∞ i is the assymetric activity coefficient. Hence, it follows by the definition of γi that

lim x k →1 γk = 1, (3.1.8) lim x i →0, i =k γi = 1. (3.1.9) K h,∞ i
comprises the Henry's law constant K h i , representing the solute-solvent relations, and the Poynting factor, providing the pressure effect on K h i : 1.10) where v∞ i (T, P ) is the partial molar volume of species i at the infinite dilution in the solvent k (e.g., water). Both approaches (ϕϕ and γϕ) are connected: γ i = φ i /φ 0 at given T and P , where φ 0 is the fugacity coefficient at a reference state. The detailed derivation can be found in [START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF].

∂ ln K h i ∂P T = v∞ i (T, P ) RT , (3.

CHAPTER 3. COUPLING MMF AND RT

The phase equilibria can be performed by flash calculation. The K-value/equilibrium ratio is defined as K i = y i /x i . By applying either ϕϕ or ϕγ approach, one derives .1.11) The K-values can be calculated from the EOS, activity models, simplified and tabulated as P T dependent coefficients. Under the assumption of ideal mixture, the K-values are composition independent.

K i = ϕ l i ϕ v i = K h,∞ i γi P ϕ v i . ( 3
Supposing that V and L are the phase mole fractions, V + L = 1 and that (3.1.12) where z i is the overall mole fraction of species i, the fractions can be expressed as

z i = x i L + y i V,
x i = z i 1 + V (K i -1) , ( 3
y i = z i K i 1 + V (K i -1) . (3.1.14)
The constraint on mole fractions Ns i=1

(y i -x i ) = 0 (3.1.15)
evolves into the Rachford-Rice equation (Smith and van Ness, 1987)

F (V ) = Ns i z i (K i -1) 1 + V (K i -1) = 0. ( 3

.1.16)

There are 2N s + 1 unknowns, V (or L), x and y and the system of 2N s + 1 equations, Eqs. (3.1.3), (3.1.12) and (3.1.15). The successive substitution (SS) method and Newton's method are usually employed to solve the system. At each SS iteration, we solve the Rachford-Rice equation on V at given P , T and z by Newton's method; then the mole fractions are calculated by Eqs. (3.1.16), (3.1.13) and (3.1.14); the K-values are updated through Eq. (3.1.3). When the estimate is close to the solution, Newton's method is applied to the system, Eqs. (3.1.3), (3.1.12) and (3.1.15), that accelerates the convergence. The analysis of the Rachford-Rice function, its derivative, the convergence and the stability of flash calculation are given in [START_REF] Firoozabadi | Thermodynamics of hydrocarbon reservoirs[END_REF], [START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF]. Different approaches of the flash calculation exist but the subject is outside the scope of this study.

Solving phase equilibria can be either integrated in the Jacobian of compositional flow system or separately coupled with the conservation system. For each cell of single phase, the phase stability test (Michelsen, 1982a) is performed to capture the switch to different phase region, single or two-phase. Then, for each cell of two-phase, the phase composition is defined through the thermodynamic flash.

The important difference between using the EOS models and flash is the way of calculation: the majority of EOS models can be performed analytically, while the flash requires the iterative procedure if a real mixture is considered. Note that the EOS approach is equivalent to the K-values definition and one iteration of flash calculation.

Primary variables

Solving the compositional flow system (3.1.1) by FIM or IMPEC, several formulations of primary variables can be applied to deal with the phase (dis)appearance. The commonly used technique is to choose the natural variables (NVF), pressures, saturations, mole fractions, then the primary variable switching is applied (Coats, 1980): the variables are changed according to the present phases, that is defined by means of the phase stability test. Different combinations are also employed; e.g., pressure (water), saturation(s) and depending on the system state, liquid or gas mole fractions, S g y i or S l x i with the constraint equations (on saturation and on mole fractions), V or L, z i , ln K i , for an overview of molar and volume formulations see [START_REF] Voskov | Comparison of nonlinear formulations for two-phase multi-component eos based simulation[END_REF].

To avoid the variable switching, several alternative choices of fixed primary variables were proposed, on the example of two-phase {g, l} two-component k = {w, n} or/and compositional problem.

• Gas pressure and total concentrations through the concept of extended saturation:

Abadpour and Panfilov ( 2009) extended the gas saturation out of range [0, 1] that indicates single-phase fluid, the authors demonstrated the equivalence of the twophase and pseudo two-phase formulations.

Remarks: The widely used negative flash method [START_REF] Whitson | The negative flash[END_REF] allows to calculate the gas phase molar fraction. Its value can also be out of range [0,1] and indicates the number of phases.

• Liquid pressure and component mass: [START_REF] Acs | General purpose compositional model[END_REF] proposed a general purpose compositional formulation by choosing one pressure and the component masses. The system of N s components comprises a volume balance equation and the N s mass balance equations.

• Liquid pressure and mass concentration/density of component n in the liquid phase ρ l,n = ρ l X l,n (two-phase two-component problem): the reformulation of solubility or molar fraction of non-wetting component in the liquid phase as a function of pressure using Henry's law was proposed [START_REF] Ippisch | Coupled Transport in Natural Porous Media[END_REF]. [START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF] used the Henry's law to impose a criteria of phase (dis)appearance,

ρ l,n = ρ l,n (p g ) = ρ l,n (p l + p c (S g ))
, so then S g = S g (ρ l,n , p l ) using the inverse capillary function.

• Liquid pressure, liquid saturation and X l,n (two-phase two-component problem):

Jaffré and Sboui (2010) applied the phase transition criteria of [START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF] in the problem system as nonlinear complementarity constraints that implies the third variable.

• Gas pressure and liquid pressure (two-phase two-component problem): Angelini et al.

(2011) chose p g instead of ρ l,n , also making use of the capillary pressure extension and the solubility function.

• Capillary pressure and gas pressure (two-phase two-component problem): Neumann et al. ( 2013) also employed the inverse capillary pressure and the relation x l,n =

x l,n (p g ) but with the nonlinear solubility function.

• Pressures, saturations, fugacities (PSF): [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] proposed the extended but fixed set of primary variables, pressures, saturations, fugacities, with the complementarity condition, solving the problem by the semi-smooth Newton method;

the method avoids the flash calculation.

• Tie-line based (gamma) variables: [START_REF] Entov | On approximation of phase equilibria of multicomponent hydrocarbon mixtures and prediction of oil displacement by gas injection[END_REF] demonstrated that the thermodynamic equilibrium problem could be parametrized by projecting the compositional system on the tie-line space Γ. Using this compositional space parametrization (CSP) approach, [START_REF] Voskov | Compositional space parameterization: theory and application for immiscible displacements[END_REF] introduced a preprocessing step that provides tabulated solutions on the tie-lines of interest. The composition is linearly interpolated within the simplex if it lies close to one of the known tie-lines; otherwise, an equation of state based procedure is applied. The compositional space adaptive tabulation method can be employed in order to avoid the phase stability test and accelerate/replace the flash calculations: the tie-lines are adaptively constructed, the Γ-space is discretized by Delaunay tesseletion and the composition is interpolated inside the simplex. [START_REF] Zaydullin | Nonlinear formulation based on an equationof-state free method for compositional flow simulation[END_REF] have presented an adaptive CSP method, using a new set of N s variables for subcritical region: pressure, liquid phase volume fraction and tie-line parameters γ. These fixed variables, called gamma variables, are switched to pressure and overall composition in the supercritical space. The new gamma variables entail a different interpolation of the composition within the simplex and a modified flash calculations without the Rachford-Rice equation (Smith and van Ness, 1987). Masson et al. (2014b) compared the NVF (Coats, 1980), PSF [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] and extended pressure-pressure formulation (PPF) (Angelini, 2010) formulations on the numerical convergence in the context of phase (dis)appearance. The authors demonstrated the equivalence of formulations and revealed the similarly better results of NVF and PSF (in terms of the convergence) than those of PPF, especially modeling the gas appearance by transport. Ben [START_REF] Ben Gharbia | Study of compositional multi-phase flow formulations with cubic eos[END_REF] presented the NVF and PPF performance on the multicomponent system. The authors showed that the numerical difficulties of PPF appearing in the transition zones caused the time step to decrease and then made it less efficient than the NVF; however, the nature of this issue was not defined in the work.

By applying the flash calculation, described above, one can solve the system (3.1.1) on P and the overall composition z, Eq. (3.1.12), that are consistent during the phase transition. The comparison between the overall composition and NVF was given by [START_REF] Lu | Evaluating variable switching and flash methods in modeling carbon sequestration in deep geologic formations using pflotran[END_REF], where the formulation using the total mole fractions displayed a higher convergence over the NVF formulation for CO 2 injection problem. [START_REF] Voskov | Comparison of nonlinear formulations for two-phase multi-component eos based simulation[END_REF] provided a comparative study of the NVF, molar and volume formulations and illustrated that the NVF is a preferable choice for isothermal immiscible and near-miscible problems of gas displacement in heterogeneous medium, especially when the phase state changes.

The implicit methods, especially AIM, are efficient methods for solving the compositional problem of the advective dominant flow. Despite the changing of the Jacobian structure, the NVF (Coats, 1980) is still one of the best performing and most reliable approaches. However, the PSF and the overall composition formulations have recently demonstrated a good convergence in Masson et al. (2014b) and [START_REF] Lu | Evaluating variable switching and flash methods in modeling carbon sequestration in deep geologic formations using pflotran[END_REF], respectively, and need more investigation.

Reactive transport

The numerical methods of RT (Yeh andTripathi, 1989, 1991;[START_REF] Steefel | Reactive transport in porous media[END_REF][START_REF] Saaltink | On the behavior of approaches to simulate reactive transport[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF] are divided into the operator splitting (OS) and global implicit approach (GIA): the former reduces the degree of freedom by splitting the RT into the transport and chemistry modules and by solving them sequentially (non-)iteratively, while the latter handles the full system implicitly. Omitting the source terms, porosity and variable saturation, the reactive transport system can be presented as

dT dt = T(C), (3.1.17) R(X, T ) = 0, (3.1.18)
where the vectors T and C are the total and mobile total concentration vectors and are supposed to be transposed in the transport block (3.1.17), T is the transposed transport operator vector, R(X, T ) is the chemistry block solved on the unknowns X (e.g., basis species) with the input total concentration vector T . Once X is known, the total and mobile total concentration vectors T and C can be defined, S(X) = C. Next, we point out the concepts of classic approaches from the sequential non-iterative approach (SNIA)

to the global differential algebraic equations (DAE) approach.

SNIA

Assuming that the fluid is first transported and the reactions occur just after the displacement, the transport and chemistry blocks are separately (by OS) and sequentially solved.

The scheme of SNIA for the time step ∆t n+1 is expressed as (3.1.19) where the explicit time discretization is employed for the transport, then the chemistry module is applied for each cell separately. It is well known that the SNIA is associated with the persistent error for the decay problem that was demonstrated by [START_REF] Valocchi | Accuracy of operator splitting for advection-dispersionreaction problems[END_REF]. To increase the accuracy, one can make use of the strang splitting approach, proposed by [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]. The whole problem is divided into three parts:

T tr = T n + ∆t n+1 T(C n ), R(X n+1 , T tr ) = 0,
the advective transport for the half time step, then the chemistry step and finally, the dispersive transport for the second half time step.

SIA

The sequential iterative approach (SIA) allows to preserve the OSA and to solve the system globally by Picard's method. The implicit scheme is applied for the transport and then the chemistry module is used for each cell as in SNIA. Therefore, the nonlinear system arising from the chemical operator is solved inside Picard's method. The SIA scheme can be written as

T n+1,2k+1 = T n + ∆t n+1 T(C n+1,2k+1 ), R(X n+1,2k+2 , T n+1,2k+1 ) = 0. (3.1.20)
The transport and chemistry modules iterate until the convergence, that reduces the errors in comparison with those of SNIA [START_REF] Carrayrou | Operator-splitting procedures for reactive transport and comparison of mass balance errors[END_REF].

DSA

In the global direct substitution approach (DSA), the RT problem is solved for the chemistry variables X, substituted in the transport, that leads to

(T (X)) n+1 = (T (X)) n + ∆t n+1 T((C(X)) n+1 ). (3.1.21)
We are left with the nonlinear system on X. By applying Newton's method, the derivatives of T (X) and C(X) should be calculated, so the chemistry is incorporated in the transport equations. The DSA system can be reduced by the elimination of mineral reactions at equilibrium [START_REF] Lichtner | Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems[END_REF][START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF].

DAE

The RT system Eqs. (3.1.17) and (3.1.18) can be generalized as a system of DAE:

A dY dt = B(Y ). (3.1.22)
A is the block matrix, where the non-zero elements form 2009). The RT system (3.1.22) embraces both the transport and the chemistry modules and is then aimed at making use of the accuracy properties of the DAE solvers. Newton's method is globally applied. The size of the final linear system is at least twice as large as that of the OSA and of the DSA, but it is sparser than that of DSA. [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the momas benchmark case[END_REF] displayed the accuracy advantages of DAE solvers.

A 11 = I, dim A 11 = dim T , Y = (T, X, C) T is the unknown vector, B = (T(C), R(X, T ), S(X) -C) T , that is detailed in de Dieuleveult et al. (
All presented above methods in Sec. 3.1.2 have pros and cons. The SNIA is computationally efficient due to its structure, although it can provide the consistent mass balance errors, depending on the problem (see above). The SIA retains the modular framework and decreases the splitting errors by iterative procedure; that makes it more computationally expensive than the SNIA by definition and than the DSA with the reduced system for certain problems [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the momas benchmark case[END_REF]. The DSA proposes a more robust method at the cost of solving the nonlinear transport system. Using reduction techniques allows to minimize the system. The DAE method entails the largest system solution for which the DAE solvers are supposed to be employed.

Existing approaches

A list of decoupled and sequential multiphase multicomponent flow (MMF) and reactive transport (RT), provided below, illustrates the diversity of solution schemes.

• CodeBright [START_REF] Olivella | Numerical formulation for a simulator (code bright) for the coupled analysis of saline media[END_REF] represents a mass balance for water, air and salt species. This simulator is coupled with solved by DSA/SIA RT, code CHEPROO, [START_REF] Bea | Cheproo: A fortran 90 objectoriented module to solve chemical processes in earth science models[END_REF][START_REF] Saaltink | Retraso, a code for modeling reactive transport in saturated and unsaturated porous media[END_REF]. • HYDROGEOCHEM, based on FEM, applies a "strong coupling" which consists of iterating saturated/unsaturated (modified Richard's equation) flow, heat and biogeochemical transport until convergence is reached [START_REF] Yeh | HYDROGEOCHEM 5.0: A three-dimensional model of coupled fluid flow, thermal transport, and hydrogeochemical transport through variably saturated conditions[END_REF][START_REF] Yeh | Hydrogeochem: A coupled model of variably saturated flow, thermal transport, and reactive biogeochemical transport[END_REF]. HY-DROBIOGEOCHEM provides the reactive transport under multiphase flow where the transport equation, called modified hydrologic transport, is adapted to FEM, [START_REF] Fang | Reactive Chemical transport under multiphase system[END_REF]. However, there is no coupling with the multiphase flow.

•
• iCP [START_REF] Nardi | Interface comsol-phreeqc (icp), an e cient numerical framework for the solution of coupled multiphysics and geochemistry[END_REF] is an interface between the multiphysic simulator COMSOL, based on FEM, and the geochemical code PHREEQC [START_REF] Parkhurst | User's guide to PHREEQC (Version 2): A computer program for speciation, batch-reaction, onedimensional transport, and inverse geochemical calculations[END_REF]. The (un)saturated flow step is provided by COMSOL, then RT is handled by SNIA where the chemistry is solved in each cell by PHREEQC. The communication between two codes (transfer of updated variables) is user-defined, thus the consistency of coupling is variable. The coupling with two-phase flow is also developed 1 .

• IPARS [START_REF] Wheeler | Modeling of flow and reactive transport in IPARS[END_REF] The species are divided into those in the flowing phases and in the stationary phases.

The fluid and rock properties are composition independent during the RT, nevertheless the fluids and rock can be compressible. During the TRCHEM step, the phase mass transfer does not impact on the flow.

In [START_REF] Peszynska | Reactive transport model coupled to multiphase flow models[END_REF], the multiphase flow formulation was presented under the hypothesis of slightly compressible fluids and of density and composition independent flow, loosely coupled with the SNI reactive transport divided into advection, reaction and diffusion steps.

• MIN3P [START_REF] Mayer | Reactive transport modeling in variably saturated media with min3p: Basic model formulation and model enhancements[END_REF] combines (un)saturated flow with the RT solved by the direct substitution approach (DSA). MIN3P-Bubble involves the sequential iterative solution of (1) unsaturated flow, ( 2) reactive transport and (3) gas bubble formation step that comprises the objective function on the trapped gas saturation, considering that the gas concentrations are expressed by partial pressures through the hypothesis of ideal gas, and finally (4) the relative permeability update [START_REF] Amos | Investigating the role of gas bubble formation and entrapment in contaminated aquifers: Reactive transport modelling[END_REF].

Furthermore, MIN3P-Dusty embodies the reaction-driven advection, hence, the gas velocity term appears in the transport [START_REF] Molins | Coupling between geochemical reactions and multicomponent gas and solute transport in unsaturated media: A reactive transport modeling study[END_REF].

• MoReS is an EOS and K-value compositional simulator, based on FVM, FIM, IM-PES, adaptive implicit method (AIM) are available. The full SNIA of multiphase flow, transport for primary (aqueous) species and chemistry modules is proposed in [START_REF] Farajzadeh | Detailed modeling of the alkali/surfactant/polymer (asp) process by coupling a multipurpose reservoir simulator to the chemistry package phreeqc[END_REF] and Wei (2012). In the flow step, the pressure equations and mass conservation for each component are solved, where the components may be hydrocarbons or the total dissolved solids (TDS). The batch reactions are modeled by PHREEQC. The fluid properties depend on the pressure, temperature, the amount of dissolved gases, TDS and salinity. This method is limited to slow chemical reactions and requires small time steps (Wei, 2012).

• NUFT sequentially solves compositional multiphase flow and then the GI system of basis species transport in equilibrium combined with mass balance equations of kinetic reactions (TRANS module) [START_REF] Hao | Overview of nuft: a versatile numerical model for simulating flow and reactive transport in porous media[END_REF].

• PFLOTRAN considers a sequential coupling of compositional flow (water/CO 2 ) in MPHASE mode and reactive transport (CHEMISTRY) by applying fully implicit solution of each module [START_REF] Lichtner | Pflotran user manual[END_REF][START_REF] Lu | High resolution numerical investigation on the e↵ect of convective instability on long term co2 storage in saline aquifers[END_REF].

• STOMP chooses a "strong coupling" by iterating flow and SI loop of transport and chemistry. The detailed description is given in [START_REF] White | STOMP Subsurface Transport Over Multiple Phases, Version 4.0, User's Guide[END_REF] • TOUGHREACT [START_REF] Xu | Coupled modeling of non-isothermal multiphase flow, solute transport and reactive chemistry in porous and fractured media: 1. model development and validation[END_REF][START_REF] Xu | TOUGHREACT user's guide: A simulation program for non-isothermal multiphase reactive transport in variably saturated geologic media[END_REF] introduces the geochemical modeling in TOUGH2 [START_REF] Pruess | TOUGH2 User's Guide[END_REF] by coupling sequentially two-component (water/air) two-phase flow and SI RT. TOUGH2 allows to model different systems (e.g., water/hydrogen) by applying the corresponding EOS module.

• UTCHEM is a chemical-flooding multiphase compositional simulator [START_REF] Delshad | Utchem version 9.0 technical documentation[END_REF], considering the equilibrium reactions and the polymer/gel and biodegradation kinetics. The compositional flow is carried out by a decoupled formulation which consists of the pressure equation and mass conservation for each component (IMPEC).

Then the phase equilibria calculations are performed for the concentrations and saturations and EQBATCH for equilibrium reactions in batch (under hypothesis of ideal aqueous solution).

The following codes are based on GIA to solve fully coupled compositional flow and geochemical reactions.

• AD-GPRS, that initially solved a compositional formulation with molar/natural primary variables by FIM/AIM (Cao, 2002;[START_REF] Zhou | Automatic differentiation framework for compositional simulation on unstructured grids with multi-point discretization schemes[END_REF], has been extended to geochemical modeling by means of the element-based method [START_REF] Fan | A fully-coupled flowreactive-transport formulation based on element conservation, with application to co 2 storage simulations[END_REF], including the phase equilibria and flash calculations. Both equilibrium and kinetic reactions can be handled.

• COORES implicitly solves the compositional flow by means of the active set method.

The active set of variables is established at each iteration, the phase transition is given by flash calculation. The geochemical block is built as a separate module. • GEM-GHG provides EOS geochemical compositional modeling [START_REF] Nghiem | Modeling co2 storage in aquifers with a fully-coupled geochemical eos compositional simulator[END_REF].

The AIM is employed to reduce the size of nonlinear system. After eliminating the explicit chemical equilibrium rates from mass conservation equations, the set is solved by Newton's method. The neighbour independent processes are also decoupled.

Mineral dissolution and precipitation can only be performed as kinetic reactions.

• GPAS initially used FIM for solving EOS compositional flow and was developed for chemical flooding [START_REF] Pope | A new generation chemical flooding simulator[END_REF]. The chemistry modeling capabilities were added similar to those of UTCHEM. The chemical species are present only in the aqueous phase. Two approaches are performed: hybrid and fully implicit. The former proposes to solve implicit balance equations for hydrocarbon (like dominant components), then explicit equations for the minor components.

CHAPTER 3. COUPLING MMF AND RT

• STARS applies the different level of implicitness (AIM) to solve the k-value based compositional flow with geochemical reactions, developed by CMG as well as GEM.

Mineral reactions are considered kinetic, similar to GEM (Computer Modelling Group, 2009).

Discussion

Both the RT and the MMF formulations had been extensively studied, the development of MMF and RT coupling started later, inheriting the same OSA and GIA as in RT. In the OSA, the general idea is to provide fluxes, pressures and saturations calculated by MMF to RT, similar to the single phase RT modeling. The flow is supposed to be miscible, so the compositional formulation is used, that implies solving the nonlinear system, whose dimension evidently depends on the number of components. Therefore, the efficiency of the method relies on the linear solvers and preconditionners (or/and on the multiscale methods) as well as for the GIA. New methods continue to appear, introducing different levels of coupling between flow and RT, from the weak sequential to the strong SIA.

Many variations of coupling approaches exist (Sec. 3.1.3), some works give the method comparisons, suggesting that the coupling issue still needs to be investigated. The detailed conclusion on the methods listed in Sec. 3.1.3, is given in the next section, where we present a new method of integrating the MMF in the RT simulator. In contrast to the existing approaches, the phase formulation is employed to model the multiphase multicomponent system by means of our coupling method (Sec. 3.2), that allows to preserve the flow structure independent of the number of components.

Integrating a compressible multicomponent two-phase flow into existing reactive transport simulator

Abstract

This work aims at incorporating compressible multiphase flow into the conventional reactive transport framework by an operator splitting approach.

A new approach allows to retain the general paradigm of the flow module independent of the geochemical processes and to model complex multiphase chemical systems, conserving the versatile structure of conventional reactive transport. The number of mass conservation nonlinear equations arising from the flow module is minimum, the phase flow formulation is employed. Applying the appropriate equation of state made it possible to precisely describe the compressible multicomponent phases, their thermodynamic properties and relevant fluxes. The proposed flow coupling method was implemented in the reactive transport software HYTEC. The whole framework preserves its flexibility in the further physical and numerical purpose developments. The verification of the coupling is obtained by modeling a problem with a selfsimilar solution. The simulation of a 2D CO 2 injection problem shows the pertinent physical results and the computational e ciency of the method.

1. Introduction

Background/Motivation

Human activity in the subsurface has been expanding and diversifying (waste disposal, mining excavation, high-frequency storage of energy), meanwhile the public and regulatory expectations have been growing. Assessment of each step of underground operations including environmental impact evaluation, relies on elaborate simulators and gives rise to a strong necessity of developing the multiphysics modeling. Reactive transport, a geochemical research and engineering tool, deals with multicomponent systems and sophisticated chemical processes (activity and fugacity correction according to di↵erent models, mineral dissolution and precipitation, cation exchange, oxidation and reduction reactions, isotopic fractionation and filiation), in addition to gas evaporation and dissolution [START_REF] Mayer | Reactive transport modeling in variably saturated media with min3p: Basic model formulation and model enhancements[END_REF][START_REF] Parkhurst | User's guide to PHREEQC (Version 2): A computer program for speciation, batch-reaction, onedimensional transport, and inverse geochemical calculations[END_REF][START_REF] Steefel | Crunchflow software for modeling multicomponent reactive flow and transport. user's manual[END_REF][START_REF] Van Der Lee | Moduleoriented modeling of reactive transport with hytec[END_REF][START_REF] Yeh | HYDROGEOCHEM 5.0: A three-dimensional model of coupled fluid flow, thermal transport, and hydrogeochemical transport through variably saturated conditions[END_REF]. Multiphase flow inherits a strong experience from reservoir engineering researches including thermodynamic modeling of complex phase behavior. In particular, the equations of state were deployed to simulate and study interfacial tension, gas, steam or alkaline injection in oil reservoir and enhanced oil recovery [26,28,[START_REF] Nghiem | Modeling co2 storage in aquifers with a fully-coupled geochemical eos compositional simulator[END_REF][START_REF] Wang | A new generation eos compositional reservoir simulator: Part i-formulation and discretization[END_REF].

This work aims at incorporating a compressible multiphase flow module into an existing reactive transport simulator. Our coupling method should therefore meet the following requirements:

1. the new approach should handle the di↵erent complex multiphase chemical models and retain the general paradigm of a multiphase flow module independent of the geochemical system conserving the conventional reactive transport structure; 2. the number of mass conservation nonlinear equations arising from the flow module should be minimum such that the reduced flow system preserves the matrix structure in order to minimize the computational intensity; 3. the entire framework should preserve its flexibility towards possible further non-isothermal, geomechanical and domain decomposition developments.

The subject of reactive transport methods was extensively investigated in the last two decades (see below). This work focuses on the coupling between multicomponent multiphase flow (MMF1 ) and reactive transport (RT) modules, starting from the survey of the existing approaches.

A review of multiphase multicomponent flow and reactive transport codes 1.2.1. Operator Splitting Algorithms between MMF and RT

The strength of the operator-splitting (OS) approach (sequential iterative SIA, or sequential non-iterative SNIA) stems from the framework flexibility allowing to develop and verify each module independently. These are major reasons to choose the OS approach for the coupling between MMF and RT, especially when a hydro-geochemical code is aimed at the extension from the single to two phase flow. The following codes apply the OS approach: CodeBright [START_REF] Olivella | Numerical formulation for a simulator (code bright) for the coupled analysis of saline media[END_REF], DuMu X (based on DUNE) [3,[START_REF] Vostrikov | Numerical simulation of two-phase multicomponent flow with reactive transport in porous media[END_REF], DUNE [33], HYDRO-GEOCHEM (unsaturated) [START_REF] Yeh | HYDROGEOCHEM 5.0: A three-dimensional model of coupled fluid flow, thermal transport, and hydrogeochemical transport through variably saturated conditions[END_REF][START_REF] Yeh | Hydrogeochem: A coupled model of variably saturated flow, thermal transport, and reactive biogeochemical transport[END_REF], iCP [START_REF] Nardi | Interface comsol-phreeqc (icp), an e cient numerical framework for the solution of coupled multiphysics and geochemistry[END_REF], IPARS [START_REF] Peszynska | Reactive transport model coupled to multiphase flow models[END_REF][START_REF] Wheeler | Modeling of flow and reactive transport in IPARS[END_REF], MIN3P (the bubble model) [START_REF] Mayer | Reactive transport modeling in variably saturated media with min3p: Basic model formulation and model enhancements[END_REF][START_REF] Molins | Coupling between geochemical reactions and multicomponent gas and solute transport in unsaturated media: A reactive transport modeling study[END_REF], MoReS [28,[START_REF] Wei | Sequential coupling of geochemical reactions with reservoir simulations for waterflood and eor studies[END_REF], NUFT [31], PFLOTRAN [START_REF] Lichtner | Pflotran user manual[END_REF][START_REF] Lu | High resolution numerical investigation on the e↵ect of convective instability on long term co2 storage in saline aquifers[END_REF], STOMP [START_REF] White | STOMP Subsurface Transport Over Multiple Phases, Version 4.0, User's Guide[END_REF][START_REF] White | STOMP Subsurface Transport Over Multiple Phases: STOMP-CO2 and STOMP-CO2e Guide: Version 1.0[END_REF], TOUGHREACT [START_REF] Xu | Coupled modeling of non-isothermal multiphase flow, solute transport and reactive chemistry in porous and fractured media: 1. model development and validation[END_REF][START_REF] Xu | TOUGHREACT user's guide: A simulation program for non-isothermal multiphase reactive transport in variably saturated geologic media[END_REF], UTCHEM [26], also see [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF].

The simulators are chiefly based on the finite volume method (FVM) due to its conservative properties. The general tendency of coupling is first to solve (non-)iteratively flow that provides the velocities. The compositional formulation is usually chosen or alternatively, one solves the conservation equations for the dominant components (e.g. water/air). When the twophase flow system involves two components only (e.g. {H 2 O, CO 2 }), the final set of equations is obviously identical and reduced to the mass conservation for each component.

Once the flow is established, the RT part can be solved by the OS (SIA, SNIA, predictor-corrector) or global implicit approach (GIA) such as the ordinary di↵erential equations based method (ODE, the chemistry module is used as a black box) or the direct substitution approach (DSA) [START_REF] Saaltink | On the behavior of approaches to simulate reactive transport[END_REF], or the di↵erential algebraic equations based method (DAE) [20], [19]. The comparison of OS and GI approaches applied to RT were given in [12], [20], [START_REF] Saaltink | On the behavior of approaches to simulate reactive transport[END_REF], [START_REF] Steefel | Reactive transport in porous media[END_REF]. In 2001, Saaltink et al. exhibited that SIA was more favorable in case of large grid and low kinetic rates than DSA due to the computer storage increase and slow linear solvers. A decade later, the MoMaS benchmark displayed the reliability of both approaches and the enhanced computational potential of the DSA with the reduction technique [38].

1.2.2. From Global Implicit Algorithm in RT to Global Implicit Algorithm in MMF&RT The advantage of GIA is its accuracy at the cost of computational resources. However, this method is becoming more competitive with the increasing computer capabilities and with the advances in the methods that allow to reduce the system of equations. Also, reservoir simulators typically make use of the global formulation for the MMF problem in combination with the fully implicit method (FIM) or the adaptive implicit one (AIM), that unifies the FIM and implicit pressure/explicit concentration (IMPEC). They are therefore naturally extended to the RT using GIA, e.g. COORES 2 , GEM-GHG [START_REF] Nghiem | Modeling co2 storage in aquifers with a fully-coupled geochemical eos compositional simulator[END_REF], GPAS [START_REF] Pope | A new generation chemical flooding simulator[END_REF], GPRS [11,27].

For the global implicit solution of the RT, several techniques propose to reduce the initial mass balance system by linear combination and to eliminate the reaction terms in the equilibrium reactions [38,[START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF][START_REF] Molins | A formulation for decoupling components in reactive transport problems[END_REF]. Such modifications of the DSA for RT lead to mathematical decoupling of the whole system and consequently make it interesting for global implicit coupling of multiphase flow and reactive transport as was recently demonstrated in works of Saaltink et al. [71] and of Fan et al. [27]. Saaltink et al. proposed a method to introduce the chemistry calculations in conventional multiphase simulator with the idea to keep a minimum number of mass conservation equations taking the fluid phase pressures and porosity as primary variables and expressing all secondary variables as a polynomial function of gas pressure (concentrations of the components, fugacity, pH, salinity etc), that requires a pre-processing to use it for each application. Fan et al. employed the element balance formulation [START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF] by means of the reduction techniques [38,[START_REF] Molins | A formulation for decoupling components in reactive transport problems[END_REF] and the decoupled linearized system on the primary and secondary equations. Saaltink et al. demonstrated that the GIA and OSA gave similar results in the context of CO 2 storage and concluded that the full coupling was not necessary for MMF&RT modeling, from a physical point of view (although the OS can be more computationally expensive).

Meanwhile, Gamazo et al. [29] highlighted the significant impact of geochemical reactions on the phase fluxes by modeling the gypsum dehydration when the non-isothermal flow is chemically restrained, i.e. the anhydritegypsum paragenesis controls the water activity and hence the evaporation process. The importance of DSA, implicitly connecting the equilibrium het-erogeneous reactions and the phase flow, was emphasized in comparison with the decoupled formulation of the flow and the RT: the decoupled formulation overestimates the evaporation, whereas its computational time was reduced by 22.5% (compared to DSA). However, the definition for water activity and liquid density were di↵erent for the GI and the OS models. Furthermore, the decoupled flow system comprised the conservation equations for dominant components (water and air), but the coupling between flow and RT was not detailed. Given the nature of the formulations described in this work, it appears likely that using the strong OSA (the SI connection between flow and RT with the precise reaction terms in flow equations) would yield results similar to those of DSA at the cost of additional iterations.

Preamble to a new approach of MMF&RT coupling

Only dominant components can be considered (the dominant components flow formulation) in order to reduce the number of nonlinear equations arising form the flow system. This is e cient when the impact of the other species are negligible. In geochemical problems, the speciation can vary largely over time and space. Ideally, the dominant components should be adapted locally to preserve the accuracy, but that would require specific treatment for the global flow solution.

If the flow system expands such that many species play significant role in the thermodynamic state and phase displacement, the number of nonlinear equations of compositional flow increases with the number of components, that is computationally expensive.

Despite the reduction techniques, the calculation of the primary nonlinear system is unavoidable. In the geochemical modeling, which encompasses the complex homogeneous and heterogeneous reactions, the primary species usually changes both in space and in time. The system is therefore redefined entailing the modification of Jacobian structure inside Newton's method. Meanwhile, the solution of the transport linear system can be several or even hundred times faster than that of multiphase flow.

Can we then replace the compositional formulation by the phase formulation in order to retain a minimal size of non-linear system? The phase flow coupled to the RT was employed in the work of Peszynska and Sun. Given their problem conditions (slightly compressible fluids and density independent flow), the authors divided the RT into three stages (advection, reaction, di↵usion) within the internal time steps by interpolating the fluxes and saturations. Later, Hron et al. simulated Escherichia coli growth and transport under aerobic and anaerobic conditions propounding the sequential coupling of phase flow and the S(N)I RT. The liquid fluid was supposed incompressible, while the gas density depended on phase composition. For advective dominant regime, the transport was solved explicitly in time that led to the time step restriction; CFL = 0.4 and the split transport time step were thus applied. Notice that the same issue arises from the IMPEC scheme where the component equations are solved explicitly.

Our alternative method

This work proposes a new approach for incorporating the compressible two-phase flow in a conventional reactive transport simulator. It employs the phase flow formulation and preserves all the sustainability and facilities of reactive part. The method inherits the unconditional stability of the fully implicit scheme and can handle the modeling of advective dominant and density driven regimes. The following simplifications are considered in this work: isothermal flow, no geomechanics. The method was applied to the (SIA) reactive transport simulator, HYTEC [START_REF] Lagneau | Operator-splitting-based reactive transport models in strong feedback of porosity change: The contribution of analytical solutions for accuracy validation and estimator improvement[END_REF][START_REF] Van Der Lee | Moduleoriented modeling of reactive transport with hytec[END_REF], that was widely evaluated in several benchmarks [22,12,[START_REF] Lagneau | Hytec results of the momas reactive transport benchmark[END_REF][START_REF] Trotignon | Design of a 2-d cementation experiment in porous medium using numerical simulation[END_REF] and in numerous applications, such as cement degradation [24], radioactive waste disposal [25,23], geological storage of acid gases [18,35,39], uranium in situ recovery processes [START_REF] Regnault | 3d reactive transport simulations of uranium in situ leaching : Forecast and process optimization[END_REF].

In Sec. 2, a concise description of the governing equations of multicomponent multiphase flow and reactive transport is given. Next, the proposed methods of coupling are detailed in Sec. 3. Sec. 4 then demonstrates the method's applicability and computational performance first on a benchmark problem with a self-similar solution and then, on a 2D CO 2 injection modeling. The conclusions and discussions are finally made in Sec. 5.

Governing processes and mathematical formulation

The entire isothermal MMRF problem is composed of phase ↵ mass and component c ↵,k mol conservation, mass balance, mass-action laws and other constitutive relations. The N f fluid phases ↵, liquid and gas N f = 2, are mobile, while solid phase is supposed to be immobile, N p = 3. The chemical system of N s chemical species and N r linearly independent chemical reactions relies on Morel's method [START_REF] Morel | Principles and Applications of Aquatic Chemistry[END_REF] of N c primary species that form a basis of all N s species: N c = N s N r .

Mass conservation for each phase

We firstly introduce the standard multiphase compositional flow problem, the physical parameters and then present the alternative and its advantages.

Let us begin with the general mass conservation in terms of mass fraction X ↵,k of species k in fluid phase ↵ that forms N c N f equations:

@( S ↵ ⇢ ↵ X ↵,k ) @t + r • (⇢ ↵ X ↵,k u ↵ ⇢ ↵ D ↵,k rX ↵,k ) = R ↵,k + q ↵,k , (1)
where ! is the porosity, S ↵ is the saturation of fluid phase ↵, ⇢ ↵ is the mass density of fluid phase ↵, u ↵ is the Darcy-Muskat velocity of fluid phase ↵, D ↵,k is the dispersion tensor of fluid phase ↵, R ↵,k is the reaction term and q ↵,k is the external source term of species k in fluid phase ↵. The Darcy-Muskat law [START_REF] Muskat | Flow of gas-liquid mixtures through sands[END_REF] states:

u ↵ = k r↵ µ ↵ K (rp ↵ ⇢ ↵ g) , (2)
where k r↵ and µ ↵ are the relative permeability and the viscosity of fluid phase ↵, respectively, K is the intrinsic permeability tensor, p ↵ is the pressure of fluid phase ↵, g is the gravitational acceleration vector. The pressures are connected by N f 1 capillary pressure relations, the saturations and mass fractions sum to 1:

X ↵ S ↵ = 1, (3)
X k X ↵,k = 1, (4)
that results in additional 2N f constitutive relations. The conventional PDE system of isothermal compositional multiphase flow can be expressed from ( 1) by deriving a conservation of each species k in all phases P ↵ . These N c nonlinear equations, 2N f constitutive relations and N c (N f 1) phase equilibrium relations (sec. 2.3) give 2N f + N c N f equations for 2N f + N c N f unknowns p ↵ , S ↵ , X ↵,k . According to the Gibbs phase rule, the number of intensive properties is N c N p +2. So, if the number of phases is locally known, N p = 2 for fixed temperature, at least N c nonlinear equations (for pressure and composition) must be solved in order to establish the thermodynamic state. The geochemical system can be abundant in species that implies the solution of considerably large nonlinear systems.

In this work, we propose to handle the phase mass conservation system of N f PDE whose size is independent of the number of chemical species N s . We apply

P k to (1) to obtain @( S ↵ ⇢ ↵ ) @t + r • (⇢ ↵ u ↵ ) = R ↵ + q ↵ , (5)
taking into account that the sum of di↵usive fluxes for each phase ↵ is equal to 0. The external source term q ↵ can be also presented as

q ↵ = ⇢ ↵ ↵ , ( 6
)
where ↵ is the source velocity. The system (5) of N f equations, N f 1 capillary pressure relations and ( 3) are assembled for 2N f unknowns p ↵ , S ↵ which are natural variables of multiphase flow problem.

When one of the phases disappears, the corresponding equation of Sys. ( 5) degenerates and the natural variables are inappropriate to describe the system. In this case, we pass to the single flow problem. Numerous formalisms are dedicated to the two-component two-phase flow and associated phase disappearance/appearance [1,6,10,[START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF][START_REF] Pruess | TOUGH2 User's Guide[END_REF]. However, the liquid phase can be supposed to be present throughout the system, even if it remains in small amounts. Therefore, solving the system of N f nonlinear equations ( 5) can be beneficial when N c > N f .

Mole conservation for each gas component and primary species

The liquid and gas phases consist of N c primary species and N g gas species, respectively, N g < N c , for which the transport must be solved. We chose the transport formulation in terms of concentrations c ↵,k / ⇢ ↵ X ↵,k /M k . Deriving the transport equations in mole/mass fractions as in Eq. ( 1), one neglects the density deviation r⇢ ↵ /⇢ ↵ that arises from the di↵usive part of flux. Based on the primary species formalism [START_REF] Lichtner | Reactive transport in porous media[END_REF][START_REF] Steefel | Reactive transport in porous media[END_REF][START_REF] Van Der Lee | Thermodynamic and mathematical concepts of chess[END_REF][START_REF] Yeh | A model for simulating transport of reactive multi-species components: model development and demonstration[END_REF], the liquid transport is defined for total liquid mobile concentration c l,k of primary species k and the gas transport for gas concentration c g,m of gas species m, that evolves into N c + N g transport (linear) equations:

@ S l c l,k @t = T l (c l,k ) + R l,k (c l,k , c s,k ) + q l,k , (7)
@ S g c g,m @t = T g (c g,m ) + R g,m (c g,m ) + q g,m , (8)
where c s,k is the immobile concentration of species k, operator T ↵ includes the advective flux presented by Darcy-Muskat law (2) and, supposing non-Knudsen di↵usion, Fick's law for di↵usive-dispersive flux that yields [21,[START_REF] Lichtner | Reactive transport in porous media[END_REF]:

T ↵ (c ↵,k ) = r • (D ↵,k rc ↵,k c ↵,k u ↵ ), (9)
where D ↵,k involves the molecular di↵usion D ↵,k , the tortuosity ⌧ ↵ [START_REF] Millington | Permeability of porous solids[END_REF] and the longitudinal and transverse dispersion.

Mole balance for each primary basis species

Here, we introduce concisely the reactive part presented by CHESS [START_REF] Van Der Lee | Thermodynamic and mathematical concepts of chess[END_REF], that yields the mass action law and phase equilibrium relations. As mentioned above, the reactive transport code is based on the primary species formulation, so we denote the concentration of species S j , j = 1, . . . , N s , that can be expressed as a function of basis species C i , i = 1, . . . , N c :

S j ⌦ Nc X i=1 ↵ ij C i , (10)
where ↵ ij is the stoichiometric coe cient. At equilibrium, the mass action laws provide the reaction a nity:

S j = K j j Nc Y i=1 ( i C i ) ↵ ij , ( 11
)
where j is the activity coe cient and K j is the thermodynamic equilibrium constant of reaction j. All aqueous reactions are treated at equilibrium. For each basis species, the mole balances can be written by the total concentration

T i , i = 1, .., N c , T i = P Ns j=1 ↵ ji S j , T i C i Ns X j=1,j6 =i ↵ ji K j j Nc Y i=1 ( i C i ) ↵ ij = 0, (12)
that constitute a system of N c equations on C i , i = 1, .., N c .

Liquid mixtures

In the non-ideal liquid mixtures, the activity coe cients i are not trivial and can be calculated by di↵erent models whose complexity increases with the concentration of solution, from less to more concentrated solution: truncated Davies formula [17], B-dot [32], SIT [30], Pitzer [START_REF] Pitzer | Ion interaction approach: theory and data correlation[END_REF].

Gas-liquid equilibrium

The chemical potentials and also the corresponding fugacities of species i in mixture f ↵ i are equal under equilibrium conditions. Considering the fugacity-activity (' ) approach, it yields

P y i ' g i = f i = K h i i x i , (13)
where y i is the mole fractions of species i in gas, ' g i is the fugacity coe cient of species i, K h i = K h i (T, P ) is the corrected Henry's constant of species i, i is the asymmetric activity coe cient of species i [START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF]. K h i involves the Poynting factor which corrects the reference fugacity in regard to the pressure and which is near the unity at low to moderate pressures. The fugacity coe cients can be calculated by cubic EOS, we use Peng-Robinson EOS [START_REF] Robinson | The characterization of the heptanes and heavier fractions for the gpa peng-robinson programs[END_REF] which precisely reflects the gas mixture properties for high pressure/temperatures. By given P , the cubic equation should be solved for the compressibility factor Z = P v/RT , where v is the molar volume. The fugacity coe cient is then calculated [START_REF] Robinson | The characterization of the heptanes and heavier fractions for the gpa peng-robinson programs[END_REF], ' g i = ' g i (T, P, Z, T c , P c , Z c , ⌦, ), where T c , P c , Z c are sets of the critical temperature, critical pressure and compressibility factor of species in mixture, respectively, ⌦ is the acentric factor set, is the matrix of empirical binary interaction parameters for each pair of species in mixture.

The EOS also provides the mass density for gas mixtures:

⇢ g = P Ng i=1 y i M i v = M v , (14)
and by analogy for liquid mixtures with the M = P Nc i=1 x i M i . It is worth mentioning that Raoult's law (P y i = P sat x i ) is derived form Eq. ( 13), supposing the low pressure condition and ideality of solution.

Liquid-solid equilibrium and kinetic relations

Considering that the activity of mineral equals unity and omitting index j, one derives from the mass action laws (11) the solubility product constant K s,j for solid at equilibrium

K s = Nc Y i=1 ( i C i ) ↵ ij = 1/K. ( 15
)
However, to describe the non-equilibrium state of mineral, the ion activity product Q s should be defined similarly to K s but for the actual activities.

Next, the saturation index SI yields the state of saturation by definition:

SI = log ✓ Q s K s ◆ = 8 > < > : < 0 undersaturated ! dissolution, = 0 saturated ! equilibrium, > 0 oversaturated ! precipitation. ( 16
)
The kinetics of mineral precipitation/dissolution is usually modeled by the transition state theory rate laws [START_REF] Lasaga | Chemical kinetics of water-rock interactions[END_REF]:

dS dt = A s k ✓ sign(SI) ✓✓ Q s K s ◆ x 1 ◆◆ y Y k a n k k ( 17
)
where A s is the specific surface area, k is the kinetic dissolution/precipitation rate constant, K s is the equilibrium solubility constant for solid S, x, y, n k are the fitting parameters, a k is the activity of the potential catalyzing or inhibiting species.

The mole balance equations with the mass action laws for the species at equilibrium, the rate laws for kinetically limited solids and phase equilibrium relations constitute a complete system of algebro-di↵erential equations whose solution yields the concentrations of basis and then secondary species. Other chemical reactions can also be handled by the formalism of basis species, e.g. cation exchange, surface and organic complexation.

Numerical solution

Applying OS concept in subsurface environmental modeling permits the independent development of separated parts of code and rigourous solution of each of them, that determines the majority choice of RT simulators [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF]. The assessment of di↵erent coupling methods was studied on the MoMaS benchmark of RT codes [13,12], during which the reliability both of SI and of GI approaches was revealed, the detailed results of HYTEC was demonstrated in [START_REF] Lagneau | Hytec results of the momas reactive transport benchmark[END_REF]. Following the general strategy of integrating the (un)saturated flow in RT codes by OS, we solve the compressible two-phase flow block first, that involves the flow system, the gas transport equations and the EOS with the fluid properties models. Then the reactive transport coupling is applied. We propose to employ SIA for each of them, and hence there are, indeed, two internal SI couplings. Let us describe the applied methods for flow and transport discretization and the coupling methods afterward.

Discretization of flow and transport

The discretization of mass phase conservation ( 5) and mole species transport ( 7), ( 8) is built on a Voronoi-type finite volume method. The time approximation of two-phase flow ( 5) is fully implicit, the fluxes are handled by TPFA. It was proven since the 1960s, that the transmissibility discretization a↵ects the numerical stability and accuracy [4,9,[START_REF] Settari | Treatment of nonlinear terms in the numerical solution of partial di↵erential equations for multiphase flow in porous media[END_REF]; in this work, the interface coe cients of flow between adjacent cells are then evaluated implicitly and upstream. For the relative permeability (k r↵ ) ij , the upstream space approximation is widely used, since its convergence was demonstrated on the Buckley-Leverett problem [7,8]. A detailed comparison of temporal discretization presented in [7,9] showed stability advantages of the implicit upstream treatment against explicit ones but increase of truncation errors at the same time. We will see the impact of truncation errors in Section 4.1. So, the relative permeability k r↵ , intrinsic permeability K, phase density ⇢ ↵ , phase viscosity µ ↵ for interface coe cient between the adjacent cells i, j at iteration k + 1 are defined as

(•) k+1 ij = ( (•) k i if u k ↵ > 0, (•) k j else ( 18
)
where u ↵ = u ↵ n ij , n ij is the normal vector from i to j.

The discretization of density (⇢ g ↵ ) ij in the gravity term ⇢ ↵ g needs to be treated di↵erently: it is weighted with regard to the e↵ective phase volume. If the gas phase is absent in one of the cells, the upstream treatment is applied [16]:

(⇢ g ↵ ) ij = 8 > < > : ⇢ ↵,i + ⇢ ↵,j (1 ) if (S ↵,i > 0) ^(S ↵,j > 0) : = V ↵,i /(V ↵,i + V ↵,j ); (⇢ ↵ ) ij else ( 19
)
where V ↵ = S ↵ V tot , V tot is the volume of the cell. The resulting nonlinear system is solved by Newton's method with an analytical Jacobian, that will be displayed below in Sec. 3.2.

The space discretization of transport operators ( 7), ( 8) is performed by the upstream weighting for advective flux and harmonic for e↵ective di↵usion, while one can apply the implicit or semi-implicit method in time. Applying a semi-implicit method for time discretization makes it possible to attain an accurate solution: implicit Euler scheme in time for di↵usive-dispersive and Crank-Nicholson for advective part.

Coupling 1: Compressible two-phase flow

Since the phase flow system ( 5) is nonlinear, the classic Newton's method is applied. We denote the discretized equations (5) as F (x) = 0, where x = (p l , S g ). When the fluids are highly compressible (e.g. gas phase), the density properties should be precisely evaluated at each deviation of intensive variables, P, V, n, where n is the quantity of matter. The gas appearance and disappearance involve the density update in Newton's loop similar to the other flow parameters to ensure the implicit treatment of interface coe cients. However, the gas density can be strongly dependent on its composition. Therefore, the gas composition has to be calculated by employing the gas transport (8), denoted by T g (c g ; x) = 0. The flow coupling algorithm for time step n + 1 is presented in 1 and below its parsing.

Algorithm 1 Newton's method for Flow

1: " Nf = 1 ⇥ 10 6 2: " qss = 1 ⇥ 10 24 3: k fl max = 9 4: k = 0 5: while kF (x k )k 1 " Nf kF (x 0 )k 1 ^ kF (x k )k 1 " qss ^ k  k fl max do 6: find x k+1 : J (x k ) x k+1 = F (x k ) 7:
x k+1

x k + x k+1 8:

find c k+1 g : T g (c k+1 g ; x k+1 ) = 0 9:

update EOS and physical parameters 10:

k k + 1 11: end while
The user defined or default tolerance " Nf 2 [10 8 , 10 6 ] and maximum number of Newton iteration k fl max for the flow coupling are firstly initialized. Next, the linearized system is solved for the increment x, line 6. All the partial derivatives involved in the Jacobian J of discretized flow system are analytical. For ex. the derivative of gas density interface coe cient is expressed by

@(⇢ g ) ij @p g,i = @(⇢ g ) ij @⇢ g,i @⇢ g,i @p g,i , @⇢ g,i @p g,i = Mi v 2 i (@p g /@v) i , ( 20
)
where @p g /@v is the analytical derivative arising from the corresponding EOS.

One can note that the density derivatives are composition dependent and proportional to the average of molecular weight M as well as the density function Eq. ( 14). Various solvers exist for solving multiphase system of linear equations [14], we apply GMRES [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], that is one of the most prevalent and e cient methods, and ILU(0) [START_REF] Van Der Vorst | Guidelines for the usage of incomplete decompositions in solving sets of linear equations as they occur in practical problems[END_REF] as a preconditioner.

With the velocities and saturations given by x k+1 from step 7, the linear transport system T g (c k+1 g ) = 0 is solved for c k+1 g at line 8 by using GMRES and ILU(0). Due to the modified composition, the EOS parameters must be updated in order to evaluate a new molar volume v by solving the EOS analytically. Then, the physical properties can be calculated, line 9. Thereafter three stop criteria must be checked, line 5: two for the flow system residual and one for the number of iterations.

The proposed coupling 1 corresponds to Newton's family for the flow system in regard to x k+1 , while it can be seen as Picard's method (fixed point method) for the transport equations on c k+1 g . Inasmuch as the gas phase is supposed to be compressible and the significant di↵erence in gas density may occur over all modeled domain, one can include an additional criterion for gas quantity n g deviation:

max n k+1 g n k g n k+1 g  " g , (21)
where max is maximum over the modeled domain. After numerous tests, we deduce, that this can be not necessary condition but su cient to finish the loop that depends on the complexity of gas dynamics. In spite of neglecting the criterion (21), there is no lack of solution accuracy. In addition, the reactive transport coupling 3.3, that goes after the flow coupling 1, entails the convergence (stop) conditions for gas and solid phase and guarantee the conservation of the whole system. An adaptive time-stepping is implemented in regard to the relaxed CFL condition, number of Newton iterations and maximum saturation/pressure changes. Extremely large time steps are not reasonable, even if the scheme is unconditionally stable. Moreover, the following reactive transport can require smaller time step. In this case, the smaller inner time stepping is usually used. Nevertheless, when the RT coupling is finished, the agressive chemistry would provide the important changes returning the larger reaction terms and the stronger modified fluid and rock properties to the flow block at the next time step. So, the flow solution would be further from the first guess and consequently, the convergence rate would decrease as it is quadratic only near the root. This subject will be discussed in Sec. 3.4.1. Notice that there is no calculation of geochemical reactions in this coupling, but the reaction term R ↵ , that can be simply expressed by:

R n ↵ = (V ↵ ⇢ ↵ ) n,rt (V ↵ ⇢ ↵ ) n,f l t n+1 , ( 22
)
where t n is time step n, (•) n,f r and (•) n,rt denote the values obtained by the flow and reactive transport couplings, respectively. So, it is then estimated a posteriori at each time step and remains the mass conservation throughout the calculations. It is also worth mentioning, that when three phases, g, l, s, are in the system, the term R g reflects not only the gas and liquid connection but also the implicit gas-liquid-solid mass transfer, which is not negligible in conventional geochemical modeling. Precipitation and dissolution of some minerals can lead to gas depletion or formation, that impacts on saturation and pressure. These phenomena were numerically demonstrated in [29]. Consequently, the value of reaction term can be of great importance, depending on the physical and chemical problem statement and on the aimed accuracy. E.g. in Vostrikov's work, the reaction terms between the flow for dominant components and the reactive transport were neglected, given the next assumptions: "only very small amounts of minerals are transferred to the liquid form" and "minor components do not have a significant impact on the physical parameters of system". Algorithm 1 constitutes simple solving the two-phase system, where the gas concentrations are assumed to be the secondary variables just like the relative permeability, the capillary pressure and others. In this coupling, the flow can be fully unreactive supposing that the gas displacement is the dominant mechanism and the phase equilibrium calculations are taken explicitly from the reactive transport coupling, in this case.

Pursuing the idea of solving the phase flow formulation (5) and not the compositional problem, other options of coupling were tested, for example:

Step 1. the phase conservation system to trace the fist estimation of velocities and saturations; Step 2. the reactive transport coupling and finally Step 3. the phase flow problem, solving Newton's method for the second time, but taking into account the phase transfer due to the geochemical reactions. This method demonstrated disadvantageous e↵ects on the convergence and performance because of the delayed treatment of the fluid properties at Step 1. Thus, the gas properties need to be implicitly estimated in the iterative loop of flow to guarantee the convergence and unconditional stability, especially if in the aimed application, the mechanical displacement of gas front occurs, so the advective dominant flow takes place. As a result, the flow system is fully implicit and inherits its advantages in contrast to the time-step restriction of IMPES method.

Coupling 2: Reactive transport

Relying on Algorithm 1 allows to find the phase velocities and saturations and consequently, to manage the reactive transport problem which consists of the linear transport and nonlinear chemical equations. So, analogous to the standard coupling in saturated porous media, Picard's method can be employed. On the one hand, the gas and liquid transports can be summed and therefore written for each species similar to GIA for multiphase multicomponent flow [15,[START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF][START_REF] Lichtner | Reactive transport in porous media[END_REF] by applying phase appearance and disappearance criteria and K-values K i = y i /x i , that would be K i = K h i i /(P ' g i ) in this work, derived from ( 13). In the non-ideal approach, i is calculated regarding electrostatic state of ions through the ionic strength, furthermore, ' g i is the function of the whole mixture properties including y i itself, then K i depends both on aqueous and on gaseous composition and should be updated at each iteration as secondary variable. On the other hand, remaining liquid (7) and gas ( 8) operators apart, one can apply the gas concentration from Algorithm 1 as a first guess and avoid K-value calculation by introducing the reaction terms. By denoting transport operators by T ↵ (c ↵ ) = 0 and geochemical reactions, described in Sec. 2.3 by R(c), where c = {c l , c g , c s }, we propose the SI reactive transport algorithm 2.

Algorithm 2 Picard's method for Reactive Transport

1: " rt = 1 ⇥ 10 5 2: k rt max = 60 3: k = 1 4: while ✓ max |c 2k+2 g c 2k g | c 2k+2 g " rt ◆ _ ✓ max |c 2k+2 s c 2k s | c 2k+2 s " rt ◆ ^(k  k rt max ) do 5: k k + 1 6: T g (c 2k+1 g ) = 0 7: T l (c 2k+1 l ) = 0 8: c 2k+2 = R(c 2k+1 ) 9: end while
Let us set out the specific aspects of the SI reactive transport algorithm 2. The reaction term, devoted to replicating dissolution and evaporation rate, is introduced in gas transport operator at line 6. Omitting the porosity and the external sources, it can be formulated for gas species i as

S n+1 g c n+1,2k+1 g,i S n g c n g,i t n+1 = (T g,i ) h + S n+1 g t n+1 r n+1,2k+1 gl,i , (23)
r n+1,2k+1 gl,i = d n+1,2k i c n+1,2k+1 g,i + e n+1,2k i , (24)
where (•) h denotes the discretization, d n+1,2k

i and e n+1,2k i reflect the dissolution and evaporation processes, respectively. The terms d n+1,2k and e n+1,2k are calculated at the previous iteration 2k, while r n+1,2k+1 gl corresponds to the iteration 2k + 1. Considering the concentrations obtained on the previous iterations {1, . . . , 2k}, these parameters are defined as

d n+1,2k i = k X m=1 R >0 ( c n+1,2m g,i ) c n+1,2m g,i c n+1,2m g,i , (25)
e n+1,2k i = k X m=1 R 0 ( c n+1,2m g,i ) c n+1,2m g,i , (26)
where c n+1,2m

g,i = c n+1,2m g,i c n+1,2m 1 g,i (27)
and R >0 (•) is the indicator function of the set of strictly positive real numbers. Note that equation ( 28) is still linear and can be directly solved by one of the linear solvers, GMRES with ILU(0) for example. This new distribution of gas species i entails di↵erent mass transfer that should be associated with the corresponding total liquid mobile concentration of primary species j. Then, the liquid transport operator, line 7, is defined by

S n+1 l c n+1,2k+1 l,j S n l c n l,j t n+1 = (T l,j ) h + S n+1 l t n+1 r n+1,2k+1
lg,j + R ls,j (c s,j ), ( 28)

r n+1,2k+1 lg,j = X i ↵ ij r n+1,2k+1 gl,i , (29)
where R ls,j is the reaction term of liquid-solid interaction whose detailed description and variable porosity management can be found in [START_REF] Lagneau | Operator-splitting-based reactive transport models in strong feedback of porosity change: The contribution of analytical solutions for accuracy validation and estimator improvement[END_REF]. Following the transport, the reactive part aims at finding a new local equilibrium state, line 8. The nonlinear chemical system for basis species is solved at each cell by Newton's method with the line search procedure. The Jacobian is analytically evaluated. To improve the convergence, a new basis can be chosen independently in each cell, hence its set of primary equations also changes. Then, the gas and solid concentrations changes are subjected to stop criterion verification, line 4, if the phase interactions, liquid-gas, liquid-solid, take place.

Coupling between compressible two-phase flow and reactive transport 3.4.1. Type of coupling

The sequential coupling of the flow 1 and reactive transport 2 modules has a specific limitation arising from the former, which is supposed unreactive. The flow method was initially devised for modeling the gas appearance due to the mechanical displacement but not due to the phase transition, from liquid to two-phase state. Nevertheless, within the reactive transport module, the evaporation, dissolution and other chemical processes are permitted and consequently, it changes the occupied volume and mass of each phase. Also, the gas phase can vanish or appear. The impact of such modifications varies with the chemical system of problem. We manage the changes arising from the RT in the flow block a posteriori, by updating the fluid and rock properties and the explicit reaction terms R ↵ , Eq. 22. This approach is valid if chemical reactions rate and time steps remain small; otherwise, the severe changes in the fluid and rock properties hinder convergence of Newton's method, as initial guess is far from the solution. Moreover, the explicit treatment of geochemical impact on the flow can be significantly underestimated. To overcome this problem and to model the reaction driven advection, the tight coupling, similar to that of STOMP [START_REF] White | STOMP Subsurface Transport Over Multiple Phases, Version 4.0, User's Guide[END_REF], should be applied in the case of the high change rate. It permits to update the velocities and the fluid/rock properties implicitly: iterative procedure of flow and SI reactive transport. However, when the mechanical force is dominant, it is computationally e cient to use the sequential coupling between the flow and the RT module.

Applying GIA avoids this issue by definition, since the equilibrium state is calculated at each Newton iteration by using the thermodynamics flash (Rachford-Rice equation) [START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF] or the complementary conditions. Instead, the compositional flow system consists of a larger set of nonlinear equations than that of the proposed method. Furthermore, in case of variable switching, the primary variables should be adapted to the local equilibrium state, and it requires to rebuild the Jacobian matrix.

Numerical assessment

By means of this approach the whole problem is divided into two subsystems: the compressible flow coupling 1 composed by N f nonlinear equations of flow and N g linear equations of gas transport, then the reactive transport coupling 2 of N c + N g linear transport equations and N c + N kin nonlinear equations from chemical system. Let us estimate the computational impact of each part. Supposing the uniform time step both for flow and for reactive transport, we denote the calculation time of the whole system t, the flow coupling per iteration t flc , the flow operator t fl , the gas transport t gt and for the reactive transport per iteration t rtc , then

t = N Nit t flc + N P it t rtc = N Nit (t fl + t gt ) + N P it t rtc , (30)
where N Nit is the number of Newton iterations, N P it is the number of Picard iterations. By taking into account t gt ⌧ t fl and N P it t rtc 2 (t fl /10 3 , N Nit t fl ], we deduce that the calculation time per time step for this work can be expressed by following

N Nit t fl < N Nit t fl + N P it t rtc  2N Nit t fl , (31)
where the lower border corresponds to the problem with a low geochemical complexity (e.g. CO 2 and H 2 O) and the calculation time t is mostly based on the flow operator part, making t fl decisive factor, as a result. It can be reduced by performant linear solver and preconditioner [37] or/and by using the AIM for the flow system. The upper limit can be reached, when the agressive geochemistry is modeled with the dis-/appearance of solids, gases, the redox reactions etc. Although, the calculation time t is still strongly dependent on t fl . Thus, we conclude that the more the chemical system is abundant in components, the more the proposed flow coupling is advantageous in comparison with the methods based on the component flow formulation, considering that the larger time stepping in GIA can be compromised by the required reactive time step.

The inner couplings have their own tolerance, another benefit of the OSA: [10 8 , 10 6 ] for flow, 10 5 for reactive transport and 10 12 for the chemistry, that results in the mass balance error in the range [10 6 , 10 5 ] according to the test cases carried out for di↵erent types of geochemical and hydrodynamic complexity.

Numerical simulations

The proposed method is first tested by modeling a problem with the selfsimilar solution. Next, we apply it to the simulation of CO 2 injection that permits to evaluate the capabilities of the simulator to represent the physical behavior and its computational e ciency.

1D axisymmetric problem: radial flow from a CO 2 injection well

The axisymmetric problem of constant injection in a saturated, horizontal, infinite reservoir admits the self-similar variable R 2 /t, whose property makes it an ideal way to validate the numerical code. In addition, this advantage was availed to model a constant 100 kg/s CO 2 injection in a long aquifer 100000 ⇥ 100 m 2 , in the workshop "Intercomparison of numerical simulation codes for geologic disposal of CO 2 " initiated by Lawrence Berkeley National Laboratory [START_REF] Pruess | Intercomparison of numerical simulation codes for geologic disposal of co2[END_REF]. We use the parameters presented in [START_REF] Pruess | Intercomparison of numerical simulation codes for geologic disposal of co2[END_REF], adapting some of them: as the pore compressibility is neglected, we intensify the intrinsic permeability, K = 2 ⇥ 10 13 [m 2 ]; the water disappearance is prohibited; the fluid properties and the solubility are treated di↵erently. In HYTEC, the gas and liquid densities are provided by the Peng-Robinson models [2,36,[START_REF] Robinson | The characterization of the heptanes and heavier fractions for the gpa peng-robinson programs[END_REF], the viscosities are predicted by [5,34]. The partial molar volume at the in-Figure 1: 1D axisymmetric problem: gas saturation S g . Results of HYTEC (in color) and of TOUGH2-ECO2, [START_REF] Pruess | Intercomparison of numerical simulation codes for geologic disposal of co2[END_REF] (black) finite dilution in the water is averaged over the relevant pressure range; the aqueous activities and gaseous fugacities are simulated according to the bdot and the PR78 [START_REF] Robinson | The characterization of the heptanes and heavier fractions for the gpa peng-robinson programs[END_REF] models, respectively. In [START_REF] Pruess | Intercomparison of numerical simulation codes for geologic disposal of co2[END_REF], TOUGH2-ECO2 module was used, whose description is given in [START_REF] Pruess | Multiphase flow dynamics during co2 injection into saline aquifers[END_REF].

Despite the mentioned deviation in parameters, the gas saturation front is similar in both cases, especially over the fist 1000 d, Fig. 1. The increased discrepancy in saturation front position at 10000 d arises from the water disappearance zone modeled in [START_REF] Pruess | Intercomparison of numerical simulation codes for geologic disposal of co2[END_REF]. Notice that the the saturation curves of HYTEC are almost perpendicular to the R-axis in contrast to those of TOUGH2-ECO2, their slope grows with time that is just the truncation error e↵ect. The results of gas saturation as a function of R 2 /t are illustrated in Fig. 2 and demonstrate a high accuracy of the numerical code. The next 2D problem of CO 2 injection in a fully water-saturated domain was proposed by [START_REF] Neumann | Modeling and simulation of two-phase two-component flow with disappearing nonwetting phase[END_REF]. The injection of constant rate 0.04 kg/m 2 /s is modeled through the left bottom boundary (whose length was not provided in [START_REF] Neumann | Modeling and simulation of two-phase two-component flow with disappearing nonwetting phase[END_REF]) of the rectangle reservoir 600 ⇥ 100 m 2 . In this work, the injection border is 10 m long, then the debit is 0.4 kg/s. The Dirichlet boundary conditions are set at the right border of the reservoir: hydrostatic pressure and S g = 0. By neglecting the dispersion and assuming the Millington-Quirk tortuosity model [START_REF] Millington | Permeability of porous solids[END_REF], the e↵ective di↵usion for this problem takes the following form:

D e ↵ = ! 4/3 S 10/3 ↵ D ↵ . ( 32
)
The geometry and grid dimension are taken after Neumann et al., while the fluid properties and the solubility di↵er from the reference. The modeling of fluid properties and of phase equilibria is analogous to that given in Sec. 4.1.

The fluid dynamics is represented exactly as in the reference [START_REF] Neumann | Modeling and simulation of two-phase two-component flow with disappearing nonwetting phase[END_REF]: the CO 2 forms a bubble that grows and rises upward, then the current is distributed gas saturation from [0, max(S g )] and contours of CO 2 mole fraction x l,CO 2 at 0.005, 0.011, 0.016 along the top of the aquifer by gradually developing its area, Fig. 3. The gas saturation is lower than that of the reference performed by DUNE, the mole fraction of CO 2 x l,CO 2 is higher. This might be due to an uncertainty on the injection rate used in [START_REF] Neumann | Modeling and simulation of two-phase two-component flow with disappearing nonwetting phase[END_REF] and due to the di↵erent fluid properties and phase equilibria models.

The same grid dimension as in [START_REF] Neumann | Modeling and simulation of two-phase two-component flow with disappearing nonwetting phase[END_REF] is chosen for this simulation, 240 ⇥ 40 cells. The maximum number of Newton iterations is set to 9 in HYTEC for this problem. The tolerance of Newton's method (Algorithm 1) is " Nf = 1 ⇥ 10 7 and the gas quantity criterion Eq. ( 21) is " g = 1 ⇥ 10 5 . The mass conservation error is of order 10 6 .

The initial time step is set to 156 s, which is also the minimum time step min(dt). The time stepping for HYTEC and DUNE is presented in Table 3. The HYTEC's imposed maximum time step max(dt) is taken slightly higher (by 8%) that results in a larger average dt and in a faster execution of HYTEC, by at least 9%.

In Table 4, the execution time of HYTEC includes the grid construction, initialization, output printing, solvers, secondary property modules. During our simulation, the total number of time steps (successful and unsuccessful) is 39% less than that of DUNE, that can be explained by slightly higher average time step and rate of Newton's method convergence in HYTEC.

By employing our method, the unsuccessful steps appear right after the moment when the gas current reaches the right boundary, hence it is restricted to the further CO 2 (g) propagation. Nevertheless, the dissolved CO 2 can be released; the depressurization thus increases the physical complexity of the system, hence the unsuccessful Newton iterations.

Conclusion

A new solution method of compressible multiphase flow was proposed, that can be integrated as a module in conventional reactive transport frameworks, based on the operator splitting or global implicit approach. The flow method is composed by the phase conservation formulation, the gas transport and the equations of state. The versatile structure allows to conserve the constant number of nonlinear equations of the flow problem, independent of the modeled geochemical system, while chemistry basis can change during one time step to enhances the convergence. This characteristic makes it advantageous for modeling the multicomponent problems. The flow phase system can be also presented by other formulations of primary variables, e.g. pressure/capillary pressure. It is worth mentioning that it is possible to generalize the variable switching in rigorous way. The whole flow coupling preserves the fully implicit advantages of the multiphase flow discretization. The present method was implemented in the SIA based reactive transport simulator HYTEC. The numerical code was verified by modeling the benchmark problem admitting the self-similar solution, the computational e ciency was showed by the simulation of CO 2 injection and compared with that of DUNE.

Verification

Once the numerical code is built, it should pass the prerequisite verification steps:

• mass balance conservation,

• modeling an analytical solution,

• benchmarking, reference models.

In Sec. 3.2, we have demonstrated the code capabilities to represent the semi-analytical and reference models. Here, we give additional information on the numerical results presented above, Sec. 3.3.1 and 3.3.2. Then we propose a benchmark study, Sec. 3.3.3. The mass conservation was calculated for all modeled problems in this work and the order of relative error is less then 10 -6 .

3.3.1 1D axisymmetric problem: radial flow from a CO 2 injection well, Pruess et al. (2002) The parameters of the model can be found in Table 2.1. The fluid densities are given by the PR78 and T-VTPR models (see Eq. (4.1.19) and Sec. 4.2) for which the binary interaction parameter are calculated by the PPR78. The gas viscosity is predicted by Altunin75, Eq. (4.3.8); we use the Islam-Carlson model of the liquid viscosity, Eq. (4.3.12). Modeling the 1D radial problem with self-similar solution by the full coupling reveals the growing impact of solubility and compressibility, compare with Figs. 2.4(a), 2.5 and 2.6. 3.1 -Parameters of the 2D problem of CO 2 injection [START_REF] Neumann | Modeling and simulation of two-phase two-component flow with disappearing nonwetting phase[END_REF].

2D problem: CO

Initial Conditions Ω = [0, 600] × [0, 100] [m 2 ] Ω : p l0 = p atm + ρ l g (900 -z) [P a], S g0 = 0 Temperature 40 • C Boundary Conditions Γ N = [(0, 0); (0, 10)] [m] Γ D = [(600, 0); (600, 100)] [m] Γ N : ψ l = 0, ψ g = 4 × 10 -2 [kg/m 2 /s] Γ D : p l0 = p atm + ρ l g (900 -z) [P a], S g = 0 ∂Ω\(Γ D ∪ Γ N ): ψ l = 0, ψ g = 0 Matrix Properties Porosity φ 0.2 Permeability K 10 -12 [m 2 ]
Fluid-Matrix Properties Brooks-Corey capillary pressure and relative permeability Eq. 2.1.7, (2.1.9) and ( 2

Benchmark proposition

In the previous sections, we implemented our coupling method for MMF and RT (3.2).

Some validation work was performed either based on analytical solutions or literature results. A full description of the thermodynamic laws was then presented, including some applications to the CH 4 -CO 2 -H 2 S-H 2 O system (4.2).

Based on these developments, a new application is proposed. The application is built from the work by Bonnaud (2012); Bonnaud et al. (2012): the authors attempt to simulate a mechanism possibly responsible for the creation of heterogeneities in the composition of gas reservoirs. In their scenario, an initially sour reservoir is leached by an underlying aquifer; the differential solubility of the gas components leads to local modifications of the gas composition. The system includes a tight connection between transport and geochemistry. The evolution of the mass/volume due to the reactive transport affects the pressure field, which leads to a strong feedback of chemistry on flow. An accurate modeling of a multicomponent gas phase is also required for this problem. Finally, the complexity of the CO 2 and H 2 S geochemistry in the host rock is of particular interest.

Result analysis based on physical meaning also reach its limits due to the numerous couplings and feedbacks. We therefore decided to use this problem as the basis for a numerical benchmarking exercise, open for the reactive transport modeling community.

The full description of the context and the increasingly complex exercises proposed in the benchmark are detailed in the following paper Sin et al. (submitted), provided below. The numerical results are partly reported in Sin et al. (2015).

Introduction

Reactive transport codes are quantitative tools that allow to explore the behavior of complex geoscience systems: coupling and feedbacks between flow and transport of fluids in the pore space and reactions between aqueous species, mineral surfaces, rocks and a possible gas phase. They rely on the resolution of several sets of equations describing the hydrodynamics and the geochemistry of the system [1]. Reactive transport simulation are now used in a wide variety of applications, in numerous domains and at different scales. E.g., the HYTEC code, developed at MINES ParisTech [2,3], has been applied to cement degradation [4], radioactive waste disposal (performance assessment, near-and far-field processes [5,6]), geological storage of acid gases [7,8,9], uranium in situ recovery processes [10].

For the gas-water-rock interaction (e.g. for the simulation of enhanced gas or oil recovery, gas storage, CO 2 sequestration), the models are based on coupled, non-linear equations for multiphase multicomponent reactive flow and transport coupled with chemistry and water-gas equilibrium. Different ways of coupling are used to handle the mathematical and numerical difficulties [11,12,13,3,14,15]. The significant mathematical and computational challenges impose constant verification of the accuracy of the codes [16]. Unfortunately, due to the non-linearity and numerous couplings, analytical solutions are limited to very specific and overly simplified or de-coupled test cases [e.g 2,17]. Alternatively, validation on experimental data (from laboratory or field)

focuses on the accuracy of the physical and chemical laws included in the codes, rather than the accuracy of the numerical schemes. Inter-comparison between codes is therefore a key for code verification.

Several such benchmarking exercises have been promoted in the recent literature [18,16,19,2,20]. Also noteworthy is the subsurface environmental simulations benchmarking initiative [21,22], which aims at promoting an international community around reactive transport benchmark exercises definition and resolution.

This paper presents such an reactive transport benchmarking exercise. The aim is to specifically target the numerical capacity of the codes in the context of multiphase (water and gas) flow, with several coupling between water and gas transport, multicomponent gas solubility, and gas-water-rock interaction.

The exercise stands half-way between pure numerics and direct application: the overall description should be geologically realistic, the multiple couplings should be preserved so that the verification is demonstrative, and finally the complexity should be moderate in order to keep the computational intensity as reasonable as possible.

The exercise is based on the description of the long term evolution of a sour gas reservoir. The geometry was simplified and several simplifying hypotheses are taken to allow for different complexity applications. A general overview of the geological context is proposed in the following section: the focus is on the physical and geochemical processes at stake in the system. Then, the usually-accepted equations describing these processes are presented, with a view to identifying the mathematical bases of the system and highlighting the couplings. A description of the exercises is then given: the actual geometry, parameters, initial and boundary conditions are detailed. Finally, some insights to the results (using the HYTEC code) are given: however, the paper does not provide a comprehensive description of the results to allow some liberty for potential benchmarkers.

Geological context of the exercise

Compositional heterogeneities of hydrogen sulfide (H 2 S) have been noticed in many sour gas reservoirs. H 2 S occurrence is an important factor of economic depreciation. Thus, the knowledge of its content and distribution is a critical parameter when planning field development.

A possible scenario for H 2 S segregation in gas reservoir stems from the preferential leaching of H 2 S by underlying aquifers [23,24]. Under conditions of pressure and temperature of typical reservoirs, H 2 S and CO 2 are far more soluble than hydrocarbons, Fig. 1. A preferential leaching of H 2 S (and to a lesser extent CO 2 ) versus CH 4 over time takes place. In the case of a deep and trapped reservoir of acid gases in contact with a bottom aquifer, Fig. 1, this differentiation is controlled by several mechanisms [23].

1. The differential solubility of gases changes the relative amounts of each gas component near the gas-water contact.

2. An active aquifer can effectively export the dissolved gases, thus enhancing dissolution on the long-term. 3. Higher ratio between the length of the gas/aquifer interface and the height of the reservoir maximizes the dissolution, leading to larger amounts leached.

4. Diffusive and advective transport in the gas phase transfers the compositional anomalies farther from the gas-water contact.

This exercise is a simplified simulation of the behavior of such reservoir: the gas phase is mostly composed of methane (CH 4 ), with significant contents in H 2 S and carbon dioxide (CO 2 ). In order to accurately represent the system, the simulation of several phenomena is required. A two-phase (gas, water) system description is obviously a pre-requisite. The description of a multicomponent gas phase is also compulsory to follow the evolution of the composition in the reservoir. The description of the aquifer implies the resolution of the water flow. The dissolution of the gas in the aquifer is bound to modify the local (gas) pressure, so that gaseous flow should also be represented. From a physicochemical point of view, accurate water-gas interactions can be described by several equations of state: ideal gas is of course the simplest, but high pressure high temperature conditions require the use of more complex equations. Finally, the description should include some water-rock interaction, involving the dissolved acids (CO 2 (aq), H 2 S(aq)) and host rock minerals.

These phenomena are strongly coupled. Beyond the (two-phase) flow, transport and chemistry coupling, several interesting feedbacks should be considered.

The evolution of the gas composition modifies the local gas density (H 2 S is more dense than CH 4 ) which impacts the (compressible) gas flow. The physicochemical dissolution of gas components in the water phase removes some mass from the gas phase; this should be accounted for by resetting the pressure-volume in the flow equations, using mass-volume relationships heavily dependent on gas composition and pressure. Assuming incompressible flow and no mass exchange in the flow system for the majority of the benchmark cases allowed to simplify the system for benchmark participants and to focus on the impact of different thermo-and hydrodynamical factors [25].

Several mathematical frameworks exist to represent each of the processes described above. The most generally used governing equations and numerical methods for these problems are given in the following section.

Mathematical Model

Two-Phase Flow

We consider an isothermal two-phase flow with compressible gas phase and reactive multicomponent transport in porous medium. Let us introduce a twophase flow model, which states a system of partial differential equations of mass conservation for each phase α, liquid and gas,

∂ωρ α S α ∂t + ∇ • (ρ α u α ) = q α + R α (1)
where the fluxes can be expressed for each phase by the Darcy-Muskat law [26,27,28]

u α = - k rα µ α K (∇p α -ρ α g) , (2)
where K = KI and K is a scalar coefficient of intrinsic permeability under the isotropic conditions. The source term q α can be also presented as

q α = ρ α φ α . (3)
To close the system we get the following constitutive relations:

S l + S g = 1, (4)
p g -p l = p c . (5)
The capillary pressure p c and relative permeability k r are presented by empirical laws, e.g. Brooks-Corey model [29]:

S α = S α -S αr 1 -S lr -S gr , (6)
p c (S l ) = p e S -1 λ l , (7)
k rw (S l ) = S 2+3λ λ l , (8)
k rn (S l ) = (1 -S l ) 2 (1 -S 2+λ λ l
). ( 9)

The gas density is calculated using the (cubic) equations of state that will be discussed next, see section 3.3.

Numerous formalisms are proposed in the literature to solve the flow equation [30,31,11,13,32]. They can rely on substitution of several unknowns [32], so that the system is described by pressure/saturation and pressure/mole fraction in liquid phase, or the formalism of pressure/pressure [31] or pressure/mass density [11] can be applied. In the HYTEC code, the resolution is built on a water pressure/gas saturation formalism, on a Voronoi-based finite volume discretization and on a fully implicit method with analytical Jacobian. An adaptive time-stepping is controlled by the number of Newton iterations. We neglect the liquid compressibility and viscosity variation.

Multicomponent Transport

Based on the principal components formalism (see 3.4), we define a limited number of basis species N tr p and denote their total liquid mobile, immobile and gas concentration by c l i , ci and c g i . The liquid and gas transports are separated, the former combines liquid mobile and immobile part of species i:

∂ωS l c l i ∂t = T c l i + q l i + R l i c l i , cl i , ( 10
)
∂ωS g c g i ∂t = T (c g i ) + q g i + R g i (c g i ) ( 11
)
where operator T includes the advective flux presented by Darcy's law and the diffusive-dispersive flux [26,33,12] T (

c α i ) = ∇ • (D ef f α ∇c α i -c α i u α ). ( 12
)
Assuming the equality of longitudinal and transverse dispersion coefficients, the effective diffusion can be written:

D ef f α = ωS α D α + a α u . ( 13
)
The Millington-Quirk tortuosity model [34] is often used, in this case the effective diffusion takes the following form:

D ef f α = ω 4/3 S 10/3 α D α + a α .u (14)
In HYTEC, the transport resolution is performed by an upwind flux discretization with the harmonic weighting of diffusion. Applying a semi-implicit method for time discretization makes it possible to achieve an accurate solution:

125 implicit Euler in time for diffusive-dispersive and Crank-Nicholson for advective flux. Other authors prefer to unify both phases and resolve implicitly fully coupled mass or mole balance for each component [35,13,32].

Equation Of State

An equation of state (EOS) is required to define the volumetric properties of the system: it establishes a relationship between pressure, molar volume and temperature (P V T ) and allows to evaluate density, viscosity, fugacity, solubility.

For any pure homogeneous fluid in the equilibrium state and for the region of single phase, this relation can be expressed by the following equation [36]:

f (P, V, T ) = 0. ( 15
)
The ideal gas is the simplest form of equation of state, although it is limited to low pressure/temperature regions:

P V = RT (16)
More complex EOS were developed for higher temperatures and pressures, notably the cubic equations of state which offer an accurate description of both liquid and gas phases over a wide pressure/temperature range. E.g. the Peng-

Robinson EOS [37]:

P = RT V -b - a(T ) V (V + b) + b(V -b) (17)
which uses the binary interaction parameters and can handle multicomponent gases. So, applying one of the EOS provides mass density for gas mixture:

ρ g = Np i=1 y i M i V . ( 18
)

Chemistry 130

Most reactive transport codes rely on the primary species formulation [15,12,38], where all species S j , j = 1, .., N s can be expressed as a function of primary species C i from established basis of N p species, where N p < N s :

S j Np i=1 α ij C i (19)
Their concentrations are constrained by the Mass-Action Law: 20) using stoichiometric coefficients α ij and activity coefficients γ j . For each basis species, the total concentration T i , i = 1, .., N p can be written through

S j = K j γ j Np i=1 (γ i C i ) αij (
T i = Ns j=1 α ji S j , (21)
T i - Ns j=1 α ji K j γ j Np i=1 (γ i C i ) αij = 0, ( 22
)
the latter constitutes a system of mass conservation on C i .

The kinetics of mineral precipitation/dissolution can be treated by a transition state theory type rate law [39]: where A s is the specific surface area, k is the dissolution/precipitation rate constant and a i is the activity of the potential catalyzing or inhibiting species.

dS j dt = A s k Q K s -1 a ni i ( 23
)
Q s is the ion activity product defined by

Q s = γ j S j Np i=1 (γ i C i ) αij . ( 24
)
Mass-Action Law for all derived species at equilibrium, kinetic low for kinetically limited minerals, and mass conservation achieve a complete set of algebrodifferential equations. Its resolution, at each node of the domain, yields the concentrations of all the species in the system.

Benchmark Description

The benchmark exercise is based on a simplified description of the geological example. The geometry is shown in Fig. 2: the 2D model is similar simplification (or a crop-out) of the general geometry Fig. 1. There are two main zones, reservoir and aquifer; the third acts as a contact/transition zone. The simulation time is 10 5 y.

φ l = (4.8e-9, 0) [m/s] ( 25
)
Γ D : p l = 1e5 + ρ l g(3.3e3 -z) [P a] (26)
The geochemical and hydrodynamic description can be found in Tables 3 and4. The initial condition gas is calculated for temperature 100 Brooks-Corey, (p e = 1e5, λ = 2, S lr = 0.01), Eq. ( 7)

k rα Brooks-Corey, (λ = 2, S lr = 0.01), Eq. ( 9)

p l , [P a] 1e5 + ρ l g(3.3e3 -z) T , • C 100 isothermal α l g D α , [m 2 /s] 8e-9 5e-8 a α , [m] 1 - µ α , [P a • s] 1e-3 3.9e-05
ρ α , [kg/m 3 ] 997 215.5 Several cases were set up with increasing complexity:

1. Geochemistry: the gas composition is set for all cases, but calcite, then calcite+siderite are added for the more complex cases;

2. EOS: ideal gas or Peng-Robinson; Table 5 summarizes the scenarios: GCase 1 to 3, increasingly complex chemistry, ideal gas and simplified flow, FCase 1 to 7 with simple chemistry, cubic EOS and increasing flow complexity. For sake of simplicity, the assumption of incompressible flow excluding mass exchange in flow module is accepted in all cases except for FCase 7, that comprehends a tight coupling with the chemistry feedback of gas-water-rock interaction.

Numerical Results

We present some simulation results obtained using HYTEC, for scenario FCase 1 (Table 5): advection in the whole domain, ideal gas and incompressible flow assumptions. The focus in this paper is more on the physical description of the simulation results: the numerical and accuracy aspects will come in a second time when comparing the results from different codes.

During first 2e4 y, the system is dominated by the gravity driven gas flow:

the gas saturation is redistributed in the reservoir, with a capillary fringe at the bottom, and gas concentrating at the top of the formation with saturation as high as 0.98 (Fig. 4). The system approaches to a stationary state, so the capillary fringe motion decelerates afterwards. Enhanced H 2 S dissolution in the aquifer modifies the H 2 S content, with a sharp decrease at the contact with the aquifer (Fig. 3). H 2 S dissolution is the strongest upwards in the aquifer, due to fresh water renewing by the aquifer, creating a visible horizontal H 2 S(g) gradient (Fig. 3, 70 ky), until nearly all the H 2 S is evacuated (100 ky).

The effect of the EOS can be observed Fig. 5: the figure illustrates results

from FCase 1 and FCase 3. Both scenarios are identical (no mineral, two-phase flow, incompressible gas), and differ only by the EOS: ideal gas for FCase 1,

Peng-Robinson for FCase 3. The preferential leaching of H 2 S is observed in both models, but the mole fraction for FCase 3 stays significantly higher than that for FCase 1. This is highlighted by the evolution of the total mass in the system, Fig.

Calcite

Calcite

K R = 0 K R = 0 K R = 0 Diff.
Diff. initial content and higher solubility), and starts decreasing after 2×10 4 y, due to decreasing accessibility. More remarkable is the large difference of behavior due to the choice of EOS: ideal gas assumption largely overestimates H 2 S solubility.

CF

The effect is minor on CO 2 and CH 4 , although they are not in the same direction (slight overestimation for CH 4 , slight underestimation for CO 2 ).

Conclusion

The multiphase multicomponent benchmark is proposed to the reactive transport modeling community. The exercise is based on geologically realistic foundations, and allows to explore a complex coupling involving gas-water-rock interactions. The system includes a tight connection between transport and geochemistry. The evolution of the mass/volume due to the reactive transport affects the pressure field, which leads to a strong feedback of chemistry on flow.

An accurate modeling of a multicomponent gas phase is also required for this problem. Finally, the complexity of the CO 2 and H 2 S geochemistry in the host rock is of particular interest.

The benchmark problem was formulated and two sets of cases were constructed to investigate the geochemical, themo-and hydrodynamical factors.

The sensibility study showed a strong quantitative impact of the choice of equation of state, diffusion model and advection in reservoir on gas leaching.

The choice of observables for comparison between codes is not finalized yet.

It will involve a combination of color maps, profiles at several rendezvous times, local evolutions at prescribed positions for gas saturation, pressure and selected species. The evolution of the total mass is also an interesting (although vastly integrative) observable. A more detailed description can be found in [25].

Overview of equations of state

Let us give a brief survey of some existing EOS models and their applicability to different type of systems. The ideal gas law Eq. (3.1.6) is the simplest form of equation of state, although it is limited to low pressure/temperature regions. More complex EOS were developed for higher temperatures and pressures, notably the cubic equations of state which can offer an accurate description of both liquid and gas phases over a wide pressure/temperature range. Since the Peng-Robinson EOS provides a better estimation of compressibility factor, vapor pressure and liquid density than SRK, we take it as a base model for further modeling in this study.

Peng-Robinson model

The Peng-Robinson model (PR78) [START_REF] Robinson | The characterization of the heptanes and heavier fractions for the gpa peng-robinson programs[END_REF] for a pure component i takes the following form:

P = RT v -b i - a i (T ) v(v + b i ) + b i (v -b i ) , (4.1.1)
where the repulsive and attractive terms are expressed using the critical pressure, temperature and acentric factor, P c,i , T c,i , ω,

a i (T ) = 0.253076587 R 2 T 2 c,i P c,i α i (T ), (4.1.2)
α i (T ) = (1 + m i (1 -T r,i ) 2 , T r,i = T /T c,i , (4.1.3)
m i =      0.37464 + 1.54226ω i -0.26992ω 2 i if ω i ≤ 0.49, 0.379642 + 1.48503ω i -0.164423ω 2 i + 0.016666ω 3 i else, (4.1.4)
b i = 0.0777960739 RT c,i P c,i . (4.1.5)
The difference between PR76 and PR78 is m i parameter, namely the ω condition. In the mixture, the parameters a and b of PR78 are calculated through the combining rule for the cross parameter a ij and the quadratic mixing rules for a and b:

a = i j x i x j a ij , (4.1.6
) where k ij are the binary interaction parameters. By using Eqs. (4.1.6)-(4.1.8), the cubic EOS yields the phase behavior of mixtures of nonpolar and slightly polar compounds and it gives a less adequate VLE description for the strongly polar systems, such as alcohols and water. A large spectrum of methods and modifications now exists to improve EOS models. We first discuss the impact of modified α and b functions on the VLE prediction.

a ij = √ a i a j (1 -k ij ), (4.1.7) b = i x i b i , (4.1.8)
Then the Volume Translation (VT) method and the derivation of mixing rules from the excess Gibbs energy will be concisely described below. The generalization of two/three constant EOS and also the derivation of main EOS except can be found in [START_REF] Schmidt | A modified van der waals type equation of state[END_REF], [START_REF] Ji | A systematic study of cubic three-parameter equations of state for deriving a structurally optimized pvt relation[END_REF], [START_REF] Michelsen | Thermodynamic Models: Fundamentals & Computational Aspects[END_REF].

α(T r ) = α (0) (T r ) + ω α (1) (T r ) -α (0) (T r ) , (4.1.9)
α i (T r ) = T N i (M i -1) r exp L i 1 -T N i M i r , i = 0, 1, (4.1.10)
where L, M , N are the constants that can be either defined for each compound [START_REF] Twu | A new generalized alpha function for a cubic equation of state part 1. peng-robinson equation[END_REF] or generalized for hydrocarbon fluids by using two points P r = P/P c = 1 and T r = 0.7. Besides, the generalized parameters differ according to the critical temperature point, see Table 4.1. The vapor pressure for light and heavy hydrocarbons was predicted with 3.3% ADD compared to 12.1% of PR76 (using Eq. ( 4.1.4)) from the triple point to the critical one. The L, M , N parameters were later adjusted to a wider range of components Table 4.1 -Parameters for the TPR EOS model: L, M , N parameters of Eq. ( 4.1.10) for hydrocarbon fluids [START_REF] Twu | A new generalized alpha function for a cubic equation of state part 1. peng-robinson equation[END_REF]; (L), (M ), (N ) parameters, generalized for alkanes, aromatics, ketons, alcohols, refrigerants and gases at T r ≤ 1 (Ahlers and Gmehling, 2002b).

T r ≤ 1 T r > 1 α (0) α (1) α (0) α (1)
L 0.125283 (1.511442)

(T ) =      exp [c 1 (1 -T r )] 1 + c 2 1 - √ T r 2 + c 3 1 - √ T r 3 2 if T r ≤ 1, exp [c 1 (1 -T r )] else, (4.1.12)
where c i , i = 1, 2, 3 are the fitted parameters that have the polynomial form of ω:

2 k=0 c i,k ω k , c i,k
are the generalized constants, calculated over the selected compounds through the objective function of pressure. Using c i,k estimated for the system of 22 pure compounds, the vapor pressure was predicted with 1.2% AAD ( 5% AAD for water) in contrast to 2.1% of PR76, while the fitting over 8 pure compounds provided 0.8% AAD (0.4% AAD for water).

α(T ) = exp 1 -c ln Tr 1 , (4.1.13)
b = 0.0778 RT c P c β(T ), (4.1.14)
β(T ) = 1 + c 2 (1 -T r ) , (4.1.15)
where c 1 and c 2 are the fitted parameters in the polynomial form of ω:

2 k=0 c i,k ω k , c i,k are
the given constants, that were calculated using the saturated vapour pressure experimental data of 53 pure compounds and the liquid density data of 49 pure compounds.

The PR78, SRK, GPR, HPR models were compared on vapor pressure, liquid density and vapor volume of pure compounds. In Table 4.2, we present the results of the selected pure components from [START_REF] Haghtalab | A modified peng-robinson equation of state for phase equilibrium calculation of liquefied, synthetic natural gas, and gas condensate mixtures[END_REF], where the overall AAD is given for 49 pure components. As expected, PR78 provides slightly better estimation of vapor pressure than SRK, then the refined results are obtained by GPR and HPR, since both of them modified 120 CHAPTER 4. EQUATIONS OF STATE AND FLUID PROPERTIES the α function and fitted it to the vapor pressure calculation. The additional adjustment on the co-volume parameter was made in HPR in contrast to GPR and HPR, that yields to the improved liquid density. However, the water density AAD is still about 18%.

Volume translation

The volume translated (VT) method is another way to impact on the volumetric properties. Péneloux et al. introduced the pseudo volume .16) where V is the volume, n i is the number of moles of component i, c i is the constant translation parameter or the correction of molar volume, which can be applied as an adjusted parameter for the volumetric properties. If the VLE is solved on pseudo volume, the pseudo partial molar volume should be introduced

Ṽ = V + i c i n i , (4.1
ṽi = ∂ Ṽ ∂n i T,P,n j ,j =i = v i + c i , (4.1.17)
and consequently, the pseudo fugacity coefficients is derived from the fugacity coefficient definition .1.18) If one writes the φφ phase equilibria by using the pseudo fugacity coefficients, the second term of the RHS of Eq. (4.1.18) appears symmetrically, that confirms the preserved consistency of the method. In the φ-γ system, it can be directly integrated in the Poynting factor Eq. (3.1.10) inasmuch as the impact of the translation parameter is expressed by the separate term.

ln φi = P 0 ṽi RT - 1 P dP = ln φ i + c i P RT . ( 4
P = RT ṽ -b i - a i (T ) ṽ(ṽ + b i ) + b i (ṽ -b i ) , (4.1.19)
where ṽ = v + c is the pseudo molar volume defined in Eq. (4.1.17). Similar to [START_REF] Péneloux | A consistent correction for redlich-kwong-soave volumes[END_REF], the correction volume function c is defined through the critical compressibility factor: .20) where δ and χ are the constant parameters fitted for each compound at the temperature range T r < 0.8, the data for 44 pure compounds can be found in [START_REF] Ahlers | Development of an universal group contribution equation of state: I. prediction of liquid densities for pure compounds with a volume translated peng-robinson equation of state[END_REF]. After optimization procedure, the average δ and the linear correlation between χ and z c were defined as δ = 0.252 and χ = 2.5448z c -0.4024, that leads to the generalized formulation of c (gVTPR):

c = δ RT c P c (χ -z c ), (4.1
c(z c ) = 0.252 RT c P c (1.5448z c -0.4042). (4.1.21)
For T r < 0.8, the molar volume v is overestimated by VTPR, it is therefore adjusted by c(T ) function which is derived from its value at the critical point and the temperature adjusting term (T-VTPR):

v = v V T P R -c(T ) = v V T P R -c c ζ(T r ), (4.1.22)
c c = (0.3074 -z c ) RT c P c , (4.1.23) ζ(T r ) = 0.35 0.35 + (η|T r -α(T r )|) γ , (4.1.24)
where η and γ are the constant parameters fitted for each compound. Then they were also generalized as the functions of z c on the same compounds except water, chloroform and methanol, that do not lie on the η curve as other compounds (gT-VTPR): In order to establish the generalized models gVTPR and gT-VTPR, the critical compressibility factor z c is thus required in addition to the standard P c , T c and ω.

The comparison of PR78, SRK, VTPR, gVTPR, T-VTPR, gT-VTPR models on the liquid density in the sub-and near critical regions is presented in Table 4.3, the results are taken from [START_REF] Ahlers | Development of an universal group contribution equation of state: I. prediction of liquid densities for pure compounds with a volume translated peng-robinson equation of state[END_REF]. The dominant position of (g)T-VTPR is distinct and these models can be advantageously applied for the liquid density estimation.

Nevertheless, the further usage of these models at high pressure and at the critical point leads to unphysical values, because isotherms intersect due to the introduced additional volume shift. Hence, the reliability strictly depends on the PT region and on the purpose of their usage. The (g)VTPR models provide less accurate results than (g)T-VTPR, but still better in total than PR78 and SRK.

Excess Gibbs energy and mixing rules

The mixing rules for a and b defined by Eqs. (4.1.6) and (4.1.8) can be differently derived.

Assuming that the excess volume is zero V E = 0 and consequently, the excess Gibbs energy equals the excess Helmholtz energy G E = A E , the latter can be derived from the equation of state and from the mixing rules, that was firstly made for van der Waals EOS. Conversely, applying the excess Gibbs energy model and introducing it in the EOS, the mixing rules can be found, called the G E mixing rules. It improves the VLE estimation of mixtures in sub-/critical regions and makes the parameters a and b temperature dependent. The G E models, assuming G E at infinite pressure or at zero pressure, imply the use of UNIFAC method [START_REF] Fredenslund | Group-contribution estimation of activity coefficients in nonideal liquid mixtures[END_REF] for G E and the activity coefficients calculations. This approach was applied for VTPR model (Ahlers and Gmehling, 2002a,b) with the GCM and the generalized α function of TPR [START_REF] Twu | A new generalized alpha function for a cubic equation of state part 1. peng-robinson equation[END_REF].

Another method was discovered by [START_REF] Peneloux | Excess functions and equations of state[END_REF] through the van Laar equation and the standard mixing rules, that yields

a b = i x i a i b i - A E C EOS , (4.1.27
) .30) where to be rigorous, A E is the excess Helmholtz, E ij are the energy binary parameters.

A E = 1 2 i j b i b j x i x j E ij i x i b i , (4.1.28)
k ij = E ij -(δ i -δ j ) 2 2δ i δ j , (4.1.29)
δ i = a i (T ) b i , (4.1
This method was employed to PR78 in [START_REF] Jaubert | Vle predictions with the pengrobinson equation of state and temperature dependent kij calculated through a group contribution method[END_REF], called the predictive PR78 (PPR78). This evolves into the temperature dependent binary energy interaction parameters E ij (T ) and hence, k ij (T ) that are expressed by the GCM. Considering that each molecule i is composed of N g,i groups (e.g., CH 4 , CH 3 , CH 2 ) and that the group k is presented twice, then one can define α ik as a fraction of group k, α ik = 2/N g,i . Knowing the distribution {α ik }, {α jk } of two molecules i and j, the binary interaction parameters k ij are given by .31) where N g is the total number of groups, A kl and B kl are the group constants, defined in [START_REF] Qian | Predicting the phase equilibria, critical phenomena, and mixing enthalpies of binary aqueous systems containing alkanes, cycloalkanes, aromatics, alkenes, and gases (N 2 , CO 2 , H 2 S, H 2 ) with the PPR78 equation of state[END_REF] for 20 groups, covering the petroleum fluids and natural gases. {A kl } and {B kl } form the symmetric matrices with the zero diagonals. a i (T ) and b i correspond to those of PR78, Eqs. (4.1.2) and (4.1.5), respectively. The method allowed to improve the prediction of the phase equilibria, the dew, bubble and even critical curves for binary mixtures and synthetic fluids.

E ij (T ) = - 1 2 Ng k=1 Ng l=1 (α ik -α jk )(α il -α jl )A kl 298.15 T B kl A kl -1 , (4.1

Cubic plus association

In order to take into account the hydrogen bonding and not only spherical molecules, the cubic EOS are combined with the Wertheim's term [START_REF] Kontogeorgis | An equation of state for associating fluids[END_REF][START_REF] Shinta | Equation of state representation of aqueous mixtures using an association model[END_REF][START_REF] Hajiw | Hydrocarbons-water phase equilibria using the cpa equation of state with a group contribution method[END_REF] that allows to describe the thermodynamic properties of pure associating components (e.g., water, alcohols) and its mixtures. Unifying the cubic PR78 and associative term results in the CPA-PR [START_REF] Hajiw | Hydrocarbons-water phase equilibria using the cpa equation of state with a group contribution method[END_REF]:

P = RT v -b i - a i (T ) v(v + b i ) + b i (v -b i ) - 1 2 RT v 1 + ρ ∂ ln g ∂ρ N i=1 x i A i (1 -X A i ) , (4.1.32)
where g is the modified radial distribution function of [START_REF] Kontogeorgis | An equation of state for associating fluids[END_REF], ρ is the molar density, x i is the mole fraction of the component i, X A i is the mole fraction of the molecule i not bonded at a site A. X A i depends on two parameters additional to PR78the association energy ε A i B i and the association volume β A i B i , that can be neglected for the non-associative components. The model demonstrated promising results of the volumetric properties on the light, medium hydrocarbons-water systems and also CO 2 and H 2 S water mixtures1 but at cost of computational time making it not suitable for the reservoir modeling.

The choice of the EOS model should be made carefully when one looks both for the phase equilibria of pure compounds and complex mixtures at sub-/critical regions and for the volumetric properties of both phases. The conventional PR78 can be applied for the solubility and gas density evaluation. The models with modified α function (TPR, CPR, GRP, HPR) provide rather better solubility results than those of PR78. The PR78 model with the binary interaction parameters predicted by PPR78 is able to reproduce the phase equilibria of mixtures, its dew, bubble and critical points. The improvement of liquid density can be expected by modifying the mixing rule of the covolume parameter.

Although, the volume translation methods ((g)VTPR, (g)T-VTPR) significantly increase the liquid density accuracy both of pure compounds and of mixtures, Sec. 4.2, the usage at high pressures should be avoided [START_REF] Ahlers | Development of an universal group contribution equation of state: I. prediction of liquid densities for pure compounds with a volume translated peng-robinson equation of state[END_REF]. The CPA-PR predicts the water systems, but it is not acceptable near the critical point. In order to improve the vapor pressure estimation, the α function can be modified. In MMF modeling, applying the different EOS models for each sub-problem is possible but entails inconsistencies. Therefore, depending on the principal goal of the modeling and PT region, the corresponding equations of state should be chosen.

EOS comparative study on volumetric properties

Pure compounds

The EOS volumetric properties study of pure compounds was fulfilled for

• CO 2 , H 2 S, CH 4 at 100 • C, 100 -350 bar (T r > 1);

• CO 2 , Ar, O 2 , N 2 , SO 2 at 50 • C, 100-250 bar (T r > 1 except SO 2 for which T r = 0.8.
The PR78, TPR, VTPR, T-VTPR, GERG-2008 models and NIST database were chosen at this stage of the study. The GERG-2008 EOS [START_REF] Kunz | Technical monograph: The gerg-2004 wide-range equation of state for nature gases and other mixtures[END_REF] was employed by the thermodynamic properties software REFPROP [START_REF] Lemmon | Nist reference database 23: reference fluid thermodynamic and transport properties-refprop, version 9.1[END_REF]. The GERG-2008 results correspond to the NIST data for the pure compounds, so they are omitted. In the T-VTPR model [START_REF] Ahlers | Development of an universal group contribution equation of state: I. prediction of liquid densities for pure compounds with a volume translated peng-robinson equation of state[END_REF], H 2 S was not taken into account.

Thus we used the generalized formulation for the temperature dependent VT parameter c(T ) Eq. (4.1.22), resulting in up to 19% AAD. The H 2 S density can be well described (0.4% AAD) by adjusting η = 5 and γ = -0.2, Eqs. (4.1.25) and (4.1.26) respectively.

Nevertheless, considering the generalized form of γ > 0 ∀z c Eq. (4.1.26), H 2 S is out of the model conception. At 100 • C and 100 -150 bar, the VTPR model reaches 7% AAD, hence the PR78 is preferable. On the whole pressure range 100 -350 bar, the PR78 is recommended for the volumetric properties evaluation of the selected pure compounds. The comparative study of different EOS on the liquid density of water binary systems at high T P conditions was carried out, considering the CPA-PR model [START_REF] Hajiw | Hydrocarbons-water phase equilibria using the cpa equation of state with a group contribution method[END_REF] as a reference. The molar volume of CPA is located between VTPR and T-VTPR with the deviation of 4% and 3%, respectively, Fig. 4.9. TPR provides a higher v m of 21% deviation and consequently a lower density of 17% deviation that is about of the same degree as AAD for water density given by the PR78, CPR, GPR and HPR (the α modified PR EOS), compare with Table 4.2. The TPR and volume translation EOS, based on the PR78 with the same mixing rules and binary coefficients, yield identical excess molar volume v E m , less than that of CPA-PR for all three binaries, see Fig. [START_REF] Mccain | Reservoir-fluid property correlations-state of the art[END_REF]. [START_REF] Adams | Equations of state for basin geofluids: algorithm review and intercomparison for brines[END_REF] displayed the reliability of McCain method providing the comparative study of five analytical density models for the water with the NaCl content up to 30 wt% at high TP conditions and for the formation waters of the Alberta Basin at S < 27 wt% at STP conditions; even though the underestimation of the measured densities increases with S, reaching about 2% at S = 30 wt%.

Viscosity

Modeling CO 2 Altunin (1975) developed the thermodynamic models for CO 2 and in particular, the viscosity model as a function of ρ and T for the range 220 -1300 K and up to 1200 bar (Altunin, 1975, Eq. 9.33): where ρ r = ρ/ρ c and T r are the reduced gas density and temperature, a ij are the constant parameters given in Table 4.13, the density can be provided by EOS. The correlation of Altunin viscosity model with the database of [START_REF] Vargaftik | Handbook of physical properties of liquids and gases-pure substances and mixtures[END_REF] presented in [START_REF] Garcia | Fluid dynamics of carbon dioxide disposal into saline aquifers[END_REF] results in the relative deviation < 2% at P > 150 bar, while it reaches 8% at 50 -70 formulation showed 2% of deviation in the gas and supercritical regions [START_REF] Garcia | Fluid dynamics of carbon dioxide disposal into saline aquifers[END_REF], the latter can be chosen for the CO 2 viscosity modeling at corresponding T P conditions.

µ CO 2 = µ 0 exp   4 i=1   1 j=0 a ij T j r   ρ i r   , ( 4
Modeling H 2 O, {H 2 O + CO 2 } and {H 2 O + NaCl + CO 2 }
The [START_REF] Islam | Viscosity models and e↵ects of dissolved co2[END_REF] (4.3.12) 4.3.13) where a i , b i , c i , d i and a r i , b r i are the constant parameters, see Tables 4.14 and 4.15. The viscosity is strongly temperature dependent and slightly pressure dependent, the CO 2 content impact decreases with the temperature, see Fig. 4.12.

µ {H 2 O,CO 2 } = µ H 2 O µ r ,
µ r = 1 + 2 i=1 a r i x CO 2 i 1 i=0 b r i T 1 , (
For the {H 2 O+NaCl+CO 2 } viscosity, Islam and Carlson derived the following equation 4.3.15) (4.3.16) where µ M D r is the relative viscosity, m is the molality of salts, A, B, C are the temperature dependent functions: The constant parameters A i , B i , C i for {H 2 O, NaCl} system are given in Table 4.16. NaCl},Eq. (4.3.15). problem. Then, we compare the numerical results of a 2D radial CO 2 injection with some analytical solutions. Finally, the model is used to investigate the effect of the presence of impurities in the injected stream.

µ {H 2 O,NaCl,CO 2 } = µ {H 2 O,NaCl} 1 + 4.65x
µ {H 2 O,NaCl} = µ H 2 O µ M D r , (
ln µ M D r = Am + Bm 2 + Cm 3 ,
A = 2 i=0 A i T 2 , (4.3.17) B = 2 i=0 B i T 2 , (4.3.18) C = 1 i=0 C i T 2 . (4.3.19)
A 0 = -0.21319213 B 1 = -0.27292263 × -3 A 1 = 0.13651589 × 10 -2 B 2 = 0.20852448 × -6 A 2 = -0.12191756 × 10 -5 C 0 = -0.25988855 × -2 B 0 = 0.69161945 × 10 -1 C 1 = 0.77989227 × -5

Methods

The numerical method described in Sec. volved in CCS, Sec. 4.2. In Al-Siyabi (2013), the PR78 was studied on eight CO 2 binary systems including O 2 , N 2 , SO 2 , during which the PR78 ADD for these three systems on the saturation pressure and on the vapor composition was evaluated as 3.2% and 1.9%, respectively. These EOS models were fitted to the experimental data acquired on binary systems. They can be applied on ternary or more systems, even though their accuracy could then be questioned, and would require datasets, not yet existing.

5.2. GAS DYNAMICS AND CHROMATOGRAPHIC PARTITIONING 149

Gas dynamics and chromatographic partitioning

The chromatographic partitioning of injected gas was investigated experimentally in the work of Bachu and Bennion, in which the authors provided the sensitivity study on the impact of gas impurity concentration (2, 5, 30% of H 2 S), PTS conditions (61 • C and 13.5

MPa; 25 • C and 6 MPa; 119 ppm) and of admixture type in the injected stream (H 2 S, N 2 , SO 2 , CH 4 ). In the last stage, the laboratory results revealed a key factor of the impurity solubility order regarding to that of CO 2 . The breakthrough time of the gas compounds and their profiles highly depend on the solubility ratio between gases. In Fig. 5.1, the mole fraction of the impurity gas H 2 S, N 2 , SO 2 and C 1 (methane) are normalized by the relevant initial mole fraction of the injected gas (5% for all four tests) in order to emphasize the scale of the gas partitioning effect due to the different solubility relations. The less soluble the admixture gas is, the earlier and higher breakthrough is expected. The impurity gases present in the injection stream after the capture processes can then be classified by the solubility. For example, we can write them in the decreasing order of the solubility (Henry's ( Meanwhile, convective fluxes start to form over the whole length of the capillary fringe.

However, the gas motion is dominant and prevents density fingering. Once the gravity perturbations in the CO 2 rich zone are large enough, the system becomes unstable to the convective fluxes on the transition zone. Its impact can be seen in Fig. 5.11, where starting from the onset time point the mass transfer increases. The onset time, t c = 92 d, can be adjusted assuming c 0 = 6317, Eq. (5.3.5). Due to the equilibration stage, the onset time is delayed in comparison with the previous results for the single phase problem [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Pau | Highresolution simulation and characterization of density-driven flow in co 2 storage in saline aquifers[END_REF]. For the relative onset time, t ct eq = 53 d, c 0 = 3623. This value can be overestimated due to the grid dimension. Additional grid refinement would be necessary to define the onset convergence point that was numerically demonstrated by high resolution modeling of single-phase problems [START_REF] Pau | Highresolution simulation and characterization of density-driven flow in co 2 storage in saline aquifers[END_REF] and also of twophase problems with transient and stagnant capillary fringe [START_REF] Elenius | Interactions between gravity currents and convective dissolution[END_REF]. These results correspond to the grid with h x,z = 1/16 m. Using the coarser grid provides high underestimation. For example, choosing the grid of 1/16 allow to model 18 fingers during the early convective time, while the grid of 2 m can produce only 1 finger after 5 y.

The CO 2 (aq) rich fingers grow down and the mass transfer flux reaches a stabilization stage, approximately at 157 d, Figs. 5.11 and 5.12. During this period, the average mass transfer flux is 9.73 × 10 -7 kg/m 2 /s that is of the same order of magnitude as the results

of [START_REF] Pau | Highresolution simulation and characterization of density-driven flow in co 2 storage in saline aquifers[END_REF].

The convective mass transfer is 174 kg of CO 2 during the stabilization period, Fig. 5.12.

The density deviation is ∆ρ l = 18.28 kg/m 3 . The analytical estimation Eq. ( 5 Hence, in addition to the convective motion in the aqueous phase, the instabilities growing on the rough contact surface between two phases imply corresponding perturbations in the gas phase, evolving the gas mixing, Fig. 5.15. It is worth mentioning that the convective mixing in both phases was also observed at the reservoir scale, Sec. 3.3.3.

We believe that the moving uneven transition zone affects the convection both in the aqueous and in the gas phases. The constant gas source also impacts to the process:

supplying fresh CO 2 (g) allows for the continued development of new fingers, compare

Initial

Conditions The modeling of CO 2 injection at the reservoir scale provides an analysis of the gas plume evolution and undergoing processes. However, using the coarse grid results in lower convective fluxes and consequently, in lower dissolution rate. We study the convective underestimation on the radial geometry applying numerical modeling and analytical formulations.

Ω = [0, 10] × [0, H] [m 2 ] Ω w = [0, 10] × [0, H -h z ] [m 2 ] Ω g = [0, 10] × [H -h z , H] [m 2 ] h z = 1/16, 1/8, 1/4, 2 [m], see Sec. 5.3.1 Ω w : p l0 = p atm + ρ l g(z D -z) [P a] S g0 = 0 Ω g : p l0 = p atm + ρ l g(z D -z) [P a] S g0

Problem description

The supercritical CO 2 is injected through a well in a vertically confined long aquifer with a height of 100 m at a depth of z D = 1100 m. The injection rate is constant, 10 kg/s; the well radius is 0.3 m; the injection is performed over the whole area of the well during 30 y. Supposing a homogeneous medium, the 2D axisymmetric geometry is chosen. The fluid and matrix parameters are identical to those of the convective mixing problem, Table 5.2.

The length of the aquifer is 10 km. Variable cell length from 2 to 1000 m allows to reduce the grid dimension to 33, 750 cells.

Numerical results

The injected CO 2 (g) rises up and then spreads along the top of the aquifer, Fig 5 .16. The maximum expansion can be followed by the radius of this plume at the top of the aquifer, Fig. 5.17. The evolution is proportional to √ t, consistent with the axisymmetric geometry.

Notice that gas displacement is the dominant mechanism: the dissolution rate is weak, so that the dissolved CO 2 area expansion is limited by the gas displacement.

In the absence of convection, the amount of dissolved CO 2 is proportional to the amount of injected CO 2 (g), Fig. 5.18. Indeed, the dissolution occurs only where CO 2 (g) is available, and at equilibrium. On the contrary, the convective flux increases CO 2 dissolution by The CO 2 rich fingers directed downward provoke the counter water stream. The induced circulation of fresh water enhances the dissolution. In this period, the convection impact is linear with time. Then at about 16 y, we can observe the later convective mixing period, Fig. 5.18, when the liquid under the current contains the dissolved CO 2 . At 30 y, 52% of injected CO 2 remains in the gas phase. So, almost half of the injected CO 2 is trapped by dissolution. Without the convective dissolution, the rate of CO 2 (g) is nearly constant, 79%.

Gas current spreading

The numerical results provide the lower convective dissolution at the reservoir scale, so we give the order of this underestimation. The radius of the gas current at the top of the aquifer neglecting dissolution can be expressed by following [START_REF] Nordbotten | Injection and storage of CO 2 in deep saline aquifers: Analytical solution for CO 2 plume evolution during injection[END_REF]: 5.4.1) where λ g and λ w are the mobilities of gas and water, Q is the volumetric injection rate.

R = λ g Qt φπλ w H , (
Assuming λ α = 1/µ α , the analytical estimation Eq. 5.4.1 corresponds well to the numerical results of the plume radius that include the dissolution and diffusion processes, Fig. 5.23.

This good agreement entails that neither dissolution nor diffusion has significant impact on the radius of the gas spreading.

Since the CO 2 (aq) rich current follows the gas spreading (Fig. 5.17 grid. We thus deduce that the 2D radial modeling can be efficiently used reducing the grid dimension and enhancing the convective mixing in comparison to the 3D modeling.

CO 2 with impurities and convective mixing

Preliminary observations of the current spread and the impact of impurities:

• From the chromatographic partitioning analysis Sec. 5.2, we conclude that species less soluble than CO 2 promote the advance of the plume; more soluble species reduce the radius. These results have been observed at the field scale, on the pilot in Tongliao, [START_REF] Wei | Geochemical impact of aquifer storage for impure co 2 containing o 2 and n 2: Tongliao field experiment[END_REF].

• The difference between the injected fluid density and that of the aquifer water (∆ρ ph = ρ f wρ inj ) influences the distance of the plume spread: the smaller difference implies smaller buoyancy forces, hence the smaller spread radius and vice versa, by Archimedes' principle.

• The difference between the increased density of the aquifer/formation water with the solutes and the initial aquifer water density (∆ρ l = ρ dissρ f w ) is responsible for the density driven convective flow:

-The higher the density difference, the higher the rate of convective fluxes should be and vice versa, see the experimental and analytical results of [START_REF] Kneafsey | Laboratory flow experiments for visualizing carbon dioxideinduced, density-driven brine convection[END_REF][START_REF] Macminn | Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers[END_REF][START_REF] Neufeld | Convective dissolution of carbon dioxide in saline aquifers[END_REF].

-The higher the density difference, the smaller the wavelength and the earlier the onset time are expected [START_REF] Ennis-King | Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations[END_REF][START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF]. The convective processes can be also observed in the gas phase. The impact of H 2 O(g) is negligible.

The higher convective mass transfer was expected in the SO 2 scenario [START_REF] Li | A numerical study of the impurity effects of nitrogen and sulfur dioxide on the solubility trapping of carbon dioxide geological storage[END_REF][START_REF] Ziabakhsh-Ganji | Sensitivity of the co 2 storage capacity of underground geological structures to the presence of so 2 and other impurities[END_REF] as the density difference ∆ρ l is almost two times bigger than that of the pure CO 2 and Air scenarios. However, due to the chromatographic partitioning a CO 2 rich front develops. So the convective transfer is similar to that of the pure CO 2 injection. As a result, the curves of convective impact in the pure CO 2 and SO 2 cases are very close, Fig. 5.30.

For the Air scenario, even though the liquid density deviation ∆ρ l is smaller then that of the pure CO 2 case, the gas density and chromatographic partitioning of Air scenario are favorable to the gas current advancing, hence to the pure dissolution. Also, it increases the spread area slightly enhancing the convective impact. That explains a low convective impact and high dissolution rate at the same time.

In the light of these results, we conclude that the high interaction between the chromatographic properties of the mixture on the one hand and its convective potential on the other, requires careful modeling two-phase compositional flow and thermodynamics of mixtures. In order to establish further tendency of convective contribution, the injection time should be extended, the grid dimension and refinement is necessary. Les nouvelles fonctionnalité d'HYTEC ouvrent un large éventail d'applications des géosciences ainsi que de possibilités et questions relatives aux développements futurs du code.

Conclusions

The present work provides a study of reactive transport coupling with multiphase multicomponent compressible flow. A new coupling method was proposed, relying on a operator splitting approach, so that its integration into an existing reactive transport code is facilitated. The method was implemented in the HYTEC code.

Before this work, HYTEC possessed saturated or unsaturated (Richards' equation) flow capabilities, coupled with aqueous and gaseous (diffusive) transport and chemistry.

A two-phase flow module was developed, using a fully implicit (with an analytical Jacobian) pressure/saturation algorithm. Then, the operator splitting methodology was implemented, so that multiphase multicomponent compressible flow with complex fluid properties is now coupled with multiphase reactive transport. Applying the formulation of phase conservation allowed to minimize the dimension of flow problem and make its structure composition independent. The coupling method conserved the existing functionality and flexibility of the code. The equations of state were used for the phase equilibria and prediction of volumetric properties, whose comparative analysis highlighted advantages of several models according to the composition and system conditions.

The verification of the code was performed by modeling problems with analytical and literature results; in addition, a benchmark exercise for the multiphase flow and reactive transport modeling community was proposed to verify and validate the coupling methods.

The developments of this work were applied for modeling the CO 2 underground storage.

The proposed coupling allows to investigate a complex gas injection (CO 2 and impurities) and to reveal new behaviors of gas and liquid dynamics -high interaction between density convective flux and chromatographic partitioning. We investigated the impact of the injected gas composition on the velocity of gas saturation front, density driven fingering and its effect on the overall gas dissolution. This simulation paves the way for future work, including hydrodynamic heterogeneities or pronounced effect of the host-rock mineralogy. 6.2. FUTURE WORK 183

Future work

The new functionality in HYTEC opens a wide array of applications, as detailed above.

It is also obvious that the present work unveils a whole world of questions relative to the future development of the code.

Applications

• Behavior of complex gas for geoscience applications: CO 2 storage but also new forms of storage for the energy sector, like H 2 or high frequency natural gas storage. The behavior of residual air during the saturation stage of nuclear waste facilities is also interesting: although the pressure temperatures are less extreme, O 2 consumption, CO 2 production or consumption over acidic or basic fronts in the storage, H 2 production by the corrosion of steel containers, make for complex gas evolution in the system.

• The impact of several important parameters has already been addressed here: differential solubility, density. The role of capillary pressure should also be investigated in several relevant applications, particularly in heterogeneous media. Capillary pressure discontinuity indeed raises some numerical developments [START_REF] Van Duijn | The effect of capillary forces on immiscible two-phase flow in heterogeneous porous media[END_REF].

Simplifying hypotheses

Along the development of our method and its implementation, several simplifying hypotheses were chosen. An obvious reason is overall limitations of the code: for instance the chemical module CHESS relies on the presence of water, so that a robust algorithm for water disappearance is pointless. Some hypotheses were considered because they simplify the independent modules without effect on the coupling method (e.g., isotropy).

• An isotropic formulation was chosen, consistent with the usual finite volumes Voronoï scheme in HYTEC. A new vertex centered scheme is now available in HYTEC for the saturated mode. Its extension for the two-phase flow is relatively straightforward, although averaging between adjoining vertices should be carefully addressed (particularly for capillary pressure and relative permeability). This implementation would have no impact on the coupling scheme, due to the operator splitting approach.

• Including the feedback of chemistry on the hydrodynamics is still a challenge for the reactive transport community. The impact of mineral volume change on fluid flow can be addressed by a volumetric source term that can be incorporated in the coupling equations of OS structure. Another problem is the impact of the porosity change on the hydrodynamic properties, intrinsic permeability, relative permeability and capillary pressure. No global satisfying method is available so far, apart from specifically designed empirical laws.

• The creation of a gas phase, due to chemical reaction, in an initially saturated zone should also be addressed with care. A comparison between potential equilibrium pressure and actual fluid pressures is a first indicator that can be applied. The communication with the flow module must again be addressed with care. This problem is particularly important in the propagation of acidic plumes in more or less carbonated environments (e.g., in situ recovery), or in non-isothermal modeling.

Extension of the physical problems

• Dual porosity leads to different discretization, homogenized models of flow system, (Masson et al., 2014a;[START_REF] Jurak | Fully homogenized model for immiscible incompressible two-phase flow through heterogeneous porous media with thin fractures[END_REF], and opens the possibility for future researches in coupled multiphase flow and reactive transport in fractured media. The chemical impact in the flow in fracture can be stirring in connection with the RT and variable porosity/permeability due to the geochemical reactions.

• A non-isothermal module is available in HYTEC: for saturated regime only, the heat and flow solvers are sequentially coupled, then the Darcy velocity and the temperature fields are used within the reactive transport coupling (e.g., modification of the solubility and kinetics constants). Non-isothermal multiphase flow raises an issue of coupling method. The next options should be studied: FIM/IMPES flow + explicit/implicit heat. Applying the explicit scheme entails the definition of associated CFL criteria. In terms of perspective coupling with the reactive transport, the chemical impact might be a key factor in choosing explicit or implicit method.

• Coupling free flow and saturated and two-phase flow in porous medium implies introducing Navier-Stokes equations and Beavers-Jeseph interface condition (Rybak 6.2. FUTURE WORK 185 et al., 2015;[START_REF] Brenner | Coupling of a two phase gas liquid compositional 3D Darcy flow with a 1D compositional free gas flow[END_REF].

Numerical methods

• Applying a strong coupling (iterative procedure of flow, transport, chemistry) needs a comparative study with global implicit method in order to establish its capabilities to reproduce the chemical impact on the flow.

• Primary variables of multiphase flow system: once CHESS can handle the drying out process, the variable switching is unavoidable and it should be generalized.

• The AIM is of great importance for the large scale modeling since it can significantly reduce the nonlinear system. (The IMPES model is applied for decoupled formulation with pressure equation Sec. 2.2.1.) The linear solvers and preconditioners such as CPR [START_REF] Cao | Parallel scalable unstructured cpr-type linear solver for reservoir simulation[END_REF] can be observed for further development. Moreover, in some applications, the coupling with single-phase flow (domain decomposition) also increases computational efficiency. The flexible structure of the proposed coupling method allows different both formulations and techniques. Abstract: Human activity in the subsurface has been rapidly expanding and diversifying (waste disposal, new mining technologies, high-frequency storage of energy), while the public and regulatory expectations keep growing. The assessment of each step of underground operations requires careful safety and environmental impact evaluations. Multiphysics modeling provides an effective way to understand and predict the behavior of such complex systems at different time and space scales. This work aims at incorporating a compressible multiphase flow into the conventional reactive transport framework by an operator splitting approach. A multiphase flow module was therefore implemented in the reactive transport software HYTEC using a new approach, developed to fully couple multiphase multicomponent compressible flow and description of the fluid thermodynamic properties with existing reactive transport codes. An alternative method based on the sequential iterative approach was then invented to separately model multiphase transport by means of numerical reaction rate terms. The coupled system was verified and compared with other coupling formulations for accuracy and computational performance. A benchmark exercise of modeling a natural gas reservoir was proposed to investigate the coupling methods of reactive multiphase simulators, as well as their numerical capacity and characteristics.

Nomenclature

The final part of this work concerns the numerical and analytical modeling of carbon capture and storage and the impact of convective mixing and gas impurities in the injected stream in particular. The proposed coupling allowed to study complex gas injection scenarios, modeling a system including but not limited to four gases, and to reveal new behaviors of gas and liquid dynamics, specifically the interaction between density convective flux and chromatographic partitioning.

9 1

 91 figures is identical, the concentration range lies in [0; 1.24]. . . . . . . . . . . 5.6 Gas density in kg/m 3 (left) and saturation (right) of CO 2 , {CO 2 + N 2 } and {CO 2 + SO 2 } models at T = 1 min (top), 55 min (center), 1 and 10 d (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Gas mole fractions of impurity gases and pure CO 2 (left top), CO 2 concentration in molal (right top), gas saturation (left bottom) and the gas density in kg/m 3 (right bottom) of CO 2 , {CO 2 + N 2 }, {CO 2 + Ar}, {CO 2 + O 2 } and {CO 2 + SO 2 } models at T = 55 min. . . . . . . . . . . . . . . . . . . .

T r ≤ 1 ( 8 1 Introduction 2 CHAPTER 1 .

 18121 abbreviations used in the document are detailed in Nomenclature, p. 187.

4 CHAPTER 1

 41 . INTRODUCTION In the reactive transport codes, the number of components in the multiphase flow system is usually minimized, in contrast to reservoir simulators. Adding a component increases the dimension of nonlinear multiphase flow system and therefore, the computational complexity of the problem. Let us suppose an isothermal multiphase flow problem of N f fluid phases and N c components. Hence, there are • N c nonlinear mass conservation equations for each component, • 2N f constitutive relations (N f -1 capillary pressure equations, 1 equation of the saturation sum, N p equations of the mole/mass fraction sum),• N c (N f -1) phase equilibria equations, resulting in a problem of 2N f + N c N f equations. Based on the Gibbs phase rule, the number of intensive variables is N c -N f + 2. Then, considering two phases, there are N c unknowns, hence N c nonlinear equations. Using a classic global approach(Coats, 1980), we solve N c nonlinear mass conservation equations. The size of matrix arising from the linearization of this system, is N c N gr × N c N gr , where N gr is the grid dimension. Any additional component introduced to a multiphase flow problem indeed amplifies significantly computational complexity for high resolution modeling. Also, it requires the development of a thermodynamic module that is applicable to multicomponent systems and describes the fluid properties of real mixtures.

  ParisTech. The biogeochemical reactions are performed by the speciation code CHESS (van der Lee, 2009). The hydrodynamic module of HYTEC accounts for saturated, unsaturated flow. The code was validated in different benchmarks (De Windt

Chapter 2

 2 presents an overview of multiphase flow formulations. Then we describe the numerical implementation of one of them in HYTEC, supposing incompressible, immiscible, isothermal flow. The verification of the flow module is demonstrated on some reference models.Chapter 3 is devoted to the coupling of the multiphase multicomponent flow and reactive transport. We give an overview on the numerical methods for each part of the problem and of the existing coupling approaches. Next, a new method of modeling compressible multiphase flow and its coupling with reactive transport are detailed. Also, an alternative way to separate multiphase transport modeling is given, by means of the reaction rate terms for the sequential iterative reactive transport framework. The first numerical results are then shown. Finally, we present an application (transformed into a benchmark proposition) that supposes a tight connection between transport and geochemistry by modeling a natural gas reservoir (with three gases), submitted to Mathematics and Computers in Simulation.

  Chapitre 2 introduit les notions de base de l'écoulement multiphasique, nécessaire pour comprendre la physique, établir le problème mathématique et ensuite le modéliser. Une fois les formulations mathématiques de l'écoulement données, nous pouvons choisir l'une d'elles et poser le problème de l'écoulement diphasique isotherme avec les hypothèses simplificatrices à cette étape -notamment les fluides sont supposés immiscibles, incompressibles. Les forces diffusives et les réactions chimiques sont négligées dans cette partie. Les méthodes numériques sont détaillées : les schémas de discrétisation en espace et en temps, la linéarisation et le jacobien analytique. La vérification est établie sur la modélisation d'un problème parabolique avec une solution auto-similaire et d'un problème hyperbolique de référence. Les résultats numériques sont comparés à ceux obtenus par d'autres codes et montrent également la convergence et la précision du module construit.

Fig. 2 .

 2 1(a) illustrates two immiscible fluids in contact with a solid phase, where θ is the contact angle between the solid surface and the fluid-fluid interface curved under the surface tension. If θ < 90 • , the fluid is called wetting, otherwise, non-wetting. Fig.2.1(b) presents the interface in case of water/oil and gas/oil. If the two fluids are miscible (for instance, gases), there is no a sharp fluid-fluid interface, the fluids represent one phase. In this case, the phase consists of components: the system is multicomponent.

Figure 2 . 1 -

 21 Figure 2.1 -Determination of the wetting and non-wetting phases.

Figure 2 . 3 -

 23 Figure 2.3 -Relative permeability functions of S g (left) and its derivatives (right). k rl -van Genuchten, Eq. (2.1.11): m = 0.457, S lr = 0.3. k rg -Brooks-Corey, Eq. (2.1.10): λ = 2, S lr = 0.3, S gr = 0.05.

  s] the liquid phase zone, whereas the third region, R 2 /t ≤ 5 × 10 -7 [m 2 /s], Fig. 2.4(b) corresponds to the dryout zone. Fig. 2.5 displays the time evolution of gas saturation of both HYTEC and of TOUGH2-ECO2 during 10000 d. The apparent shift of saturation profile takes place due

  1e-07 1e-06 1e-05 1e-04 1e-03 1e-02 1e-01

Figure 2 . 4 -

 24 Figure 2.4 -1D axisymmetric problem: gas saturation S g as a function of R 2 /t. Results of (a) HYTEC and (b) TOUGH2-ECO2(Pruess et al., 2002). The dry-out zone (R 2 /t ≤ 5 × 10 -7 ) is not simulated by HYTEC.

Figure 2 . 6 -

 26 Figure 2.6 -1D axisymmetric problem: gas saturation S g as a function of R 2 /t. Results of the workshop Intercomparison of numerical simulation codes for geologic disposal of CO 2 (Pruess et al., 2004): IFP -Institut Français du Pétrole, SIMUSCOPP; IRL -Industrial Research Limited, an in-house version of TOUGH2 and CHEM-TOUGH; LBNL -Lawrence Berkeley National Laboratory, TOUGH2/ECO2; CSIRO -CSIRO Petroleum, an in-house version of TOUGH2/ECO2; ARC -Alberta Research Council, GEM.

Fig. 2 .

 2 Fig. 2.7(a) shows the model geometry. The water (wetting phase) is injected in the left bottom well during 750 d and the host oil (non-wetting phase) is produced from the right top corner. Table 2.2 gives the problem statement; the curves of relative permeability are plotted in Fig. 2.8.

Table 2 . 2 - 20 Γ 3 ]

 22203 Parameters of the hyperbolic problem. Initial Conditions Ω = [0, 300] 2 [m] µ = 20 µ = 1/20 Ω: p l0 = 1 [bar], S n0 = 0.999 p l0 = 1 [bar], S n0 = 0 N : ψ w = 0.096/ Γ N , ψ n = 0 [kg/m 2 /s] ψ w = 0, ψ n = 0.096/ Γ N [kg/m 2 /s] Γ D : p w = 1 [bar], S n = 0.999 p w = 1 [bar], S n = 0 ∂Ω\(Γ D ∪ Γ N ): ψ w = 0, ψ n = 0 ψ w = 0, ψ n λ = 2, S wr = 0, S nr = 0 Fluid Properties Density ρ w = 1000 [kg/m 3 ] ρ n = 1000 [kg/m Viscosity µ w = 10 -3 [P a • s] µ n = 20 × 10 -3 [P a • s] (a) Homogeneous media (b) Heterogeneous media

Figure 2 . 7 -Figure 2 . 8 -

 2728 Figure 2.7 -Hyperbolic problem: geometry.

Figure 2 . 9 -

 29 Figure 2.9 -Hyperbolic problem. Homogeneous media. Saturation S w isolines at 750 d: µ = 20 (left) minimum level -0.0011, contour interval -0.05; µ = 1/20 (right) maximum level -0.9999, contour interval -0.05. 40 × 40 (top), 80 × 80 (center), 160 × 160 (bottom) elements.

Figure 2 . 10 -

 210 Figure 2.10 -Isolines of saturation S w , 160 × 160 elements: (a) results of UG, minimum level -0.0001, the contour interval -0.05; (b) results of HYTEC, minimum level -0.0011, contour interval -0.05.

  homogeneous media. This test demonstrated a high significance of grid scale, especially in case of high µ, the convergence of the built module and the capability to model the discontinuous permeability.

Figure 2 . 11 -

 211 Figure 2.11 -Discontinuous permeability problem. Isolines of wetting phase saturation S w at 375 d, 160 × 160 elements: (a) results of UG, minimum level -0.0001, contour interval -0.05; (b) results of HYTEC, minimum level -0.0011, contour interval -0.05.

Figure 2 . 12 -

 212 Figure 2.12 -Discontinuous permeability problem. Saturation S w isolines at 375 d: µ = 20 (left) minimum level -0.0011, contour interval -0.05; µ = 1/20 (right) maximum level -0.9999, contour interval -0.05. 40 × 40 (top), 80 × 80 (center), 160 × 160 (bottom) elements.

1 . 1 ) 40 CHAPTER 3 .

 11403 has to be closed by the phase equilibria conditions. The equilibrium criterion, derived from the Gibbs-Duhem equation, dictates a minimum of the Gibbs free energy at constant temperature T , pressure P and composition. Consequently, it requires the equality of the chemical potentials of species i in all phases and then the COUPLING MMF AND RT equality of the fugacities:

  that can represent VLE for pure components and mixtures. In order to improve the volumetric properties near the critical point, the attractive part was modified by Redlich and Kwong in 1949 (RK) and it became temperature dependent through the α function appeared in P A . Pursuing the idea to better predict the VLE, the α function was elaborated by Soave in 1972 (SRK). Furthermore,Peng and Robinson (PR76) modified the attraction pressure P A to attain a better estimation of critical compressibility factor Z, liquid density and vapor pressure in comparison with those of SRK.

  .1.13)

  DUNE has recently been extended to RT modeling. The work of Hron et al. (2014) demonstrated the S(N)IA of MMF and RT applied to the growth and transport of microorganisms Escherichia coli. Later, in the DuMu X simulator based on the DUNE framework, Ahusborde et al. (2015); Vostrikov (2014) proposed a sequential non-iterative coupling of two-phase multicomponent flow with RT by separating the whole problem into dominant component flow and minor reactive transport subsystems applying SIA to the last.

  , based on the mixed FEM (MFEM), handles iterative coupling of flow and RT, presented by TRCHEM module. The different flow formulations are available: air-water, oil-water, black-oil, EOS compositional models. The FI, iteratively coupled IMPES/IMPEC can be employed. The multiphase RT are connected with any flow model in IPARS and is time split into advection (Godunov), reaction (second order Runge-Kutta) and diffusion/dispersion (expanded MFEM).

  ;[START_REF] White | STOMP Subsurface Transport Over Multiple Phases: STOMP-CO2 and STOMP-CO2e Guide: Version 1.0[END_REF]. Solving flow system evolves into a nonlinear PDE set for mass conservation of water, air, CO 2 , CH 4 , volatile organic compounds and salt. The geochemical module ECKEChem uses a general paradigm of BIOGEOCHEM for modeling reactive chemicals in batch systems.

Figure 2 :

 2 Figure 2: 1D axisymmetric problem: gas saturation S g as a function of R 2 /t

  (a) 7 days, max(S g ) = 0.73, max(x l,CO2 ) = 0.0251 (b) 20 days, max(S g ) = 0.79, max(x l,CO2 ) = 0.0250 (c) 65 days, max(S g ) = 0.84, max(x l,CO2 ) = 0.0248

Figure 3 :

 3 Figure 3: 2D problem of CO 2 injection in a fully water-saturated domain: gas saturation from [0, max(S g )] and contours of CO 2 mole fraction x l,CO 2 at 0.005, 0.011, 0.016

90 CHAPTER 3 .

 903 COUPLING MMF AND RT

  2 injection in a fully water-saturated domain, Neumann et al. (2013) Neumann et al. (2013) posed the model of CO 2 injection in homogeneous media taking the gravity into account. The geometry and grid dimension are the same as those of Neumann et al., while the fluid properties and the solubility differ from the reference. Table 3.1 contains the used model parameters. The tortuosity is presented by the Millington-Quirk model. The density and viscosity are variable and their models were presented in Sec. 4.3. During 65 days, 2.25 kt of CO 2 was injected.

  .1.10): p b = 10 3 Pa, λ = 2, S lr = 0, S gr = 0 Fluid Properties Density ρ l : T-VTPR ρ g : PR78 Viscosity µ l : Islam-Carlson µ g : Altunin75 Diffusion coefficient D l = D g = 2 × 10 -9 [m 2 /s]

Figure 1 :

 1 Figure 1: Schematic representation of the differential solubility of gases in water (left) and of the longitudinal section of a deep and trapped reservoir of acid gases in contact with an underlying aquifer (the dashed line stands for the localization of the simplified modeled grid).

Figure 2 :

 2 Figure 2: Geometry: Reservoir, Contact and Aquifer zones (red, yellow and blue, respectively)

155 3 .

 3 Fluid dynamics: incompressible two-phase flow is tested first using two different saturation-diffusion laws, then the full two-phase compressible flow is proposed.

Figure 3 :

 3 Figure 3: FCase 1. Time evolution of H 2 S(g) [molal]

Figure 4 :

 4 Figure 4: FCase 1. Gas saturation. Incompressible flow.

Figure 5 :Figure 6 : 21 116CHAPTER 4 .

 56214 Figure 5: H 2 S(g) mole fraction of FCase 1 ideal gas (top) and FCase 3 Peng-Robinson EOS (bottom)

4. 1 .

 1 OVERVIEW OF EQUATIONS OF STATE 117

4. 1 . 2 α

 12 function: vapor pressureIn order to give a better description of vapor pressure and VLE, the modified α function were proposed by representing it in the polynomial or/and exponential form, by applying the different formulation for sub-and supercritical regions, T r ≤ 1 and T r > 1. Considering that the attractive forces are expressed through the a parameter and consequently, through the α function, the latter should obey the next requirements: 0 < α(T ) < C, C is the positive real constant, ∀ T ; α(T c ) = 1. Moreover, the interaction between molecules weakens with the increase of temperature, so lim T →∞ α(T ) = 0. Note that the PR78 α function Eq. (4.1.2) does not fulfill this condition. Finally, to keep the consistency of the whole thermodynamic system, α(T ) ∈ C 2 (T ). Twu et al. (TPR) generalized the α function, switched in the sub-and supercritical regions:

  Gasem et al. (GPR) proposed a modified α function to improve the vapor pressure prediction of pure components without switching it for sub-and near crititcal conditions as was given in Mathias and Copeman (1983): α(T ) = exp (A + BT r ) 1 -T C+Dω+Eω 2 E are the generalized constants. The GPR model displayed the vapor pressure with 1.1% AAD for 28 pure compounds (3.7 AAD for water) and the heat capacity with 3% AAD for methane, propane and nitrogen for sub-and supercritical zones. AAD is defined as 100/N N i=1 |X calc -X exp |/X exp . Coquelet et al. (CPR) combines both the exponential and polynomial form of α function α

  Haghtalab et al. (HPR) changed both the α and the co-volume b functions to improve the estimation of the vapor pressure and the volumetric properties:

  Péneloux et al. applied the VT method to SRK with the three formulations for c, among which the volume correction is the function of the critical compressibility factor z c : c(z c ) = 0.40768(RT c /P c )(0.29441z c ). All three models were tested on 233 compounds and exhibited the conspicuous results. Even though the c(z c ) approach was the second best due to the lack of the precise values of z c , it showed overall 3% AAD against 11% of SRK and 5% of PR76.Thereafter,[START_REF] Ahlers | Development of an universal group contribution equation of state: I. prediction of liquid densities for pure compounds with a volume translated peng-robinson equation of state[END_REF] proposed the volume translated PR model (VTPR) and temperature dependent VTPR model (T-VTPR). Basing on PR78, VTPR 4.1. OVERVIEW OF EQUATIONS OF STATE 121 model is formulated by following

126 CHAPTER 4 .

 1264 The next modified PR78 models were implemented in HYTEC: the TPR model (the modified α function) in combination with the (g)VTPR and (g)T-VTPR models (the liquid density correction) and the temperature dependent k ij parameters defined by the PPR78 model. EQUATIONS OF STATE AND FLUID PROPERTIES

CO 2

 2 Fig. 4.1 presents the density of gas and supercritical CO 2 at 50 and 100 • C and to 20 from 350 bar, provided by the PR78, TPR, VTPR, T-VTPR and NIST database. In general, the density is underestimated except the T-VTPR at P > 160 bar and the TPR at P > 260 bar. Table 4.4 shows the AAD of the relevant EOS models: for the pressure range 20 -300 bar and T = 50 • C, the PR78 and TPR yield about 3%, the VT models 4%; the accuracy of all 4 models improves at 100 • C. The AAD varies significantly with the pressure range of interest. So, at 50 • C, the PR78, TPR and T-VTPR models can be applied for the pressure range 100 -250 bar, yielding 3 -4% AAD. However, we should note that the AAD of PR78 and TPR exceed 8% at 120 -130 bar. At 100 • C, all four models give acceptable results within 4% AAD. H 2 S at 100 • C, 100 -350 bar Fig. 4.2 shows the supercritical density of H 2 S at 100 • C and 100 -350 bar modeled by PR78, TPR, VTPR, T-VTPR and GERG-2008 (NIST). There is an obvious discrepancy between the EOS results: the T-VTPR curve is close to that of NIST near P r = 1, then the deviation grows with the increasing pressure and reaches 19% ADD. The PR78 and TPR with 2% AAD are preferable for the range 100 -150 bar, while the VTPR provides 1% for higher pressures 250 -350 bar, see Table 4.2.

. 2 .

 2 CH 4 at 100 • C, 100 -350 bar All models give a relevant density of supercritical CH 4 at the chosen T P conditions,

Figure 4 . 1 -

 41 Figure 4.1 -Density of pure CO 2 at 50 and 100 • C and 20 -350 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

Fig. 4 .

 4 Fig. 4.3, within 3% of AADTable 4.6. The volume translation models maintain about 1% AAD at the whole pressure range. CO 2 , H 2 S, CH 4 at 100 • C, 100 -350 bar: At high temperature and pressure conditions, 100 • C and 250-350 bar, the VTPR and PR78 models provide the acceptable gas/supercritical density (4.6% AAD at maximum) of the pure compounds CO 2 , H 2 S, CH 4 .

N 2 ,Figure 4 . 2 -

 242 Figure 4.2 -Density of pure H 2 S at 100 • C and 100 -350 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

Figure 4 . 3 -

 43 Figure 4.3 -Density of pure CH 4 at 100 • C and 100 -350 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

Figure 4 . 4 -

 44 Figure 4.4 -Density of pure N 2 at 50 • C and 100 -250 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

Figure 4 . 5 -

 45 Figure 4.5 -Density of pure O 2 at 50 • C and 100 -250 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

Figure 4 . 6 -

 46 Figure 4.6 -Density of pure Ar at 50 • C and 100 -250 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

Figure 4 . 7 -

 47 Figure 4.7 -Density of pure SO 2 at 50 • C and 100 -250 bar provided by PR78, TPR, VTPR, T-VTPR and NIST data.

4. 2 .

 2 EOS COMPARATIVE STUDY ON VOLUMETRIC PROPERTIES

Fig. 4 .

 4 9 shows the molar volume and the corresponding excess molar volume for three bynaries {CH 4 +H 2 O}, {CO 2 +H 2 O}, {H 2 S+H 2 O} at 100 • C and 300 bar given by TPR, VTPR, T-VTPR, CPA-PR, GERG-2008 (REFPROP). The solubility of each component at the given T P was calculated by CPA-PR, then varying the component mole fraction, the liquid density was calculated up to the saturation point.

4

 4

  .

  9. Because, by definition, v E m shows the deviation of the real mixture molar volume from the molar volume of the ideal mixture, the molecules of the real mixture require less volume than those of the ideal mixture for all three binaries. The CPA-PR molar volume is closer to that of the ideal mixture, than the results of VTPR models. The v E m variation is about 5 × 10 -8 for CH 4 and CO 2 mixtures CHAPTER 4. EQUATIONS OF STATE AND FLUID PROPERTIES and slightly larger for H 2 S binary. The mass liquid density is displayed in Fig. 4.10. The CH 4 binary density decrease with the growing CH 4 mole fraction due to its molar weight, which is less than that of H 2 O. For the CO 2 binary, the density grows toward the CO 2 saturation except the GERG-2008 results. The H 2 S binary density declines by using the CPA-PR, GERG-2008 and VTPR, their v m slopes are greater than those of the T-VTPR and TPR, Fig. 4.9.In spite of 4% underestimation, the VTPR model can qualitatively predict the liquid density of the water binary systems with CH 4 , CO 2 , H 2 S at high P T conditions. For the system CO 2 + H 2 O, the T-VTPR produces a good estimation of the liquid density with an acceptable overestimation over a wide range of P T .

Figure 4 . 9 -

 49 Figure 4.9 -Molar volume (left) and excess molar volume (right) of water binaries with CH 4 (top), CO 2 (center), H 2 S (bottom) at 100 • C and 300 bar, predicted by GERG-2008 (REFPROP), CPA-PR, VTPR, T-VTPR.

Figure 4 . 10 -

 410 Figure 4.10 -Liquid density of water binaries with CH 4 (top), CO 2 (center), H 2 S (bottom) at 100 • C and 300 bar, predicted by GERG-2008 (REFPROP), CPA-PR, VTPR, T-VTPR.

  presents the available fluid density models in HYTEC, where the constant mode for both phases implies incompressible flow. The majority of density options stem from the EOS models, although the analytical formulation for brine/formation water is also possible by adjusting the measured fluid properties. The EOS description was given previously in Sec. 4.1.Table4.12 -List of the density models available in HYTEC.Gas density modelLiquid density McCain models use the total dissolved solids (TDS). The former simply yields ρ l (S, T ) = ρ w (T )(1 + S), (4.3.1) where S is the TDS in mass fraction. Despite its trivial definition, applying this model gives a rough estimation of the phase transfer that arises from the solid-liquid and additionally from the liquid-gas interaction. By assuming the primary NaCl content of formation water, McCain Jr. derived the liquid density ρ l (S, T, P ) based on fluid samples data from the numerous petroleum reservoirs: ρ l (S, T, P ) = ρ w (S) B w (P, T ) , (4.3.2) where ρ w is the density of formation water at standard temperature pressure (STP) con-140 CHAPTER 4. EQUATIONS OF STATE AND FLUID PROPERTIES ditions (25

Figure 4 . 11 -

 411 Figure 4.11 -Viscosity of CO 2 at 45, 50, 75 and 100 • C and 20 -300 bar predicted by Altunin model Eq. (4.3.8) and NIST data.

  model describes the H 2 O viscosity at the T P range 20 -105 • C, 0 -600 bar, Eq. (4.3.11), producing a maximum error 0.05 of the IAPWS data of 2009. Then, the authors proposed the formulation for {H 2 O + CO 2 } system, Eqs. (4.3.12) and (4.3.13). The error up to 1% occurs in the T P x range 20 -100 • C, 50 -450 bar and 0 -0.024 mole fraction. µ H 2 O = a 0 +

Fig. 4 .

 4 Fig. 4.13 displays the example of the HYTEC syntax used for the definition of the density, viscosity, relative permeability and capillary pressure models.

Figure 4 . 13 -

 413 Figure 4.13 -Example of the model definition of the fluid and rock-fluid properties in the HYTEC input file.

  3.2 is employed to model the two-phase flow. The density properties of gas and liquid mixtures are predicted by the PR78 and T-VTPR models, respectively; we chose the binary interaction parameters given by[START_REF] Li | Evaluating cubic equations of state for calculation of vapor-liquid equilibrium of CO 2 and CO 2 -mixtures for CO 2 capture and storage processes[END_REF], that were calibrated for calculation VLE of CO 2 with the impurities. The preliminary work on the liquid density was carried out: the comparison of PR78, T-VTPR, CPA models and GERG-2008 (REFPROP) for pure components and binary systems in-

Figure 5 . 1 -Figure 5 . 4 -

 5154 Figure 5.1 -Mole fraction of H 2 S, C 1 (methane), N 2 and SO 2 at the outlet normalized by the relevant mole fraction in the initially injected gas, as a function of the pore volume of the total injection[START_REF] Bachu | Chromatographic partitioning of impurities contained in a CO 2 stream injected into a deep saline aquifer: Part 1. effects of gas composition and in situ conditions[END_REF].

Figure 5 . 5 -Figure 5 . 6 -

 5556 Figure 5.5 -Aqueous composition of {CO 2 + N 2 } (left) and {CO 2 + SO 2 } (right) in molal at T = 1 min (top), 55 min (center), 1 and 10 d (bottom). The scale of the figures is identical, the concentration range lies in [0; 1.24].

Figure 5 . 8 -

 58 Figure 5.8 -CO 2 concentration in molal of CO 2 , {CO 2 + N 2 }, {CO 2 + Ar}, {CO 2 + O 2 } and {CO 2 + SO 2 } models at T = 1 min (top), 10 min (center) and 10 d (bottom).

  .3.2) yields 102 kg with α = 0.2. The α coefficient adjusted to the mass transfer flux is 0.342. Based on the reported mass fluxes by[START_REF] Pau | Highresolution simulation and characterization of density-driven flow in co 2 storage in saline aquifers[END_REF], the adjusted α is 0.47. The mass flux estimation of[START_REF] Elenius | Convective mixing influenced by the capillary transition zone[END_REF], that also takes into account the capillary fringe, can be fitted to the numerical results with the scaled coefficient 0.055. In[START_REF] Elenius | Interactions between gravity currents and convective dissolution[END_REF], the mass flux of the models with no transition zone or transient capillary fringe is lower than the estimation of[START_REF] Elenius | Convective mixing influenced by the capillary transition zone[END_REF] with the coefficient 0.03. The adjusted α = 0.342 and mass flux can be used to estimate and compare the convective flow at reservoir scale as the deviation of the values is low.162 CHAPTER 5. CCS, CO-INJECTED IMPURITIES AND CONVECTIVE MIXINGAfter the onset point, the contact surface between two phases undergoes high perturbations and finally becomes uneven; further, its deviation (maximum and minimum distance) increases with the time. Fig.5.10 shows that the maximum and minimum gas front vertical position diverges from the outbreak point, 141 d. Therefore, the convective development, Fig.5.13, leads to a heterogeneous gas front, Fig.5.14. The uneven gas surface was also exposed in[START_REF] Elenius | Interactions between gravity currents and convective dissolution[END_REF] by modeling a long term injection in the reservoir scale.

Figs. 5 .

 5 Figs.5.9 and 5.13. Elenius et al. (2015) have shown that both transient and stagnant capillary transition zones accelerate the onset time and convective mass transfer.

Permeability K 10 - 12 [m 2 ]Figure 5 . 10 -Figure 5 . 11 -Figure 5 . 12 -

 10122510511512 Figure 5.10 -Evolution of the gas front vertical position.

Figure 5 . 13 -

 513 Figure 5.13 -CO 2 (aq) concentration and the velocity field at 120 and 365 d, the length is 10 cm.

Figure 5 . 14 -Figure 5 . 15 -

 514515 Figure 5.14 -Gas saturation at 120 and 365 d, the length is 10 cm.

Figure 5 . 16 -Figure 5 . 17 -Figure 5 . 18 -

 516517518 Figure 5.16 -Gas saturation map in range [0, max (S g )]: max (S g ) = 0.7 at 120 d (top), max (S g ) = 0.73 at 16 y (center) and 30 y (bottom), [100, 2500] m.

5. 4 .

 4 MODELING CO 2 INJECTION 169

Figs. 5 .

 5 Figs.5.19, 5.20 and 5.21 illustrate the CO 2 (aq), HCO - 3 concentration and pH at 120 d, then 16 y, 30 y corresponding to the different stages of convective mixing. The gas and liquid densities are presented in Fig.5.22. The gas density gradually grows with increasing hydrostatic pressure. The liquid density deviation is 18.51 kg/m 3 .

Figure 5 . 19 -

 519 Figure 5.19 -CO 2 (aq) concentration map in range [0, max (CO 2 (aq))]: max (CO 2 (aq)) = 1.296 molal at 120 d (top), max (CO 2 (aq)) = 1.294 molal at 16 y (center) and 30 y (bottom), [100, 2500] m.

Figure 5 .

 5 Figure 5.20 -HCO - 3 concentration map at 120 d (top), 16 y (center) and 30 y (bottom), [100, 2500] m. The HCO - 3 scale is [0, 0.00087] molal.

Figure 5 . 21 -

 521 Figure 5.21 -pH map at 120 d (top), 16 y (center) and 30 y (bottom), [100, 2500] m.The pH scale is[3.07, 6.63].

Figure 5 . 22 -

 522 Figure 5.22 -Gas (top) and liquid (bottom) density map at 30 y, [100, 2500] m. The gas and liquid density scales are [395, 493] and [1022, 1040] kg/m 3 , respectively.

Figure 5 . 23 -Figure 5 . 24 -

 523524 Figure 5.23 -Prediction of current radius at the top of the aquifer without convective dissolution by analytical and numerical methods. The analytical estimation also neglects diffusion and dissolution, Eq. (5.4.1); the numerical results account for diffusion and dissolution.

Fig. 5

 5 Fig. 5.25.

Figure 5 . 25 -

 525 Figure 5.25 -Graphical representation of the pilot in Tongliao, Inner Mongolia, China

Figure 5 . 27 -Figure 5 . 28 -

 527528 Figure 5.27 -Injection of {CO 2 + Air}: gas density at 10 y, [100, 1800] m. (The x-axis is scaled, 4 : 1.)

Fig. 5 .

 5 31 presents the mole fraction profiles in gas phase at 10 y. The gas convective instabilities of the Air scenario are emphasized by the fact that N 2 and O 2 constitute a thin front in contrast to the advanced CO 2 rich layer in the SO 2 case. As was mentioned in Sec. 5.3.2, this gas behaviour is associated with the convective flow in the liquid phase, Figs.5.32 and 5.33.

Figure 5 . 29 -

 529 Figure 5.29 -Ratio of the total CO 2 (g) quantity to the injected CO 2 for the injection scenarios of pure CO 2 , {CO 2 + Air} and {CO 2 + SO 2 } during 10 y.

Figure 5 . 30 -

 530 Figure 5.30 -Convective impact for the injection scenarios of pure CO 2 , {CO 2 + Air} and {CO 2 + SO 2 }: ratio of total dissolved CO 2 with and without convection.

Figure 5 . 32 -Figure 5 . 33 - 6

 5325336 Figure 5.32 -Injection of {CO 2 + Air}: liquid density at 10 y, [100, 1800] m. (The x-axis is scaled, 4 : 1.)

  AbbreviationsAIMadaptive implicit method CFL Courant-Friedrichs-Lewy number CPA-PR cubic plus association Peng-Robinson equation of state (Hajiw et al., 2015) CPR constrained pressure residual (CPR) preconditioner Cao et al. (2005) CPR Peng-Robinson equation of state modified by Coquelet et al. (2004) CSP compositional space parametrization method (Voskov and Tchelepi, 2009) DAE differential algebraic equations based method

Keywords:

  Compressible two-phase flow, Reactive transport, Operator splitting based coupling, HYTEC, Equation of state, Impure CO 2 storage, Convective mixing.

  [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF]. Numerical difficulties at the extreme values of saturation also come from the Jacobian structure, p. 21. As advantage of PS, the initial and Dirichlet boundary conditions have direct physical meanings.The system (2.2.2) on (p α , S α ) consists of one parabolic equation on the pressure (that turns to elliptic type in the incompressible case) and one degenerate parabolic equation on the saturation. If the capillary forces are neglected, the latter becomes hyperbolic, called Buckley-Leverett equation, and admits discontinuous solutionsshock waves. Since for the hyperbolic problems, several weak solutions can exist, the

numerical scheme should satisfy an entropy condition

[START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]

. In Sec. 2.5.2, we demonstrate numerically that the Godunov method provides entropy-satisfying Riemann solutions.

• Decoupled pressure and saturation equations (DPS) Decoupling eliminates the temporal derivative of saturation and leads to the pressure equation by summing the equations of system (2.2.2) divided by ρ α

  n is the outward normal, we discretize the integrals of

	Eqs. (2.3.1) and (2.3.2) over the element Ω i , ∀i ∈ I\I D α :

Table 2

 2

	.1 -Parameters of the axisymmetric problem admitting self-similar variable R/ √	t.
	Initial Conditions
	Ω = [0.3, 10 5 + 0.3] × [0, 100] [m 2 ]
	Ω: p l0 = 120 [bar]	S g0 = 0
	Boundary Conditions
	Γ N = [(0.3, 0); (0.3, 100)] [m]
	Γ D = [(10 5 + 0.3, 0); (10 5 + 0.3, 100)] [m]

  0

	2.5. VERIFICATION						25
			100
			80
		[bar]	40 60				p c p c
			20
			0	0	0.2	0.4	0.6	0.8	1
						S g
	Figure 2.2 -Capillary pressure function of S g and its derivative: p c -van Genuchten,
	Eq. (2.1.8): m = 0.457, 1/α = 19.61[kP a].
				Matrix Properties
		Porosity φ 0.12
	1 0.2 0.4 0.8 0.6	10 Fluid-Matrix Properties Permeability K 10 -13 [m 2 ] van Genuchten capillary presssure Eq. (2.1.8): m = 0.457 1/α = 19.61[kP a] van Genuchten liquid relative permeability -50 0 -10 k rl k rg -40 -20 -30	k rl k rg
	0	Eq. (2.1.11): m = 0.457 0.2 0.4 0.6 0.8 Brooks-Corey gas relative permeability S lr = 0.3, S gr = 0 0 1 -60 0 0.2 S g	0.4	S g	0.6	0.8	1
		Eq. (2.1.10): λ = 2		S lr = 0.3, S gr = 0.05
					Fluid Properties
		Density ρ l = 1000 [kg/m 3 ] ρ g = 666 [kg/m 3 ]
		Viscosity µ l = 10 -3 [P a • s] µ g = 4.8 × 10 -5 [P a • s]

  Irina Sina,⇤ , Vincent Lagneau a , Jérôme Corvisier a a MINES ParisTech, PSL Research University, Centre for geosciences and geoengineering, 35 rue Saint-Honoré, F-77305 Fontainebleau Cedex, France

Table 1 :

 1 Nomenclature

	Latin Symbols
	a	attraction parameter in Peng-Robinson EOS
	a ↵	dispersivity in phase ↵, [m]
	A s	specific surface area, [m 2 /m 3 solution] or [m 2 /kg mineral]
	b	van der Waals covolume
	c k l c k s c m g	total liquid mobile concentration of basis species k, [mol/kg w] immobile concentration of basis species k, [mol/kg w] gas concentration of basis species m, [mol/m 3 ]

Table 1

 1 Continued: Nomenclature

	C i	concentration of primary species i in chemical module
	d	dissolution parameter of transport model
	D ↵	molecular di↵usion coe cient of phase ↵, [m 2 /s]
	D e ↵	e↵ective di↵usion coe cient of phase ↵, [m 2 /s]
	e	evaporation parameter of transport model
	F	residual function
	g	gravitational acceleration vector, [m/s 2 ]
	J	Jacobian
	k	kinetic constant, [mol/m 2 /s]
	k fl max	maximum number of iterations in flow coupling
	k r↵ K	relative permeability of phase ↵ intrinsic permeability tensor, [m 2 ]
	K i	K-value/equilibrium ratio
	K j	equilibrium constant of reaction j
	K s	solubility constant of solid phase
	K h i	Henry's law constant
	M	molecular weight, [kg/mol]
	max	maximum over the modeled domain
	n ↵	quantity of matter in phase ↵
	n	normal vector
	N f	number of fluid phases ↵
	N g	number of gas species
	N c	number of basis species in chemistry module
	N p	number of phases
	N r	number of independent chemical reactions
	N s	number of species in chemistry module
	p ↵	liquid/gas pressure, [P a]
	p c	capillary pressure, [P a]
	p b	entry pressure of capillary pressure model, [P a]
	P	pressure in thermodynamics, [P a]
	q ↵	mass source term of phase ↵, [kg/s]
	q m g q k l	source term of basis species m in gas phase, [mol/m 3 ] source term of basis species k in liquid phase, [mol/kg w]
	Q	ion activity product
	R	gas constant, [J/K/mol]
	R ↵	reaction term of phase ↵, [kg/s]

Table 3 :

 3 2D problem of CO 2 injection in a fully water-saturated domain:

	time stepping
		min(dt), s max(dt), s mean(dt), s
	DUNE	156.25	5000	3579.7
	HYTEC	156	5400	4474.9

Table 4 :

 4 2D problem of CO 2 injection in a fully water-saturated domain:

	information of Newton successful and unsuccessful iteration number (Ni),
	total number of steps (successful and unsuccessful) and total execution time
		Tot. of time steps mean(Ni) Tot. of Ni Tot. exec. time, s
	DUNE	2249	3.9	-	13975
	HYTEC	1380	5.2	7183	12708

Table 1 :

 1 Nomenclature

	Latin Symbols
	a	attraction parameter in Peng-Robinson EOS
	a α	dispersivity in phase α, [m]
	A s	specific surface area, [m 2 /m 3 solution] or[m 2 /kg mineral]
	b	van der Waals covolume
	c l i	total liquid mobile concentration of basis species i, [mol/kg w]
	ci	immobile concentration of basis species i, [mol/kg w]
	c g i	gas concentration of basis species i, [mol/m 3 ]
	C i	concentration of primary species i in chemical module
	D α	molecular diffusion coefficient of phase α, [m 2 /s]
	D ef f α	effective diffusion coefficient of phase α, [m 2 /s]

Table 2

 2

	details the zone geometry and decomposition; the boundary condi-
	tions represent the average aquifer velocity 0.15 m/y at a depth 3.3 km.

Γ in,out :

Table 2 :

 2 • C and pressure 300 bar and distributed homogeneously in the Reservoir and Contact zones.{x 0 , x L ; z 0 , z H }, [m] S g p l Zone decomposition

	Global domain Ω {0, 2000; 0, 300}	hydrostatic
	Reservoir	{0, 2000; 110, 300}	0.9
	Contact	{0, 2000; 100, 110}	0.2
	Aquifer	{0, 2000; 0, 100}	0
	ω	0.2
	K, [m 2 ]	1e-13
	p c

Table 3 :

 3 Hydrodynamic parametersModerate salinity, calcite, and/or siderite are added, depending on the complex-

150

ity of the scenario.

  6. CH 4 and CO 2 loss is limited, quasi linear in time, and similar for both EOS. On the contrary, H 2 S evacuation rate is much higher (due to lower

	siderite/pyrrhotite	Siderite	1 kg/l
	calcite	NaCl Calcite	0.1 molal 10 kg/l
	gas	CO (g) H S(g) CH (g)	0.05 0.15 0.8
	Geochemical solution	species/mineral	y i /concentration

Table 4 :

 4 Geochemical parameters

	Geochemistry Fl.Dynamics & EOS	GCase 1/ GCase 2 GCase 3 FCase 2 FCase 3 FCase 4 FCase 5 FCase 6 FCase 7	FCase 1
	Part		Parameter\Case

  < 1, see Table4.1. In the supercritical region, it was fitted in regard to the solubility of hydrogen and methane in hydrocarbon liquids.

	0.511614 (0.567879) 0.401219 0.024955
	M 0.911807 (2.788270) 0.784054 (0.774270) 4.963070 1.248089
	N 1.948150 (0.161590) 2.812520 (2.575067)	-0.2	-8.0
	at T r

Table 4 . 2 -

 42 Comparison of PR78, SRK, GPR, HPR. AAD% of the vapor pressure and liquid density of Ar, CO 2 , N 2 , SO 2 , O 2 , H 2 S, H 2 O and overall 49 pure compounds, extracted from Haghtalab et al. (2011).

			vapor pressure			liquid density
		PR78 SRK GPR HPR PR78 SRK GPR HPR
	Ar	0.31 1.64	1.08 1.05 10.34	3.62 10.44	3.63
	CO 2	0.71 0.39	1.84 0.90	3.16 10.88	3.12	2.42
	N 2	0.67 1.19	0.21 0.42	9.79	3.46	9.84	3.45
	SO 2	2.39 2.08	2.87 1.94	2.33 12.98	2.37	2.72
	O 2	1.57 1.54	0.41 0.45 10.61	3.37 10.69	3.60
	H 2 S	2.12 1.17	1.83 1.19	7.24	5.22	7.23	2.26
	H 2 O	3.33 5.91	3.94 4.20 18.44 27.79 18.48 17.50
	Overall	1.6	1.7	1.5	1.4	7.2	11.7	7.2	4.8

4.1.3 α and b functions: vapor pressure and liquid molar volume

Table 4 . 3 -

 43 Comparison of PR78, SRK, VTPR, gVTPR, T-VTPR, gT-VTPR. AAD% of the liquid density of Ar, CO 2 , N 2 , CH 4 , H 2 O and overall 44 pure compounds, extracted

	from Ahlers and Gmehling (2001).
		PR78 SRK VTPR gVTPR T-VTPR gT-VTPR
	Ar	8.87	5.92	4.12	4.21	2.03	2.08
	CO 2	6.03 14.26	6.36	6.85	0.87	0.88
	N 2	8.50	7.08	4.70	4.73	2.33	2.38
	CH 4	10.09	5.14	3.08	3.63	2.30	2.37
	H 2 O	21.22 29.98	10.81	12.14	1.20	3.64
	Overall	6.7	13.1	3.0	4.1	2.0	2.5

Table 4 . 4 -

 44 Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the gas and supercritical density of CO 2 at 50 and 100 • C, 20 -350 bar.

		T = 50 • C
	P, bar \ PR78 TPR VTPR T-VTPR
	20 -300 100 -120 100 -250	2.9 6.2 4.1	2.8 5.5 3.9	4.3 7.0 6.1	3.9 3.0 3.1
		T = 100 • C
	P, bar \ PR78 TPR VTPR T-VTPR
	100 -350 100 -150 250 -350	2.7 0.5 2.7	1.8 0.7 1.6	3.2 0.5 3.8	2.5 0.5 2.7
	Table 4.5 -Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the supercritical
	density of H 2 S at 100 • C, 100 -350 bar.
	P, bar \ PR78 TPR VTPR T-VTPR
	100 -350 100 -150 250 -350	3.2 2.2 4.6	3.4 2.5 4.9	3.4 7.0 1.4	13.8 6.5 17.9

Table 4 . 6 -

 46 Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the supercritical density of CH 4 at 100 • C, 100 -350 bar.

	P, bar \ PR78 TPR VTPR T-VTPR
	100 -350 100 -150 250 -350	2.6 1.1 3.1	2.5 1.1 2.9	0.8 0.1 1.2	0.8 0.1 1.2

Table 4 . 7 -

 47 Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the supercritical N denisity at 50 • C, 100 -250 bar.

	P, bar \ PR78 TPR VTPR T-VTPR
	100 -120 100 -250	1.6 2.4	2.4 3.4	1.6 2.1	1.6 2.1

Table 4 . 8 -

 48 Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the supercritical O density at 50 • C, 100 -250 bar.

	P, bar \ PR78 TPR VTPR T-VTPR
	100 -120 100 -250	1.6 2.0	1.8 2.2	0.6 0.4	0.6 0.4
	Table 4.9 -Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the supercritical Ar
	density at 50 • C, 100 -250 bar.
	P, bar \ PR78 TPR VTPR T-VTPR
	100 -120 100 -250	1.8 2.2	1.6 2.0	0.2 0.9	0.2 0.9

Table 4 .

 4 10 -Comparison of PR78, TPR, VTPR, T-VTPR. AAD% of the liquid SO density at 50 • C, 100 -250 bar. Binary system {CO 2 + H 2 O}: (T-)VTPR vs. experimental data We present the AAD of the binary system liquid density predicted by VTPR models, because only pure compounds were considered in Ahlers and Gmehling (2001). The linear mixing rule for the volume translation c was employed. The reference 29 experimental data of {CO 2 + H 2 O} presented in Li et al. (2004) covers the pressure range 30 -290 bar at 59 • C, hence the CO 2 reduced temperature T r = 1.1. The results are given in Fig. 4.8. The VTPR model underestimates the liquid density providing 3.6% AAD while the T-VTPRoverestimation is of 2.5% AAD. The accuracy of VTPR models for the binary system is of the same magnitude as that for pure compounds Table4.3. Liquid density of the binary system {CO 2 + H 2 O} at 59 • C and 30 -290 bar predicted by VTPR, T-VTPR vs. experimental data[START_REF] Li | Densities and solubilities for binary systems of carbon dioxide+ water and carbon dioxide+ brine at 59 c and pressures to 29 mpa[END_REF].4.2.3 Water binary systems with CO 2 , H 2 S, CH 4 : CPA-PR, (T-)VTPR,

	P, bar \ PR78 TPR VTPR T-VTPR
	100 -120 100 -250	1.3 1.5	1.3 1.5	0.3 0.2	1.7 1.8

Table 4 .

 4 11 presents the comparison of the liquid density, the error is given relatively to CPA-PR results.

Table 4 .

 4 11 -Comparison of CPA-PR, VTPR, T-VTPR and GERG-2008. The error, relative to CPA-PR of the liquid density of the binary systems {CH 4 +H 2 O}, {CO 2 +H 2 O}, {H 2 S+H 2 O} at 100 • C and 300 bar, in %.

		GERG-2008 TPR VTPR (TPR) T-VTPR (TPR)
	{CH 4 +H 2 O}	1.1	17.9	3.9	3.0
	{CO 2 +H 2 O}	0.8	16.7	3.6	2.8
	{H 2 S+H 2 O}	1.3	17.4	3.9	3.6

  • C, 1 atm), calculated by ρ w = 62.368 + 0.438603S + 1.60074 × 10 -3 S 2 ,(4.3.3) where S in wt% and B w is the formation volume fraction of formation water, the temper-∆V wT = -1.0001 × 10 -2 + 1.33391 × 10 -4 T + 5.50654 × 10 -7 T 2 , (4.3.7)where T and P are in • F and in psia, respectively. The validity of model was given with 2% deviation from the experiment data for P < 350 bar, T < 127 • C and S < 30 wt%,

	ature and pressure dependent function given by
	B w = (1 + ∆V wP )(1 + ∆V wT ),	(4.3.4)
	∆V wP = -1.95301 × 10 -9 P T -1.72834 × 10 -13 P 2 T	(4.3.5)
	-3.58922 × 10 -7 P -2.25341 × 10 -10 P 2 ,	(4.3.6)

Table 4 .

 4 13 -Parameters a ij for the Altunin viscosity model of CO 2 , Eq.(4.3.8).

	• C

  CO 2 is the molality of CO 2 . The {H 2 O + NaCl} viscosity formulation will be presented below, Eq.(4.3.15). Figure 4.12 -Viscosity of H 2 O and H 2 O + CO 2 with x CO 2 = 0.02 at 45, 50, 75 and 100 • C and 20 -300 bar predicted by Islam-Carlson models Eqs. (4.3.11) and (4.3.12).

	1.0134 CO 2 200 40 • C: H 2 O+CO 2 , 150 P, [bar] H 2 O 50 • C: H 2 O+CO 2 H 2 O 75 • C: H 2 O+CO 2 H 2 O 100 • C: H 2 O+CO 2 H 2 O The Mao and Duan (2009) model predicts the viscosity with the accuracy about 1% for (4.3.14) 143 2.5e-04 3.0e-04 3.5e-04 4.0e-04 4.5e-04 5.0e-04 5.5e-04 6.0e-04 6.5e-04 7.0e-04 50 100 250 300 Viscosity, [Pa • s] Modeling {H 2 O + NaCl} where x 4.3. FLUID PROPERTIES binary systems including {H 2 O + NaCl} up to 6 mol/kg at the T P range 273 -623 K and 1 -1000 bar:

Table 4 .

 4 14 -Parameters a i , b i , c i , d i for the Islam and Carlson viscosity model of H 2 O, Eq. (4.3.11).

Table 4 .

 4 15 -Parameters a r i , b r i for the Islam and Carlson viscosity model of {H 2 O+CO 2 }, Eq. (4.3.13).

Table 4 .

 4 16 -Parameters A i , B i , C i for the Mao and Duan viscosity model of {H 2 O +

  ).

	154 CHAPTER 5. CCS, CO-INJECTED IMPURITIES AND CONVECTIVE MIXING
	Concentration, [molal]	0 0.2 0.4 0.6 0.8 1 1.2	T = 1 min CO 2 N 2				Concentration, [molal]	0 0.2 0.4 0.6 0.8 1 1.2	T = 1 min CO 2 SO 2
		0	0.02	0.04	0.06	0.08	0.1	0	0.02	0.04	0.06	0.08	0.1
				x, [m]						x, [m]
	Concentration, [molal]	0 0.2 0.4 0.6 0.8 1 1.2	T = 55 min CO 2 N 2				Concentration, [molal]	0 0.2 0.4 0.6 0.8 1 1.2	T = 55 min CO 2 SO 2
		0	0.02	0.04	0.06	0.08	0.1	0	0.02	0.04	0.06	0.08	0.1
				x, [m]						x, [m]
	Concentration, [molal]	0 0.2 0.4 0.6 0.8 1 1.2	T = 1 d CO 2 N 2				Concentration, [molal]	0 0.2 0.4 0.6 0.8 1 1.2	T = 10 d CO 2 SO 2
		0	0.02	0.04	0.06	0.08	0.1	0	0.02	0.04	0.06	0.08	0.1
				x, [m]						x, [m]

  .2.2)The CO 2 concentration is lower in the presence of Ar, O 2 , N 2 than that of the pure CO 2 and SO 2 impurity cases, Figs.5.7 and 5.8. At equilibrium state, the concentration of CO 2 is the highest in the {CO 2 + SO 2 } system. However, that does not involve a better mass solubility capacity, as the initial volume is identical for the models but not the initial mass of CO 2 . Figure 5.7 -Gas mole fractions of impurity gases and pure CO 2 (left top), CO 2 concentration in molal (right top), gas saturation (left bottom) and the gas density in kg/m 3 (right bottom) of CO 2 , {CO 2 + N 2 }, {CO 2 + Ar}, {CO 2 + O 2 } and {CO 2 + SO 2 } models at T = 55 min.

	156 CHAPTER 5. CCS, CO-INJECTED IMPURITIES AND CONVECTIVE MIXING 5.2. GAS DYNAMICS AND CHROMATOGRAPHIC PARTITIONING	157
											T = 1 min
							1.4
	Gas mole fraction	0 0.2 0.4 0.6 0.8 1		CO 2 + N 2 + Ar + O 2 + SO 2		2 (aq), [molal] CO	0.2 0.4 0.6 0.8 1 1.2 0 1.4	0	CO 2 + N 2 + Ar + O 2 + SO 2	0.02	0 0.2 0.4 0.6 0.8 1 1.2 0.04 CO x, [m] CO 2 + N 2 + Ar 0.06 2 (aq), [molal] + O 2 T = 10 min + SO 2	0.08	0.1
	Gas saturation	0 0.02 0.04 0.06 0.08 0.1	0 0	CO 2 + N 2 + Ar + O 2 + SO 2	0.02 0.02	0.04 x, [m] 0.04 CO 2 (aq), [molal]	0.06 0.06 0 0.2 0.4 0.6 0.8 1 1.2	0	0.08 0.08	0.1 0.1 0.02	0 0 x, [m] 0 100 200 300 400 CO 2 + N 2 + Ar + O 2 + SO 2 0.04 ] 0.06 0.02 0.02 Gas density, [kg/m 3	0.04 x, [m] 0.04 0.08	0.06 0.06 0.1	0.08 0.08	0.1 0.1
						x, [m]							x, [m]
											T = 10 d
							1.2
						CO 2 (aq), [molal]	1
							0.8
								0		0.02	0.04	0.06	0.08	0.1
											x, [m]

Table 5 . 2 -

 52 Parameters of the 2D problem of the convective mixing.

Á. Sáinz, personal communication,

A. Michel, T. Faney, personal communication, 2015

The abbreviations are detailed in Table

C.Coquelet, M. Hajiw, personal communication,

2015: for instance, the AAD of water density is 2.7% at [273.15, 643.15] K and < 1% at [273.15, 573.15] K.

Remerciements I am grateful to École des Mines and BRGM for the financial support that made this work possible. I would also like to thank the MoMaS, GRC and SIAM organizations for their conferences and workshops. Je remercie mon directeur de thèse, Vincent Lagneau, pour chacune des rencontres, éclairantes, chacune des explications, patientes et passionnantes, pour l'esprit scientifique qu'il m'a transmis, pour la liberté confiée et le temps consacré. Merci de m'avoir fait poser les problèmes et chercher de nouvelles possibilités. Спасибо за его и мои pérégrinations, заботу и за то, что был рядом в важные моменты моей жизни. Je tiens à remercier Jérôme Corvisier d'avoir été mon maître de thèse. Je lui suis reconnaissante pour son soutien et pour m'avoir fait découvrir des visions alternatives sur des problèmes critiques. I thank Denis Voskov for reading my thesis and providing valuable recommendations. I am grateful for his perceptive questions and fruitful discussions over the last two years. Я глубоко признательна Денису за советы, знания и за душевное общение. I express my profound gratitude to Ulrich Mayer for reviewing the thesis, for his outstanding editorial contribution to this work, insightful questions and excellent discussion during the oral defense that disclosed new interesting issues to me. Je tiens à remercier Roland Masson d'avoir accepté de présider mon jury de thèse. Je lui suis très reconnaissante pour toutes les connaissances qu'il m'a apportées depuis notre toute première rencontre, ce qui a joué un rôle important dans mes recherches. Je tiens à remercier Mohamed Azaroual et Pascal Audigane d'avoir été mes encadrants. Je leur suis fortement reconnaissante pour l'accueil chaleureux au BRGM, pour les conseils, les échanges et l'expertise qui ont favorisé l'avancement de ce travail. Je remercie vivement Christophe Coquelet pour la rédaction de la partie thermodynai Acknowledgements This work was supported by MINES ParisTech and BRGM. I. Sin thanks Christophe Coquelet (Centre Thermodynamique des Procds, MINES Paris-Tech) for his perspicacious suggestions on the thermodynamics development and Claude Tadonki (Centre de Recherche en Informatique, MINES Paris-Tech) for the constructive discussions on the linear solvers.

Hence, the gases on the RHS to CO 2 of the inequality (from Ar to N 2 ) should evolve the similar profile of the normalized mole fraction with high and early peaks at the front, and vice versa, SO 2 and H 2 S yield the lower and later breakthrough. [START_REF] Bachu | Chromatographic partitioning of impurities contained in a CO 2 stream injected into a deep saline aquifer: Part 1. effects of gas composition and in situ conditions[END_REF] demonstrated experimentally this effect in the case of H 2 S, N 2 , SO 2 , CH 4 and later, numerically for H 2 S admixture Bachu et al. (2009). To complete sequence (5.2.1), we propose the simulation of the following problem where the resident gas is represented by pure CO 2 and then by CO 2 + 5% of impurities: SO 2 , Ar, O 2 N 2 .

Problem description

The homogeneous medium is 10 cm long, and 1 cm high. Its right half is fully saturated by water. The left one contains some amount of gas (S g = 0.2), the water is at equilibrium with the gas phase. To exclude the gravity forces, the geometry is simplified to 1D, also the liquid density is supposed to be constant. Supposing the P T conditions of a deep reservoir, the temperature is 50 • C and the initial pressure is 100 bar. The parameters of the problem are given in Table 5.1. We firstly present the results for the case when the gas consists only of CO 2 . Then, the impact of impurities is observed on the example of SO 2 and N 2 , as they are the most and the least soluble gases of sequence 5.2.1. Finally, all five systems are compared: CO 2 , +SO 2 , +Ar, +O 2 , +N 2 .

CO 2

Driven by the capillary forces, the gas, CO 2 in this case, moves forward from the left to the right, regressively occupying the porous medium and dissolving in the fresh water until reaching its stationary state, Fig. 5.3 (left). During the continuous gas dissolution, both pressures, liquid and gas, decrease. In Fig. 5.3 (right), the temporal evolution of the phase pressures is plotted at x = 0.005, at the left border of the domain. The shape of the pressure function is correlated with the decreasing dissolution rate as the gas approaches slower to the right border.

CO 2 + 5% impurities

Supposing that the gas is represented by 95 mol% of CO 2 and 5 mol% of impurities, we compare the numerical results for N 2 and SO 2 (SO 2 is more soluble than CO 2 and N 2 is less soluble than CO 2 ). The chromatographic partitioning of gas composition occurs arising from the different solubility rates, as expected, Fig. Nevertheless the mole fraction normalized to its initial value (5 mol%) is higher during the simulation time, until the stationary state is reached. The gas partitioning of CO 2 and SO 2 is opposite to the system {CO 2 + N 2 }, CO 2 takes a leading position of the gas front, The partitioning of the aqueous components can be found in addition to that of the gas phase. Being less soluble than CO 2 , N 2 forms a slight accumulation ahead towards the displacement, which rapidly diminishes. Meanwhile, the high solubility of SO 2 yields the pronounced concentration behind the moving front. Fig. 5.5 shows the concentration contours of aqueous components for both systems.

Bachu and Bennion reported the unexpected earlier breakthrough of CO 2 in the cases of N 2 and CH 4 suggesting that these gases pull CO 2 along. Based on the present results, we confirm the associated movement of CO 2 along with the impurities. Fig. 5.6 (right column) demonstrates the advancing front of gas saturation of the system with less soluble N 2 and consequently, a delay of gas saturation profile of the system CO 2 + SO 2 . The different position of the fronts is also conditioned by the gas density. Fig. 5.6 (left column) illustrates that the admixtures result in a high deviation of gas density: SO 2 increases the density, whereas N 2 reduces it. The gas chromatographic partitioning implies a density heterogeneity. In this problem, such important changes in gas mixture density determine the dynamics of the system and, in particular, for the gas saturation front. The models of CO 2 and {CO 2 + N 2 } achieve the equilibrium state within one day, while it takes about

Convective mixing

In the context of CO 2 geological sequestration, the dissolution is identified as an important long-term trapping mechanism. As the supercritical CO 2 is buoyant to the ambient brine, the current spreads along the reservoir or aquifer inducing a security issue of the CO 2 where H is the height, ∆ρ is the density increase due to the dissolved CO 2 . In Eq. 5. (5.3.3) where X CO 2 is the mass fraction of CO 2 (aq), Sh is the Sherwood number, α = 0.19 ± 0.01.

Applying linear stability analysis to the convective problem in the single phase only with diffusion, [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF] derived the wavelength λ c and the onset time of instability t c .

Then [START_REF] Elenius | Effects of a capillary transition zone on the stability of a diffusive boundary layer[END_REF][START_REF] Elenius | Convective mixing influenced by the capillary transition zone[END_REF]) considered a capillary transition zone that introduced the vertical flux at the instability interface and accelerated the convection: .3.4) where the left limit includes the capillary forces, the right limit is given for the diffusive dominant mode, [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF]. The onset time is expressed as follows

( 5.3.5) where c 0 is the constant: In the presence of the gas layer, the capillary forces create additional perturbations. The development of the density driven fingering can entail the uneven contact surface between the gas and liquid phases. The subject is raised in this section by 2D numerical modeling.

Problem description

We model the convective mixing problem with a thin CO 2 (g) layer at the top of the 2D rectangle domain supplying CO 2 (aq) at the upper boundary by pure diffusion process.

Table 5.2 gives the problem parameters. The fluid, rock properties are identical to those of the previous modeling, Table 5.

Numerical results

The initial conditions for saturation are out of equilibrium. The capillary forces set the initially present gas in motion so that the gas gradually spreads downwards, developing the capillary fringe and increasing the dissolution. The system passes through the equilibration stage within the first 39 d (t eq ) as the gravity and the capillary pressure are taken -The higher the convective rate, the faster the aquifer water under the plume becomes saturated by the dissolved gases that decelerates the dissolution rate afterwards and hence extends the plume lifetime if the aquifer height is limited, [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Macminn | Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers[END_REF].

The viscosity ratios also affect on the current dynamics. We point out some factors, although the viscosity is supposed to be constant in this study. The spreading increases with M ph = µ f w /µ inj , where µ f w is the viscosity of the aquifer water, µ inj is the viscosity of the injected gas. The viscosity ratio M l = µ f w /µ diss also affects the current dynamics, µ diss is the viscosity of the aquifer water with the solutes. MacMinn et al. ( 2012) showed analytically and experimentally that the current extension decreased with M l lengthening its life-time.

Problem description

Taking into account the chemical and thermodynamic properties of gas impurities, we model the co-injections of SO 2 and Air that have contrary characteristics of solubility, density, chromatographic partitioning effect, Sec. 5.2. Three gas injection scenarios are proposed • pure CO 2 ;

• 95 mol% of CO 2 with 5 mol% of SO 2 ;

• 95 mol% of CO 2 with 5 mol% of Air (i.e., N 2 and O 2 , 4 and 1 mol% of the total gas composition, respectively).

The mass injection rate is identical, 10 kg/s and performed through the whole height of the aquifer during 10 y. The problem parameters and geometry are identical to those of the pure CO 2 injection model, Sec. 5.4.1. The coarse grid is used to perceive the general dynamics of the system and its sensitivity to the injected gas composition. The cell length is variable, from 4 m to 1 km.

Numerical results

Given the equal mass injection rate and different gas composition, the gas volumetric rate q vol of three scenarios can be written by following q vol,SO 2 < q vol,CO 2 < q vol,Air , (5.5.1) as their densities are in the inverse order. As a results, the larger gas extension (maximum radius) is observed in the case of Air admixture and the smaller radius in the case of SO 2 admixture. Fig. 5.26 displays the maximum gas radius of the three cases. We stress that the order of plume extensions is the consequence of the gas composition density. The gas current density in the case of Air admixture is low in contrast to that of SO 2 . In addition, the front of the plume primarily consists of N 2 and O 2 due to the chromatographic partitioning. On the contrary, the co-injection of SO 2 entails the higher gas density and CO 2 rich gas front of the plume, Figs. 5.27 and 5.28. Since the gas current of the Air case spreads faster than those of pure CO 2 and SO 2 cases, the CO 2 dissolution rate is the highest. Fig. 5.29 illustrates the ratio of CO 2 (g) quantity to its total injection mass during 10 y. Similarly to Sec. 5.4.2, the curves decline and converge with time, indicating a growing impact of convection on the dissolution, which is more important in the case of SO 2 .