Improvement of the Calculation Accuracy of Acid Gas Solubility in Deep Reservoir Brines: Application to the Geological Storage of CO 2

Résumé — Amélioration de la précision du calcul de la solubilité des gaz acides dans les saumures des réservoirs profonds : application au stockage géologique du CO 2 — L’évaluation des conséquences à court et long terme de l’injection de CO 2 dans les aquifères nécessite à la fois des expérimentations en laboratoire et l’utilisation de modèles afin de mieux comprendre les différents processus physicochimiques induits. La modélisation de l’injection dans un réservoir où les conditions in situ sont souvent relativement élevées en termes de température (au-delà de 50 °C), pression (plusieurs centaines de bars) et salinité (supérieure à celle de l’eau de mer), implique donc d’avoir à notre disposition des outils numériques capables de prendre en compte à la fois les effets spécifiques de chaque électrolyte et la non-idéalité de la phase CO 2 gaz. Cet article propose une évaluation de la cohérence des différentes corrections (activité, fugacité, influence de la pression sur les constantes thermodynamiques) à considérer dans les modèles géochimiques afin de répondre aux besoins en termes de précision de calcul. Ces corrections ont été introduites dans le logiciel de modélisation thermo-cinétique Abstract — Improvement of the Calculation Accuracy of Acid Gas Solubility in Deep Reservoir Brines: Application to the Geological Storage of CO 2 — The assessment of the short and long term consequences of CO 2 injection in aquifers requires both laboratory experiments and numerical modelling in order to better understand the various physical-chemical processes taking place. Modelling injection in a reservoir, where relatively high temperature (above 50°C), high pressure (several hundreds of bars), and high salinity (greater than that of seawater) conditions are likely to be encountered, thus requires numerical tools able to take into account the specific effects of the various electrolytes dissolved in brines, and the non-ideal behaviour of the CO 2 gaseous phase. This study evaluates the consistency of the various corrections (activity, fugacity, influence of pressure on thermodynamic constants) to be taken into account in geochemical models to meet these calculation accuracy requirements. These corrections were implemented in the thermo-kinetic modelling software SCALE2000 (Azaroual et al. , 2004a) which was used to check their validity by comparing the calculation results with available experimental observations and other results from CO 2 solubility calculation models. An estimation of


Abstract -Improvement of the Calculation Accuracy of Acid Gas Solubility in Deep Reservoir
Brines: Application to the Geological Storage of CO 2 -The assessment of the short and long term consequences of CO 2 injection in aquifers requires both laboratory experiments and numerical modelling in order to better understand the various physical-chemical processes taking place. Modelling injection in a reservoir, where relatively high temperature (above 50°C), high pressure (several hundreds of bars), and high salinity (greater than that of seawater) conditions are likely to be encountered, thus requires numerical tools able to take into account the specific effects of the various electrolytes dissolved in brines, and the non-ideal behaviour of the CO 2 gaseous phase. This study evaluates the consistency of the various corrections (activity, fugacity, influence of pressure on thermodynamic constants) to be taken into account in geochemical models to meet these calculation accuracy requirements. These corrections were implemented in the thermo-kinetic modelling software SCALE2000 (Azaroual et al., 2004a) which was used to check their validity by comparing the calculation results with available experimental observations and other results from CO 2 solubility calculation models. An estimation of the relative weight of each of the corrections for a 237 g·l -1 brine (60°C, pCO 2 = 200 bar) showed a systematic overestimation (higher than 100%) of CO 2 solubility when either salinity (NaCl equivalent) is neglected or gas is considered ideal. The error induced by the NaClequivalent approximation compared to real brine is lower (less than 5%). The second part of this study presents an application example of a hypothetical scenario of massive CO 2 injection in a carbonated reservoir; data used for the brine composition are actual data (Moldovanyi and Walter, 1992) from the Smackover site (Arkansas, United States). The simulations performed considering a representative elementary volume of saturated bulk rock (porous mineral assemblage saturated with the Smackover brine) with a prescribed constant CO 2 pressure of 150 bar, show two distinctively different behaviours whether the system is assumed to be a closed (batch reactor) or an open reactor fed by a constant brine flow rate. In the first case, the calculations performed with SCALE2000 lead to negligible variations in the mineralogy. In the second case, more representative of the dynamical nature of an injection system, the results show major modifications in the mineralogy finally leading to a strong increase in porosity (from 20% initially to 85% after 50 y of simulated time). Further calculations were carried out with SCALE2000, now considering a 1D system constituted of a set of four homogeneous identical reactors connected in series (fluid velocity of 1 m.day -1 ). With initial and boundary conditions similar to those considered earlier, and prescribing a constant pCO 2 in the first reactor only, the results showed that significant dolomite precipitation occurred in the most-downstream reactor hence inducing some CO 2 precipitation. Mass balance calculations performed on the four reactors system finally demonstrated a global loss in total mineral carbon with respect to the amounts initially available. However, the evolution trends observed in the most-downstream two reactors indicated that possible trapping might be expected beyond the relatively limited geometrical boundaries considered in the modelled system. activity of the species i in the aqueous phase, dimensionless a i,gas activity of the species i in the gaseous phase, dimensionless f i fugacity of gas i, Pa fr fraction of the total surface area effectively reactive, dimensionless I ionic strength, mol·(kg H 2 O) -1 k kinetic constant, mol·m -2 ·s -1 K thermodynamic equilibrium constant, dimensionless m i,aq molality of the aqueous species i, mol·(kg H 2 O) -1 P total pressure, Pa p i partial pressure of gas i, Pa pCO 2 CO 2 partial pressure, Pa Q ionic activity product, dimensionless R universal gas constant, J·K -1 ·mol -1 S total surface area, m 2 T temperature, K TDS Total Dissolved Solids, kg·m -3 V net effective reaction rate, mol·s -1 x i,gas molar fraction of the species i in the gaseous phase, dimensionless partial molal compressibility variation for species involved in a reaction r, m 3 ·mol -1 ·Pa -1 partial molal volume variation for species involved in a reaction r, m 3 ·mol -1 Φ i,gas fugacity coefficient for the species i, dimensionless γ i,aq activity coefficient of the species i in the aqueous phase, dimensionless.

INTRODUCTION
CO 2 sequestration in deep saline aquifers emerges as a viable option to mitigate greenhouse gas emissions. Saline aquifers, oil and geothermal reservoirs are the most interesting targets for sequestration. Results of several experimental (e.g. Tackenouchi and Kennedy, 1964;Rumpf and Maurer, 1993;Kaszuba et al., 2003;Azaroual et al., 2004b) or modelling (e.g. Gunter et al., 1997;Czernichowski et al., 2001;Xu et al., 2004;Gauss et al., 2005) studies have been published. However, knowledge on reactivity and geochemical evolution of such CO 2 enriched systems still has to be improved. Hence, the assessment of the short and long term consequences of injection requires further laboratory experiments and numerical modelling work in order to understand the key processes taking place when CO 2 is injected in a geological reservoir. High temperature (50°C and above), high pressure (several hundreds of bar), and high salinity (greater than that of seawater) conditions in which CO 2 is likely to be injected, make numerical investigations intricate as, in particular, it is necessary to take into account the non ideal behaviour of both the mixed electrolytes dissolved in brines and the CO 2 gaseous phase.
This study aims at evaluating the impact of the various corrections for aqueous and gaseous phases (activity, fugacity, influence of pressure on thermodynamic constants) to be considered in geochemical models when attempting to calculate CO 2 solubility accurately in brines at high pressure and high temperature. These corrections were implemented in the thermokinetic modelling software SCALE2000 (Azaroual et al., 2004a). After a brief presentation of the theoretical aspects of the Pitzer's formalism, on which these corrections mainly rely, the first part of this paper is devoted to a comparison between modelling results obtained with SCALE2000 and PHREEQC (Parkhurst and Appelo, 1999), and available experimental data from the literature. On the basis of this comparison, we finally propose a quantification of the relative error induced by neglecting each of the suggested corrections. The second part of this article describes an application example of this new version of SCALE2000 to a hypothetical scenario of CO 2 injection in a saline carbonated aquifer (based on actual brine composition data from the Smackover formation, United States). This example focuses particularly on the calculation of CO 2 solubility in the aqueous phase and on the evolution of mineral assemblage under the effect of high CO 2 partial pressure (150 bar).

Gas-Brine Thermodynamic Equilibrium
A reliable assessment of gas solubility has to take into account the complex thermodynamics of the gas-brine-minerals system. Gas solubility is dependent on temperature, pressure, and detailed chemical composition of both the gas and the brine. Considering a species i likely to exist in both the gaseous and the aqueous phases, its activity in the aqueous phase (a i,aq ) is related to its activity in the gaseous phase (a i,gas ) according to the general law of mass action equation: (1) where K i is the thermodynamic equilibrium constant of the reaction i aq ↔ i gas , P is the pressure (bar), and T is the absolute temperature (K). Considering a mixture of ideal gases, the activity of i is defined as: where p i is the partial pressure of i and p 0 is a reference pressure chosen to be equal to 1 bar. In the case of real gases, fugacity (f i in bar) has to be considered instead of partial pressure, so that activity can be written as: ( 3) Fugacity is a function depending on the pressure, the temperature, and the molar fraction of the considered species i in the gaseous phase, according to the following equation: where Φ i,gas , x i,gas , and x j≠i,gas are, respectively, the fugacity coefficient (dimensionless), the molar fraction (dimensionless) of the species i in the gaseous phase, and the molar fraction of each of the other species j present in the gaseous phase.
Similarly, the activity of an aqueous species depends on the pressure, the temperature, and the solution composition. Although the solution composition effect can be treated as a whole for solutions having a salinity lower than that of sea water (≈ 35 g·l -1 ), the brines that are typically targeted for acid gas sequestration generally require a differentiation of the specific effect of individual chemical components according to the general equation: where γ i,aq , m i,aq , and m j≠i,aq are, respectively, the activity coefficient (dimensionless), the molality (mol·(kg H 2 O) -1 ) of the aqueous species i in the aqueous phase, and the molality of each of the other aqueous species j; m 0 is a reference molality chosen to be equal to 1 mol·(kg H 2 O) -1 .
Combining Equations (1, 3, 4 and 5) leads to the following general formulation of gas solubility in a complex brine: When pressure and concentration are expressed in bar and in mol·(kg H 2 O) -1 respectively, Equation (6) can finally be simplified as: (7)

Activity Correction
Estimating the activity of aqueous species requires the calculation of their activity coefficient integrating the effect of salinity, temperature, and pressure. For low salinity solutions, activity coefficients can be easily calculated by ion association approaches such as Debye-Hückel's, B-Dot, or Davies'. These approaches, implemented in most of the geochemical codes such as PHREEQC (Parkhurst and Appelo, 1999) or EQ3/6 (Wolery, 1992), are mathematically simple and well parameterised (Helgeson et al., 1981).
However, for higher salinities, the calculated results diverge from observations, even in the case of NaCl solutions (see Fig. 1). In these cases, the use of a more specific interaction model such as Pitzer's model (Pitzer, 1973) is necessary. Pitzer's model, implemented in SCALE2000 (Azaroual et al., 2004a), has been specifically designed to reproduce experimental data at high ionic strength. Mean activity coefficient for NaCl (γNaCl) vs. ionic strength of the solution, calculated at 25°C, 1 bar, using various activity models.
Compared to ion association models, Pitzer's formalism is mathematically more complex and relies on many parameters. Despite the lack of some of experimental data, notably for complex electrolytic solutions at high temperature, it was finally possible to build a consistent Pitzer's database based on a very systematic and refined interpolation, extrapolation, and verification work using the simple electrolytes data available from the literature. Hence, in its actual version, the model accurately describes the behaviour of carbonate-gaswater systems for temperature and pressure conditions ranging respectively from 25 to 250°C and 1 to 1000 bar, and for high salinity solutions (ionic strength up to about 6 M) in the system Na-K-Ca-Mg-Sr-Ba-Fe-Cl-SO 4 -HS-OH-H-Ac-CO 3 -HCO 3 -H 2 S-CO 2 -AcH-SiO 2 -H 2 O. The geochemical software SCALE2000 (Azaroual et al., 2004a) relies on this Pitzer's formulation and relevant databases. SCALE2000 was validated against numerous experimental data available in the literature (Azaroual et al., 2004c). strength I (molal) for a pure NaCl solution at 25°C and 1 bar. At a salinity beyond that of seawater (ionic strength of approximately 0.7), this figure clearly shows significant discrepancies between the results obtained with SCALE2000, relying on Pitzer's formulation, and those based on simple ion association approaches (Debye-Hückel, Davies, B-Dot). Considering that Pitzer's formulation is the most adequate for problems involving high salinity solutions (Harvie et al., 1984) as confirmed by experimental measurements, it can be concluded that Davies' formalism is not reliable for thermodynamic predictions for solutions more saline than seawater. Figure 1 also confirms that the Debye-Hückel formalism can only be applied at much lower ionic strengths than the Davies formalism, which practically makes it unusable for brackish and saline natural systems.

Thermodynamic Equilibrium Constant
The accuracy of the thermodynamic equilibrium constants is of course of major importance to estimate an actual gas or mineral solubility in a complex system. Geochemical codes (EQ3/6, PHREEQC, etc.) generally only take into account the influence of the temperature on equilibrium constants, thereby neglecting the pressure effect, which is usually assumed to be less significant (Kharaka et al., 1988). However, in CO 2 geological sequestration problems, the pressure involved may reach several hundreds of bars, and in this case, its effect on equilibrium constants becomes significant. A first formulation proposed by Monnin (1990) for calculating the influence of the pressure on equilibrium constants (valid for pressures lower than 1000 bar) can be used: where T, P, and R are, respectively, the temperature, the pressure, and the universal gas constant; index 0 refers to a reference state, and and stand respectively for partial molal volume and partial molal compressibility variations for the components involved in a reaction r. This approach assumes that and do not depend on pressure but it has been shown that these functions might however vary significantly with pressure (Johnson et al., 1992). A more powerful formalism with a more extended validity domain (0 < T < 1000°C; 1 < P < 5000 bar), explicitly taking into account the pressure effect on and was developed by Johnson et al. (1992). In order to improve the accuracy of the model results, this latter formalism was finally implemented in SCALE2000. Figure 2 illustrates this pressure effect for the reaction CO 2,aq ↔ CO 2,gas . It is interesting to notice that the slope of the curves increases with temperature. Thus, considering the temperature domain relevant for CO 2 geological sequestration, Influence of pressure on the thermodynamic equilibrium constant of the reaction CO 2,aq ↔ CO 2,gas . this result highlights the increased inaccuracy when this pressure effect is not accounted for in the models.

Fugacity Correction
The calculation of gas fugacity necessitates first to evaluate the molal volume which depends on the gas phase composition, the pressure, and the temperature. Among the various equations of state (EOS) available in the literature, those from Duan et al. (1992aDuan et al. ( , 1992bDuan et al. ( , 2003 appeared to be well adapted to the problem of acid gas sequestration as they asses phase changes occurring in the pressure-temperature domain considered. The general form of these equations, derived from those of Lee and Kesler (1975), can be found in Duan et al. (1992b, Equations (1, 9); Table 1 and 2). They permit to calculate the molar volume of a pure or mixed gas phase as well as the fugacity coefficient of each gas.
Duan et al.'s formalism was implemented in version 3.1 of SCALE2000. Figure 3 presents a 2D graph of the calculated CO 2 fugacity coefficient as a function of pressure and temperature. Also represented on Figure 3 is the critical point of CO 2 . It can be observed that in gaseous and supercritical phases, CO 2 fugacity coefficient ranges from 0.3 to 1. For example, under conditions typical for CO 2 sequestration, i.e. 200 bar and 60°C, the fugacity coefficient for pure CO 2 will be 0.45, leading to a fugacity value of 90 bar instead of 200 bar of assumed partial pressure of CO 2 without applying any correction (thus assuming a fugacity coefficient of 1). Figure 4 shows a good agreement between SCALE2000 calculation results and experimental measurements of CO 2 solubility in pure water for a pressure up to 400 bar and a temperature up to 250°C, confirming the validity of the fugacity coefficient and equilibrium constant calculations.
Chemical composition of the reference waters considered in the simulations. A: initial analytical data from the Smackover formation water (fluid sample no. 71 in Moldovanyi and Walter, 1992). P1: initial speciation calculated with PHREEQC at 25°C. S1: initial speciation calculated with SCALE2000 at 25°C. P2: initial speciation calculated with PHREEQC at the reservoir temperature (71°C). S2: initial speciation calculated with SCALE2000 at the reservoir temperature (71°C). S3: speciation calculated with SCALE2000 at the reservoir temperature (71°C) assuming equilibrium with calcite. S4: speciation calculated with SCALE2000 at the reservoir temperature (71°C) assuming equilibrium with the hypothesised mineralogical composition of the reservoir (anhydrite, calcite, dolomite, quartz, siderite). P3: speciation calculated with PHREEQC at the reservoir temperature (71°C) assuming equilibrium with pCO 2 = 150 bar. S5: speciation calculated with SCALE2000 at the reservoir temperature (71°C) assuming equilibrium with pCO 2 = 150 bar (calculated fugacity = 79.8 bar). P4: speciation calculated with PHREEQC at the reservoir temperature (71°C) assuming equilibrium with pCO 2 = 79.8 bar (CO 2 fugacity value calculated by SCALE2000). * Measured at the reservoir temperature (71°C). ** Chemical element not considered for speciation calculation but taken into account for density calculation. Figure 5 also shows a satisfactory agreement between calculated CO 2 solubility and experimental measurements in NaCl, CaCl 2 , and MgCl 2 solutions for electrolyte concentrations as high as 5 molal at 1 bar and 25°C, corroborating the validity of the activity coefficient calculation at standard conditions. Another important aspect pointed out by the results of Figure 5 is the great variability of CO 2 solubility according to the nature of the electrolyte. Hence, the detailed knowledge of the chemical composition of the brine appears to be a key point in any evaluation of the CO 2 sequestration capacity of a potential site.

A
In order to evaluate the accuracy of the model at high temperature, we performed a series of calculations for a 1 M NaCl solution at 150°C for pressures ranging from 50 to 400 bar. The calculation results presented on Figure 6 can be favourably compared to the measurements of Takenouchi and Kennedy (1964). On the basis of the experimental data from Malinin and Kurovskaya (1975), we carried out other CO 2 solubility calculations in 0.44 M and 0.88 M CaCl 2 solutions for moderate temperatures (25 to 75°C). The calculation results are also in very good agreement with the measurements (Fig. 7). Pure CO 2 fugacity coefficient calculated with SCALE2000 as a function of the pressure and the temperature.      Major parameters relative to mineral properties used in batch and open-reactor simulations. Initial concentrations were calculated assuming a water-rock ratio of 1:4 in volume (20% porosity). Last five columns refer to parameters used in Equations (9) and (10). Specific surfaces were estimated based on a classical spherical grain approach. Because of the scarcity of CO 2 solubility data for electrolyte systems other than NaCl, especially at high temperatures, we were not able to do any further validation analysis. However, the results achieved up to now lead us to believe that all the corrections implemented in SCALE2000 notably improve CO 2 solubility prediction.

Initial
Figures 4 to 7 finally highlight the importance of taking into account the effects of the temperature, the pressure, and the fugacity as well as the necessity to consider the specificity of the electrolyte when calculating CO 2 solubility, particularly in the context of geological sequestration where relatively high pressure (several hundreds of bar), temperature (above 50°C) and salinity (higher than that of seawater) conditions are likely to be encountered.

Relative Weight of the Various Corrections
As an illustration of the relative weight of the different corrections accounted for in SCALE2000 (i.e. fugacity, effect of the pressure on equilibrium constants, activity model), we performed CO 2 solubility calculations using PHREEQC (Parkhurst and Appelo, 1999) without any fugacity nor activity correction, PHREEQC with the fugacity coefficient corrected "by hand" (on the basis of the value calculated by SCALE2000), and SCALE2000 taking into account all the corrections. These calculations were applied to a 1 M NaCl solution at 150°C and for pressure varying from 50 to 400 bar. The PHREEQC input file used for these calculations is provided in Appendix 1.
The graphs represented on Figure 8 show a good agreement between the results of the calculations performed with SCALE2000 and the experimental measurements over the whole pressure range. In contrast, important discrepancies between experimental and PHREEQC results (11% at 50 bar, 127% at 400 bar) are observed. The fugacity correction notably reduces the error but the discrepancy still remains noticeable (4% at 50 bar, 47% at 400 bar). Consequently, these results prove that the main error arises from neglecting the fugacity correction. However, neglecting the pressure effect on the equilibrium constants as well as considering an inadequate activity model also contribute significantly to the inaccuracy of the CO 2 solubility calculations.
In an attempt to quantify the impact of the aqueous activity model on the estimated CO 2 solubility, we performed the following series of calculations: on one hand we considered the detailed composition of the brines and, on the other hand, an equivalent-NaCl solution expressed as total dissolved solids (TDS). The equivalent-NaCl solution is a classical approach used in many solubility calculations aiming to simulate the salting out effect (Drummond, 1981). The simulation conditions were 200 bar and 60°C and the 35 brine compositions considered were selected among data from Hitchon et al. (1971), Langmuir and Melchior (1985), Kelly et al. (1986), Leach et al. (1991), Land (1995), Stueber et al. (1998), Demir andSeyler (1999), and Martel et al. (2001). The results presented on Figure 9a globally show relatively small gaps between the two calculation methods even though a significant gap can be observed for the highest salinity. The same results, plotted in terms of relative variation between the two calculation methods (Fig. 9b), confirm that the error induced by the equivalent-NaCl approximation has a relatively low impact on the dissolved CO 2 value (less than 5%) for the brines considered in our simulations. The maximum discrepancy on Figure 9b   Comparison between CO 2 solubility value (150°C in a 1 M NaCl solution) calculated with 1) PHREEQC, 2) PHREEQC with fugacity correction, 3) SCALE2000, and experimental data of Takenouchi and Kennedy (1964), (T and K).
of calcium (Kelly et al., 1986). This observation suggests that the small differences observed on Figure 9b are probably due to the fact that the majority of the brines were of NaCl type. However, a second comparison now performed with equivalent-NaCl solutions of the same ionic strength (thus including the ions charge inpact), shows that the approximated method systematically underestimates CO 2 solubility (Fig. 10a). The relative error induced by NaCl approximation (Fig. 10b) becomes significant for the majority of the brines (between 5% and 20%).
Consequently, these results show that the use of a simplified equivalent-NaCl approach may be acceptable only for NaCl-type waters of relatively low salinity. In the general case, the use of models accounting for the complete composition of the brine appears to be suitable. When not considering this correction, the error generated may reach significant values beyond the acceptable limit of 5%.
Our attempt of quantification of the relative weight of all the corrections described above is summarised in Figure 11. We arbitrary chose to work on a 237 g·l -1 brine at 200 bar    Synthesis diagram of the relative weights of the various corrections applied to the calculation of CO 2 solubility and mineral equilibrium for a 237 g·l -1 brine at 60°C and 200 bar. Indicated percentages are calculated with respect to the reference water composition "real brine, real gas". Quantity of dissolved CO 2 is expressed in mol·(kg H 2 O) -1 . and 60°C. Under these calculation conditions, with CO 2 treated as a real gas (fugacity coefficient of 0.45) assumed to occupy the whole gaseous phase, and taking into account the real brine composition (with all activity corrections), finally led to a CO 2 solubility of 0.49 mol·(kg H 2 O) -1 . Assimilating CO 2 to an ideal gas and neglecting salinity (thus considering pure water) induces an error on CO 2 solubility which is more than 400%! Figure 11 also shows that considering either salinity (NaCl equivalent) or non-ideality of CO 2 (applying fugacity correction) notably improves the prediction even though the calculated solubility remains largely overestimated (more than 100%). These results highlight the necessity to take into account both the fugacity correction and the salinity effect in order to get good estimates of CO 2 solubility. In comparison, the improvement provided by considering the real brine composition instead of its NaCl-equivalent remains relatively small in this case (less than 5%). With this dissolved CO 2 reference value of 49 mol·(kg H 2 O) -1 , the results on calcite solubility underline the impact of the aqueous activity model used: we calculated that 1 kg of dissolved CO 2 in the brine could lead to the dissolution of 71 cm 3 of calcite or 54 cm 3 of calcite using, respectively, PHREEQC or SCALE2000. The percentages indicated in Figure 11, though to be considered as exact only for the specific brine considered in the simulations, however represent good indicators of the relative impact of each of the corrections when considered separately. When all the corrections were simultaneously taken into account (calculations using SCALE2000), the results obtained could always be favourably compared to available experimental measurements, thus indicating that the improvement in the prediction was significant.

EXAMPLE OF APPLICATION TO CO SEQUESTRATION
In order to illustrate the improvements provided by the integration of fugacity and activity corrections in the simulation of gas-brine-minerals mass exchange processes, we studied the potential effects of CO 2 injection in a deep brine reservoir typical of some of those targeted for CO 2 geological sequestration. We considered here a hypothetical scenario of CO 2 injection in a Smackover-type formation (Arkansas, United State) for which detailed brine composition data were available in Moldovanyi and Walters (1992). The Smackover site is mainly characterised by a carbonateevaporite formation, fluid temperatures ranging from 51 to 105°C, in situ pressure of 150 to 400 bar, Na-Ca-Cl brines with TDS varying from 152 to 342 g·l -1 (Moldovanyi and Walters, 1992). For this study, we selected a 196 g·l -1 brine (sample no. 74 of Tables 1 to 3 in Moldovanyi and Walters, 1992) that was collected from a well located in the northern reservoir, at a depth of 1715 m and at a bottom hole temperature of 71°C. Its detailed chemical composition is given in Table 1 (column "A"). For the water-rock interaction simulations, we considered a hypothetical mineralogical composition (in weight: dolomite, 70%; calcite, 12%; anhydrite, 10%; quartz, 5%; siderite, 3%) consistent with carbonated reservoir conditions. For the following calculations, we assumed constant in situ pressure and temperature conditions of 180 bar and 71°C. Laboratory pressure and temperature conditions were supposed to be 1 bar and 25°C.

Initial Speciation Calculations
Simulating a massive CO 2 injection scenario first requires the definition of an initial reference state that we chose to be thermodynamic equilibrium between the brine and the mineralogical assemblage defined above. Prior to that simulation, it is then necessary to calculate the formation water speciation on the basis of the analysis data, in order to check the consistency between the analytical results (with inevitable measurement inaccuracies, electric charge imbalance, etc.) and the thermodynamic model results. These initial calculations were carried out with both SCALE2000 and PHREEQC to allow for comparison and the results are presented in details in Table 1. It has to be noticed that all the data presented in Table 1 (except temperature, pressure and most of the data in column "A" which are actually used as input in the models) are calculation results. Moldovanyi and Walters (1992) indicated however that their pH analyses must be viewed as approximations because of the potential for CO 2 degassing from brine samples. Consequently, we did not use the pH value as an input data in the input files (see Appendix 2 for the PHREEQC input file) and the carbonate system was constrained prescribing the total alkalinity (equivalent to At in Moldovanyi and Walters, 1992) and the total dissolved inorganic carbon (equivalent to TIC in Moldovanyi and Walters, 1992) that we considered to be more reliable input data (Moldovanyi and Walters describing a more careful analytical protocol for these specific measurements).
The first series of calculations aimed at verifying the global consistency of the input data set. For that purpose, we calculated with SCALE2000 and PHREEQC the brine characteristics at the laboratory conditions. In both codes, electro-neutrality was reached by automatically adjusting the total amount of chloride. The results detailed in Table 1 (columns "P1" and "S1") are both in very good agreement with the analytical data (column "A") in terms of total amounts of chemical elements except for chloride which was obviously slightly modified for electro-neutrality adjustment. This result is not of major interest in itself as the total amounts in column "A" (Table 1) were used as input data in the codes; it just confirms the correct behaviour of the models algorithms for converting initial input data in mg·l -1 into molalities (mmol·(kg H 2 O) -1 ). However, it has to be noticed that, contrary to PHREEQC, SCALE2000 did not use the density value (which is a key-parameter in these unit conversion calculations) as an input but recalculated it based on a formalism evaluating the volumetric properties of the aqueous phase using Pitzer's ion interaction approach (Monnin, 1994). It is then interesting to remark the very good agreement between the measured density at 25°C (1.132 kg·l -1 ) and the density calculated by SCALE2000 (1.134 kg·l -1 ) considering moreover that some of the trace elements present in the studied sample (B, F, I, Mn, Zn, Br, Li) are not accounted for in SCALE2000 speciation calculations (however, Li and Br initial amounts are taken into account in density calculation). Due to the fact that pH measurement was performed immediately upon collection of the fluid sample (thus at a temperature close to that of the reservoir, i.e. 71°C), the values of calculated and measured pH must not be compared directly even though they are very close between PHREEQC results (Table 1, column "P1") and analytical data. Lastly, the comparison between PHREEQC and SCALE2000 results ( Table 1, column "P1" and "S1") show some discrepancies for Saturation Indexes, pH, and CO 2 fugacity; provided that gaps between thermodynamic equilibrium constant used in the two codes are very weak, this can be easily explained by the activity models considered i.e. Davies in PHREEQC, and Pitzer in SCALE2000.
In a second step, we performed calculations under reservoir conditions (71°C, 180 bar) using PHREEQC ( Table 1, column "P2") and SCALE2000 (Table 1, column "S2"). It has to be noticed that the pressure effect on equilibrium constants (see Section 1.3) is not taken into account in PHREEQC. The pH value of 6.27 calculated with PHREEQC is close to the measured pH (6.30) whereas SCALE2000 predicts a higher pH value (6.53); however, analytical uncertainties on pH measurement previously mentioned do not permit to be conclusive on calculation accuracy, referring to this only pH criteria. However, on the basis of the results obtained with both SCALE2000 and PHREEQC and even though some slight discrepancies were observed with respect to the analytical data, it can be concluded that the calculated brine composition is globally consistent with the very detailed analytical data set provided by Moldovanyi and Walters (1992).
Once this initial water composition was characterised, we performed another speciation calculation with SCALE2000 assuming thermodynamic equilibrium with calcite, which is likely to occur in a carbonated reservoir. The results are presented in Table 1 (column "S3") and it is interesting to notice that, in this case, CO 2 fugacity is notably greater than in the case where calcite equilibrium was not supposed ("S2" water) while pH, alkalinity, and dissolved carbon are lower. These results are consistent with Moldovanyi and Walters' hypothesis of degassing during sampling. The loss in total carbon and calcium is explained by over-saturation of water "S2" with respect to calcite inducing calcite precipitation to reach thermodynamic equilibrium; the slight variations observed in Cl and Na concentrations (main components of the brine), expressed in molalities, are explained by the variation in the total amount of pure water (H 2 O) in the brine composition induced by the precipitation reaction.
The next initial calculation performed with SCALE2000 aimed at establishing a reference water composition to be used as an initial condition for the CO 2 injection simulations presented Section 2.2. For that purpose, we equilibrated the "S2" water with the selected mineral assemblage: calcite, dolomite, anhydrite, siderite, and quartz. To reach equilibrium, SCALE2000 predicted precipitation of dolomite and quartz (15.739 and 0.016 g·(kg H 2 O) -1 respectively) and dissolution of anhydrite, siderite, and calcite (0.367, 0.45, and 16.585 g·(kg H 2 O) -1 respectively) leading to significantly different concentrations for the major species involved in the precipitation/dissolution reactions (C, Na, Ca, Mg, Si, Fe in Table 1, column "S4"). Compared to the previous calculation (equilibrium with calcite alone, "S3" water), pH, alkalinity, and total carbon do not vary significantly.
In this particular example where salinity (TDS = 196 g·l -1 ) is more than five times higher than that of seawater, the importance of the activity model accuracy was underlined as calculation results obtained with PHREEQC (Davies activity model) and SCALE2000 (complete Pitzer's ion interaction model), even though comparable in terms of elementary composition, might lead to significant discrepancies when comparing pH, CO 2 fugacity, or saturation indexes. In order to corroborate the results presented in the first section of this article, we have performed a last series of simulations more particularly focusing on the fugacity correction effect. Starting from the initial water composition ("P2" and "S2", for PHREEQC and SCALE2000 calculations, respectively), we prescribed a high CO 2 pressure of 150 bar (this CO 2 pressure remaining constant for the whole simulation). The speciations obtained are detailed in Table 1 (column "P3" for PHREEQC and "S5" for SCALE2000) and it can be observed that neglecting the fugacity correction leads to an overestimation of the total dissolved carbon of more than 300%! Attempting to reduce this error, we prescribed in the PHREEQC input file (see Appendix 2) the CO 2 fugacity value calculated by SALE2000 (79.8 bar) instead of the partial pressure of CO 2 (150 bar); one must remark that, obviously, this fugacity value is generally not available as an input data for all conditions of potential interest. The results presented in Table 1 (column "P4") show that, despite this artificial correction, the error made on the total dissolved carbon still reaches 50%. This demonstrates the importance of the fugacity correction that was implemented in SCALE2000 when attempting to evaluate amounts of dissolved CO 2 ; moreover, for cases where the salinity of the fluid is above that of seawater (which are likely to be frequently encountered in many CO 2 injection studies), it is also of major importance to take into account in the models the combined effect of salinity, pressure, and temperature.

Effect of High pCO 2 on the Geochemical System
The CO 2 "injection" simulation first considered a closed reactor (batch) with a prescribed initial water-rock ratio of 1 to 4 in volume (namely 1 l of fluid for a total volume of 4 l of minerals) corresponding to the assumed initial porosity of 20%. In this reactor (total pressure: 180 bar; temperature: 71°C), we applied for the whole simulation duration a constant CO 2 partial pressure of 150 bar (leading to a fugacity value of 79.8 bar, according to SCALE2000 calculation). The idea was to estimate on a representative elementary volume of the geochemical system the effects of a high CO 2 pressure likely to be encountered in massive CO 2 sequestration scenarios. These calculations were carried out using the "kinetics" module of SCALE2000 where the incorporated kinetic law for mineral dissolution/precipitation is derived from the Transition State Theory (Lasaga, 1984) for which the effective reaction rate V net (mol·s -1 ) can be expressed according to the following general equation: where k is a kinetic constant expressed in mol·cm -2 ·s -1 , S is the total surface area (cm 2 ) assumed to be a constant term when precipitation occurs and proportional to the remaining mineral mass in the case of dissolution, fr is the fraction of the total surface area effectively reactive, a H+ is the H + activity term, Q is the ionic activity product, and K is the thermodynamic equilibrium constant of the considered precipitation/dissolution reaction. The reactive surface area, which is one of the most difficult parameters to characterise, was defined here on the basis of a classical spherical grain model; in our simulations, fr was set to 1 and no pH effect on kinetic was considered (n set to 0), whatever mineral species is considered. The kinetic constant k is calculated, for each mineral species, in function of the temperature according to the following law (Aagaard and Helgeson, 1982): In which E a is the activation energy and T r is a reference temperature. The main operational parameters for this simulation are presented in Table 2 (initial fluid composition being detailed in Table 1, column "S4").
The main results of this batch simulation are illustrated on Figure 12 where mineral concentration variations are expressed as relative mass variations calculated according to the following equation: In this equation, C i,min (t = 0) and C i,min (t) are, respectively, the initial concentration of mineral i and its concentration at time t, expressed in g·m -3 of saturated bulk rock. We preferred this unit of concentration over the "classical" g·(kg H 2 O) -1 as this latter unit may sometimes be a bit confusing for minerals. For example, considering a rigorously constant number of moles of a mineral (e.g.: quartz) in thermodynamic equilibrium with the fluid phase, variation with time of this mineral concentration may however occur if expressed in g·(kgH 2 O) -1 because of other progressing chemical reactions (e.g.: precipitation/dissolution of a hydrated mineral such as gypsum) inducing fluctuations of the H 2 O amount in solution. This concentration variation could then be erroneously interpreted as precipitation or dissolution of the non-reactive mineral (quartz in this example).
What can be seen first on Figure 12 is that the response time of the system to this strong initial input in CO 2 is relatively short as within a few tens of days, a new thermodynamic equilibrium state is reached. On the other hand, the mineralogical composition changes only slightly: the relative variations induced by the presence of 0.48 mol·(kg H 2 O) -1 of dissolved carbon (to be compared to 0.0026 mol·(kg H 2 O) -1 , which is the total dissolved carbon concentration of the brine initially equilibrated with the mineralogical assemblage, before CO 2 "injection") are of the order of 0.01% or less for precipitating minerals (quartz, dolomite, anhydrite) and less than 0.08% for dissolving ones (siderite, calcite). Another indicator of this low impact on the mineralogical assemblage is the negligible porosity variation observed on Figure 12f. Lastly, it can be observed on Figure 12g that pH remains constant at a value of 4.2 which is not surprising due to the Relative mass variation (%) =100 0 Modelling results of the batch system under a prescribed pCO 2 of 150 bar. Evolution with time of mineral phases relative mass variation (a-e), porosity (f), and pH (g). buffering effect induced by the presence of carbonates; however, looking more in details at the result file show a rapid initial increase in pH from 3.1 to 4.2 in less than 0.5 day (not visible on Figure 12g). This is a response to the dramatic initial change in CO 2 partial pressure varying instantaneously from 0.382 bar (initial "S4" water in Table 1, corresponding to a fugacity of 0.203 bar) to 150 bar (fugacity of 79.8 bar) and to the simultaneous dissolution/precipitation of carbonate minerals.
However, this batch simulation has to be considered as a first convenient way of evaluating qualitatively the potential impact of a high dissolved CO 2 concentration on a geochemical system. By definition, such a closed reactor at fixed pressure and temperature conditions, has a limited capacity to dissolve CO 2 and once gas-water-mineral equilibrium is reached, nothing else is going to happen. Moreover, CO 2 injection in an aquifer is by nature a dynamical process. It is then essential to be able to consider the combined effects of both fluid dynamics and chemical reactions, which is much more intricate in terms of physical processes and numerical aspects. SCALE2000 enables the simulation of coupled chemistry and mass transfer using a specific approach based on a set of perfectly mixed reactors connected in series. Originally, this option was implemented for simulating scale deposition along production or injection wells in oil or geothermal fields. Further details on this concept which enables 1D chemistry-transport coupled modelling, and on the basic equations, can be found in Azaroual et al. (2004aAzaroual et al. ( , 2004c. A diagram of the set of reactors and the initial and boundary conditions are presented on Figure 13. For comparison with batch calculation, we first focused on a single reactor (reactor no. 1) identical to the batch reactor presented above in terms of total volume (5 l), initial water-rock ratio (1:4 in volume), and initial chemical composition (fluid + minerals). Similarly to what was done for the batch reactor, we prescribed a constant pCO 2 of 150 bar in reactor no. 1 for the whole simulation duration (50 y). However, we considered here that the system was continuously fed by an injection fluid at a constant flow rate (corresponding to an initial pore water velocity of 1 m·day -1 ), the chemical composition of this injection fluid being identical to that of the brine initially present in the reactors (i.e. fluid "S4" in Table 1).
All things being equal initially, the results presented on Figure 14 demonstrate the impact of mass transport on the geochemical evolution of the system: the continuous fluid renewal has a dramatic effect on the geochemistry of the system and no steady state can be reached after a simulated Injection of brine S4 at constant flow rate Temperature = 71°C ; total pressure = 180 bar -Temperature = 71°C -Total pressure = 180 bar -Porosity = 20% -Temperature = 71°C -4 l of {dolomite (70.04%) + calcite (12.69%) + anhydrite (9.67%) + quartz (5.41%) + siderite (2.18%)} -1 l of brine S4 (see Table 1  Diagram of the open-reactor system constituted of four identical homogeneous reactors connected in series. This configuration enables to simulate 1D chemistry-transport coupled processes. Initial and boundary conditions are detailed. Modelling results of the open-reactor system for reactor no. 1 in which pCO 2 is prescribed at a constant value of 150 bar. Evolution with time of mineral phases relative mass variation (a-e), porosity (f), and pH (g). Modelling results of the open-reactor system for reactor no. 4 (most downstream) in which pCO 2 is not prescribed. Evolution with time of mineral phases relative mass variation (a-e), porosity (f), and pH (g).
period of 50 y, whereas thermodynamic equilibrium was reached in a few tens of days only for the batch reactor case. The mineral composition now appears to be strongly modified with respect to the initial state (detailed in Table 2): all carbonated minerals (calcite, siderite, and dolomite) are exhausted from the system in less than 35 y (Figs 14b, 14c,  14e). The disappearing of all carbonates results in a loss of the buffering capacity of the system illustrated by the sudden decrease in pH from 4.2 to 3.14 ( Fig. 14g) observed as soon as the last carbonate species (dolomite) has totally dissolved; pH is then directly controlled by the high partial pressure of CO 2 (150 bar). In the meanwhile, quartz can be considered as inert even though some negligible precipitation can be observed (Fig. 14a). Anhydrite precipitates continuously at a constant rate (Fig. 14d) independently of carbonates dissolution and pH variation. This indicates that the common ion effect (Ca provided by calcite and dolomite dissolution) does not influence the precipitation rate of anhydrite. As a consequence of these dissolutions and because amounts of precipitated anhydrite remains relatively weak, porosity strongly varies from an initial value of 20% to above 85% within the first 30 y of simulation and then remains approximately constant. This very high increase in porosity is in contrast to the mineral trapping effect that one might expect to occur when injecting CO 2 in a reservoir. However, being in a carbonated reservoir under constant high CO 2 pressure conditions, this behaviour could have been predicted beforehand. Nevertheless, it is interesting to compare the behaviour of the system of reactors studied here to the results obtained in Bonijoly et al. (2003). In this work, 1D finite difference coupled simulations of CO 2 injection in the carbonated Dogger aquifer (Paris Basin, France) were performed and high rates of calcite dissolution, leading to a huge increase in porosity, were also observed in the immediate vicinity of the injection well. At the same time, calcite reprecipitation (thus inducing some CO 2 trapping) was observed a few tens of metres upstream the injection well (where it has to be noticed Modelling results of the 1D system constituted of a set of four identical homogeneous reactors connected in series. Relative variation of mineral carbon mass (mole of representative elementary volume of saturated bulk rock) in the four reactors of the system (from upstream to downstream). In reactor no. 1 (most upstream), pCO 2 is prescribed (150 bar); in the three other reactors, pCO 2 is not prescribed. that the temperature was higher). In the present study with SCALE2000, where calculations were carried out on a set of four identical reactors connected in series (pCO 2 being prescribed in reactor no. 1 only) at constant temperature and pressure, it is interesting to remark that in the mostdownstream reactor (i.e. reactor no. 4, see Fig. 13), while siderite and calcite also disappeared (with noticeable retardation with respect to reactor no. 1), significant dolomite precipitation can be observed as illustrated on Figure 15, suggesting some possible carbon mineral trapping. In reactor no. 4, porosity variations are almost negligible in comparison to what was observed in reactor no. 1. In order to investigate the possible trapping capacity of the studied system, we performed carbon mass balance calculations for the whole simulated period (50 y). For that purpose, we calculated in the four reactors the total amounts of carbon trapped in the mineral phases. Figure 16 show the corresponding results expressed in terms of relative mass variation of carbon stored in carbonates (calcite, dolomite, and siderite) with respect to the initial amounts (as detailed in Table 2). Confirming what was observed in Figure 14, all mineral carbon was removed from reactor no. 1 (variation of -100%) after 35 y. In reactor no. 2, only 50% of the initial carbon in solid was dissolved; however, the general trend suggest a similar tendency than in reactor no. 1, i.e. total dissolution of mineral carbon, with a retardation factor of about 35 y. More surprising are the results obtained for reactors no. 3 and no. 4 where mineral carbon dissolution is here considerably less important (less than 1% of the initial mass), due to significant dolomite precipitation (see Fig. 15). The net carbon mass balance over the four reactors is thus largely negative, however, the trend shown in the last two reactors may let us expect a positive mass balance (i.e. actual carbon trapping in carbonates) beyond the fourth reactor considered here. As mentioned earlier, this behaviour can be favourably compared to modelling results obtained by Bonijoly et al. (2003) where carbonates precipitation was predicted a few tens of metres downstream the injection well. The space scale considered here (a set of four 5 l total capacity reactors) is obviously incomparable to the site scale. Consequently, further modelling work should be carried out at a scale of at least a few hundred metres around the "injection point", thus requiring the use of more sophisticated numerical tools enabling 1D, 2D, or even 3D hydrodynamicstransport-chemistry coupled modelling. Relying on the coupling approach presented in Kervévan et al. (1998), this could possibly be achieved through interfacing the SCALE2000 geochemical calculation modules with a 1D, 2D, or 3D hydrodynamics and mass transport code. Provided that a detailed mineralogical description of the site would also be available, it might then be possible to draw more definitive conclusions on the actual mineral trapping capacity of a carbonated formation such as Smackover.

CONCLUSION
This work presents some improvements for brine-rock-gas interaction modelling through the inclusion of thermodynamic approaches specifically dedicated to high salinity, high pressure, and high temperature conditions. These formalisms were all implemented in the thermo-kinetic modelling software SCALE2000 (Azaroual et al., 2004a). They mainly deal with activity correction for dissolved species (25-250°C, 1-1000 bar, ionic strength up to 6 M) and density calculation involving Pitzer's formalism (Monnin, 1994), the taking into account of the influence of both temperature and pressure (particularly at high pressure conditions) on thermodynamic equilibrium constants, and fugacity coefficient correction based on Duan et al's. (1992aDuan et al's. ( , 1992bDuan et al's. ( , 2003 equations of state. They enable to calculate more accurately CO 2 solubility in brines, which is obviously a key parameter in the context of CO 2 injection studies.
A quantification of the relative weight of each of these corrections was proposed and it appears that fugacity and activity corrections have the greatest impact on the precision of the geochemical calculations applied to a CO 2 -brineminerals system under conditions of relatively high temperature, pressure, and salinity, which are likely to be encountered at potential CO 2 sequestration sites.
On the basis of these developments recently integrated in SCALE2000, we presented in the second part of this paper an example of application to a hypothetical scenario of massive CO 2 injection in a carbonated formation, considering actual geochemical data from the Smackover site (Arkansas, United State) for the brine composition (Moldovanyi and Walters, 1992). The purpose of these calculations was to evaluate the geochemical impact on both mineral and fluid compositions of a high prescribed CO 2 pressure (150 bar). The first series of simulations for a representative elementary volume of saturated bulk rock considered as a closed system (batch reactor), showed a relatively short response time (a few tens of years) but finally demonstrated a relatively low impact of the presence of CO 2 on the mineral composition. A second series of simulations performed on the same initial representative elementary volume of bulk rock, now considered as a single open reactor crossed by a constant water flow rate, led to dramatic changes in the mineral composition at the end of our 50 y simulation: carbonated minerals (calcite, dolomite, siderite), initially present in large amounts, totally disappeared while anhydrite precipitated continuously. As a consequence, a strong increase in porosity was observed (variation from 20 to 85%) thus suggesting a poor mineral trapping capacity of the system. However, further calculations carried out with SCALE2000, now considering an approximated 1D system constituted of a set of four homogeneous identical reactors, showed significant dolomite precipitation in the most-downstream reactor and hence, some CO 2 precipitation. First mass balance calculations performed on the four reactors system finally showed a global loss in total mineral carbon with respect to initial conditions. However, the trend in the most-downstream two reactors indicated that possible trapping might be expected beyond the considered geometrical limit of our modelled system. This behaviour can be favourably compared to observations made in the study by Bonijoly et al. (2003) on CO 2 injection in the carbonated Dogger aquifer (Paris Basin) where high porosity increase was also predicted near the injection well while carbonates precipitation was calculated by the 1D coupled model a few tens of metres downstream. However, drawing more definitive conclusions on the mineral trapping capacity of the Smackover formation would necessitate further simulations at the site scale considering, in particular, a more detailed mineralogical description and a sophisticated hydrodynamic and mass transfer model (1D, 2D, or 2D) more adapted to deal with complexity of the physics of gas injection and flow in aquifers.

APPENDIX 1
PHREEQC input file used for the series of calculations performed considering a 1 M NaCl solution at 150°C and pressures varying from 50 to 400 bar (see Section 1.5 for details). Corresponding results are represented on Figure 8.
Database used is llnl.dat (revision 1.11) based on the Davies activity model which is enabled adding "#" at the beginning of every "-llnl gamma" command line.