Fabrizio Gherardi

Pascal Audigane

Eric C Gaucher

Predicting long-term geochemical alteration of wellbore cement in a generic geological CO 2 confinement site: tackling a difficult reactive transport modeling challenge

Keywords: reactive transport modeling, CO 2 geological storage, cement, caprock, Paris Basin

come

Predicting long-term geochemical alteration of wellbore cement in a generic geological CO2 confinement site: tackling a difficult reactive transport modeling challenge

Introduction

In recent years, many studies in subsurface carbon dioxide (CO 2 ) storage in geological formations have been performed to evaluate the leakage risks connected to the different storage options under consideration (saline aquifers, depleted gas and oil reservoirs, unminable coal seems; e.g., among many others, Bachu et al., 2007Bachu et al., , 1994;;Gaus et al., 2005;Bachu, 2000;Hitchon et al., 1999;Freund and Omerod, 1997;Holloway, 1997). Laboratory experiments and studies on artificial and natural analogues of cement/rocks interactions indicate that long-term wellbore integrity is a crucial aspect to be investigated when assessing the impact of CO 2 geological storage (IPCC 2005(IPCC , 2008)). One of the most relevant issues is whether the petrophysical and geochemical properties of the host rocks and sealing materials (i.e. caprock and cement of wellbores) will remain constant over long periods of time (hundreds to thousands of years). Abandoned wells are a major concern because of the inherent high reactivity of cements. When CO 2 encounters a well without proper plugging, it will tend to migrate upwards due to buoyancy force. This rising plume chemically interacts with aquifer, caprock and cement minerals predominantly through the wetting aqueous phase, which has low pH and high mineral dissolution potential, due to the solubilization of high amounts of CO 2 . The interaction between these fluids and the rocks through which they migrate is expected to alter the initial minerals and lead to the formation of secondary minerals, resulting in changes in physical and chemical properties of the system. Possible leakage may occur through preferential flow pathways along rock-cement and casing-cement interfaces, or through degraded materials because of mechanical failure and/or chemical degradation of rocks, wellbore casings and/or cementitious materials. Recent laboratory experiments and field observations show evidences about the type of geochemical modifications induced in the cement after exposure to CO 2 (see among many others Duguid and Scherer, 2010;Rimmelé et al., 2008;Jacquemet et al., 2008;Carey et al., 2007;Kutchko et al., 2007). A large part of these results are specific to T,P and geochemical conditions of the experiments and/or of the field observations, and they are not directly usable outside these conditions to infer possible evolutionary scenarios in other sites under consideration for CO 2 geological confinement. In addition, laboratory experiments are unable to produce data for long periods of interaction. Numerical modeling is then the only effective means to predict the long term evolution of the cement plug and to evaluate the integrity of wellbore systems designed to prevent fluids from escaping the sequestration reservoir.

Simulating reactive fluid flow and transport under these conditions is a challenging task because of the complexity of multiphase fluid flow, water-gas-rock and fracture matrix interaction mechanisms, and difficulties dealing with the heterogeneous physical and chemical properties of the medium. Additional difficulties stem from the inherent chemical complexity of the investigated system, related to the concomitant migration of two compositionally different plumes: an acidic, Crich plume propagates from the reservoir upwards into the seal, whereas an alkaline pulse, caused by the equilibration of local groundwater with cement minerals, diffuse into the adjacent caprock and the underlying reservoir.

The objective of this paper is to investigate, by numerical modeling, the long-term durability of the cement of the pre-existing wells of a generic site for CO 2 geological storage in the Paris Basin (France). The target reservoir is hosted in the Middle Jurassic (Dogger) regional aquifer, in the centre of the Paris Basin, at a depth of 1500-2000 m. The Dogger aquifer, already exploited for geothermal and petroleum resources, is mainly composed of carbonates. Sedimentary studies (Rojas et al., 1989) have shown that the most permeable horizons of the Dogger aquifer correspond to high-porosity layers of oolitic limestones which still retain, or even increase, their initial porosity by mechanical fracturing and/or chemical dissolution. Recently, several French research projects have been conducted to investigate the potential of this system for carbon dioxide geological storage purpose. The "PICOREF" project (Pilote pour l'Injection de CO 2 dans les Réservoirs géologiques, En France, 2006France, -2008) ) aimed, in particular, at developing a methodology of site selection for CO 2 geological storage in deep saline aquifers and studying the injection of CO 2 in specific contexts of the Paris Basin (Brosse et al., 2007). The "GeoCarbone-Injectivity" and "GeoCarbone Integrity" projects have been funded by the French National Agency of Research (ANR, see http://www-anr-co2.cea.fr). The main goal of the "GeoCarbone-Injectivity" project was to develop a methodology to understand and predict the injectivity evolution of a CO 2 well during the injection operations in saline aquifers (Lombard et al., 2007), while the project "GeoCarbone Integrity" provided an integrated approach to evaluate sealing efficiency (Fleury, 2007). One work package of this project is dedicated to numerical modeling assessment of cap rock integrity (Bildstein et al., 2009).

In a very similar fashion, our study aims at providing preliminary estimates of the geochemical processes occurring at the interface between (i) reservoir, (ii) cement and (iii) caprock porous formations under isothermal, water-saturated and diffusively-controlled mass transfer conditions with a specific attention on long term well cement integrity. The aim of the simulations is to provide a sort of reference scenario for the processes likely occurring at the top of the injection target reservoir over a time span of 1000 years, which lies reasonably in the middle of the range (hundreds to thousands of years) required for long-term integrity modeling of geological CO 2 sequestration (IPCC 2005, 2008, Gaus et al., 2008). The present simulations are prototypical in nature because no specific information is available on cement mineralogy and pore water chemistry of the abandoned wells in the Paris Basin. In contrast, field and laboratory data are available to constrain the initial hydraulic and geochemical conditions for reservoir and caprock domains. The effect of highly uncertain parameters, such as tortuosity, or intrinsically variable initial conditions, such as the initial chemical composition of porewaters, has been assessed by performing a set of different simulations with selected parameters calibrated on the basis of experimental data or numerical results of previous simulations.

Modeling approach

All calculations were performed with the TOUGHREACT code, a comprehensive, publicly available reactive transport simulator for multiphase systems (Xu et al., 2004). The modeling of reactive transport is based on space discretization by means of integral finite differences (Narasimhan and Witherspoon, 1976), and is carried out by using a sequential non iterative approach. The chemical transport is solved on a component basis, and the resulting concentrations obtained from the transport are substituted into the chemical reaction model. The system of chemical reaction equations is solved on a grid-block basis by Newton-Raphson iteration. Full details on numerical methods are given in Xu and Pruess (2001) and Xu et al. (2006).

Simplified geometries with hydraulically and geochemically homogeneous media have been taken into consideration, and fully saturated (i.e. aqueous phase saturation, S l =1.0) and isothermal conditions at 75°C assumed throughout our calculations. These highly idealized conditions provide a reasonable approximation to average field conditions far from CO 2 injection areas in the Paris Basin, at the top of the Dogger reservoir (Figure 1). Previous simulations have indicated that a supercritical CO 2 bubble is likely expected to occur in the first 5 m around the injection well after 10 years injection, and that this zone is surrounded by a two-phase zone of about 650 m (André et al., 2007).

Model implementation has been carried out in stages of increasing complexity. We started from geometrically simplified systems, i.e. 0-and 1-dimensional models, to explore the sensitivity with respect to different space discretizations, initial mineralogy and chemical composition of pore waters. Nearly equilibrated, steady-state geochemical systems were created by running batch simulations, and the resulting mineralogical and pore water chemical compositions used as initial conditions for the 1-D and 2-D simulations. Simplified 1-D models have been preliminary used to investigate the processes occurring at reservoir-cement, reservoir-caprock and caprock-cement interfaces. Space and time discretization of the 2-D model results from a compromise between computation time, simulation time, grid spacing and numerical stability. Numerical stability has been achieved through von Neumann analysis (Hindmarsh et al., 1984) by enforcing small maximum time-step increments, variable between 5×10 4 to 1×10 3 s depending on initial conditions and geometry. Maximum simulation times were typically in the order of few to several days (up to a maximum of about 15 days) on a 3.6 MHz CPU. The results of a 2-D reference case will be presented in the next sections.

Geometry

In an attempt to keep the implementation of the model as simple as possible, a two-dimensional, idealized Cartesian grid has been arranged to reproduce a vertical cross-section of a heterogeneous porous medium sampled in proximity of the interfaces between cement, caprock and reservoir domains (Figure 2). This section is drawn through the axis of a hypothetical abandoned well, and because of the symmetry of the model, only half of the domain has been modeled. The grid extends horizontally 5 m from the axis of the abandoned well and vertically 10 m from the top of the reservoir. A 2 m-thick reservoir with large-volume blocks at the bottom is placed below cement and caprock domains so as to guarantee an almost infinite supply of chemicals from the carbon-rich, low-pH waters likely permeating the reservoir after massive injection of CO 2 at depth. The reservoir is further discretized into 8 layers between the interface with cement and caprock domains (at the top) and this large volume (at the bottom of the model), in order to account for possible geochemical variations induced in the reservoir by the diffusion of alkaline and near-neutral fluids coming from the cement filling and the cap rocks, respectively. A refined spatial discretization, with internodal distances of 0.05 m, is adopted near the interfaces to better resolve concentration gradients between the different rock domains. A coarser spacing is applied in the rightmost part (along the horizontal direction, far from cement-caprock interface) and in both the uppermost and lowermost parts of the grid (along the vertical direction, far from reservoir-caprock interface). In these zones, the maximum internodal distance slowly increases up to 0.75 and 1.76 m along the vertical and horizontal directions, respectively. This non-uniformly spaced cartesian grid responds to the double need to keep the highest resolution localized where the highest concentration gradients occur, i.e. near the interfaces, and to reduce the computation cost of the simulations (by reducing the total number of cells). Exploratory, simplified simulations (only unreactive aqueous tracers 1 have been allowed to propagate all over the modeled domain to reduce simulation times) with uniform and non-uniform grids confirmed that the choice of increasingly larger blocks with increasing distances from the interfaces does not cause any appreciable dispersion effect. Selected 1 Aqueous tracers are unreactive chemical components characterized by negligible molecular weight and molar volume. In our model, a fictitious chemical component, with initial concentration of 10 -6 mol kg H2O -1 , has been used as aqueous tracer for each of the three rock domains considered in the simulations. relevant parameters on grid geometry are shown in Table 1. Eleven observation points have been used to analyze the simulation results (Figure 2B). The mesh element CEM-01 is part of the cement domain and is located at the reservoir-cement interface. It is directly connected to the RES-01 element which is part of the reservoir domain. The CEM-02 element is connected to the CEM-01 element in the Z direction. In a similar fashion, elements CEM-03, CEM-04 and RES-02 illustrated the reservoir-cement interface where cement, reservoir and caprock domains are in contact.

Elements CAP-01, part of the caprock domain, and element RES-03 part of the reservoir domain are used for the understanding of the interactions at the caprock-reservoir interface. Element CAP-02 is connected to the CAP-01 element in the Z direction. Finally, elements CEM-05 and CAP-03 are used for the analysis of the geochemical reactivity occurring at the caprock-cement interface in the upper part of the simulated area. A no-flux boundary condition applies to the borders of the model. Our simulations are focused on a hypothetical abandoned well plugged with cement, and then any possible chemical effect induced by the corrosion or the degradation of the steel of the casing (e.g. Hope et al., 1985;Alonso et al., 2000;Glass and Buenfeld, 1997;Hausmann, 1967) is not taken into account.

Initial geochemical system

Geology, petrophysical and geochemical characteristics of the Dogger aquifer, in the Paris Basin, have been extensively investigated in the past (Rojas et al., 1989). The Paris Basin is a Late Tertiary intracratonic basin with a diameter of 600 km and a maximum thickness of 3000 m, bordered by Palaezoic massifs (Vosges to the east, Armorican Massif to the west, Central Massifs to the south, Ardennes in the northeast). The reservoir is hosted in limestones and dolomite-bearing formations.

The most permeable part of the aquifer is represented by oolitic and bioclastic limestones predominantly made of carbonates (up to 85% by vol.) and minor alumino-silicates (albite, Kfeldspar) and clays (Rojas et al., 1989). Several large aquifers are present at depths that are being exploited for potable supply and geothermal purposes. Maximum temperatures range up to 85 °C at a maximum depth of -1900 m below sea level in the central part of the Basin, some 20 km east of Paris.

The geochemistry of this aquifer has been extensively investigated by means of isotope, noble gas and hydrochemical studies (Michard and Bastide, 1988;Criaud et al., 1989;Coudrain-Ribstein and Gouze, 1993;Matray et al., 1994;Azaroual et al., 1997;Marty et al., 2003). Pore waters have different salinities ranging between 5 to 35 g kg H2O -1 . An aqueous solution with a chemical composition similar to that of Fontainebleau waters (Michard and Bastide, 1988), but higher total salinity (ionic strength, IS, of about 0.3 mol kg H2O -1 instead of about 0.1 mol kg H2O -1 ), has been preliminary considered as representative of the average initial conditions in the reservoir, before CO 2 injection at depth. Following André et al. (2007), zero-and mono-dimensional models have been successively run using this water and a simplified primary mineralogy to evaluate the geochemical effects induced by the injection of CO 2 at depth. The modified mineralogy and a new chemical composition of the pore water, characterized by the lowest pH values calculated during these simulations, have been then used to initialize the 2-D models (Tables 2 and3). et al., 2007;Landais, 2006) are available to constrain the initial conditions of our model.

Mineralogy and pore water chemistry of Callovian-Oxfordian argillites has been thoroughly investigated by experimental methods and numerical simulations (Sammartino et al., 2003;Vieillard et al., 2004;Gaucher et al., 2004a;2006;2009;Leroy et al., 2007;Vinsot et al., 2008).

The Callovian-Oxfordian argillite is a complex mineral assemblage predominantly composed of phyllosilicates, silicates, carbonates and S-bearing minerals (Claret et al., 2002;Claret et al., 2004;Gaucher et al., 2004a 2 and3, respectively.

An idealized composition for a high alkali (Na TOT and K TOT aqueous concentrations around 3 × 10 -2 m), ordinary Portland cement (OPC) has been used in the calculations. From a chemical standpoint, OPC can be considered as a very heterogeneous material with discrete particles typically in the size range of about a few hundred nanometers up to a few hundred micrometers mainly consisting of variable amounts of calcium silicate hydrates, portlandite, calcium aluminates, and minor reactive phases (e.g., hydrotalcite, ettringite, hydrogarnet, ferrite etc.), depending on ambient temperature and initial chemical composition of the cement paste (Paul and Glasser, 2000). A thermodynamically stable mineralogical assemblage that consists of portlandite, tobermorite, katoite-Si, hydrotalcite, ettringite and C 3 FH 6(s) (Table 3), and the associated pore water chemistry (Table 2) have been reconstructed at 75°C by numerical calculations. This mineralogical assemblage reflects the relative stability of the minerals included in the thermodynamic database and is coherent with some of the few available information on prolonged cement hydration at high Kusachi et al. (1989) indicated that afwillite is expected to be thermodynamically stable at about 100°C, rather than at 75°C. For all these reasons, the association tobermorite + portlandite, whose thermodynamic properties are well defined at 75°C, has been retained in our calculations.

Redox conditions are poorly constrained in our system, and for this reason, some redox-sensitive minerals of uncertain occurrence, such as hematite and magnetite, have been excluded from calculations. In contrast, ferrihydrite and goethite have been included because widely reported in the literature as possible secondary phases related to circulation of hyperalkaline fluids in clay-rich formations (e.g. Gaucher and Blanc, 2006). We initially set an uniform Eh value of about -210 mV for all pore waters considered in our model, calculated by batch modeling from data available for COx formation (around -175 mV at 25°C; Gaucher et al., 2006). Slightly different values have been then used for cement and reservoir redox conditions in the 2D model (Table 2), after preliminary 0D and 1D simulations carried out to achieve nearly equilibrated, steady state geochemical conditions.

Methods for geochemical calculations

Aqueous species are assumed in chemical equilibrium throughout the calculations. The detail of primary and secondary species considered in the simulations is given in Table 2. Aqueous activity coefficients are computed using an extended Debye-Hückel equation, according to Helgeson et al.

( 1981) and Tanger and Helgeson (1988), with minor adaptations (André et al., 2007) to account for the high pressure conditions of the system (around 200 bar). All the dissolution and precipitation reactions involving cement minerals are assumed to proceed at local equilibrium. This approach is justified by the observation made by previous investigators that chemical transformations occurring at the alteration front in cementitious materials are fast if compared with the time frame of hundreds to thousands of years required in the modeling of the deterioration of containment materials of radioactive or non-radioactive waste repositories (e.g. Gaucher et al., 2004b;Burnol et al., 2006;Trotignon et al., 2007), and with the slow rate mass transfer conditions assumed in the model (see sections below). The lack of reliable information on kinetic parameters for these minerals further supports this assumption. The local thermodynamic equilibrium approach has been applied also to simulate the behavior of other fast reacting minerals, such as calcite and anhydrite (Palandri and Kharaka, 2004).

Dissolution/precipitation reactions involving the remaining minerals from reservoir and caprock domains are modeled under kinetic constraints. A kinetic formulation derived from transition state theory (Lasaga, 1981(Lasaga, , 1984;;Aagaard and Helgeson, 1982;Steefel and Van Cappellen, 1990) is applied. Our model takes into account the feedback between transport and chemistry (see section 3.5 below) by allowing variations in porosity to be computed from the algebraic sum of mineral volume fractions variations.

Thermodynamic data

The primary mineralogy of reservoir, caprock and cement domains is described in terms of 7, 8 and 6 minerals (Table 3), respectively, and the complete mineralogical dataset is composed of 56 minerals. The choice of this large number of minerals responds to two major needs: accurately describe the geochemical complexity of the three media considered in the simulations and avoid the risk of converging towards an unique, biased numerical solution (e.g. Gaucher and Blanc, 2006;Marty et al., 2009). A complete list of minerals and geochemical reactions included in the simulations is given in (1992) and Damidot and Glasser (1995).

Amorphous CSH solid solution is discretized into a set of three pure phases with Ca/Si ratios varying between 0.8, 1.2 and 1.6.

Several assumptions and simplifications in the primary and secondary mineralogical assemblage have been made to account for the inherent complexity of the investigated system, and the problematic thermodynamic representation of all the phases possibly occurring during the calculations. Major approximations are related to the characterization of highly undifferentiated minerals like clays, zeolites, micas, and amorphous phases. The existence of solid solutions has not been taken into account. The following simplifications hold for these calculations.

Kinetic data

The kinetic approach is still inherently problematic in reactive transport modeling because of the partial knowledge of the large complexity of factors on which experimental rates have been proven to depend on. Major flaws of kinetic approach rely on the lack of reliable information on the values of mineral reactive surface areas and their variation with time, of mineral surface roughness, on the variation of rates over large span of Gibbs free energy of reaction, r G , speciation conditions of the associated aqueous phase and reaction temperature, and the limited use of laboratory data for predicting the kinetic reactions in natural systems (see, among many others, Yang and Steefel, 2008;Hellmann and Tisserand, 2006;Köhler et al., 2005;Zhu, 2005;Lüttge et al, 2003;White and Brantley, 2003;Yokoyama and Banfield, 2002;Gautier et al., 2001;Cama et al., 2000;Cama et al., 1999;White and Brantley, 1995;Ganor et al., 1995;Oelkers et al, 1994;Lasaga et al., 1994;Burch et al., 1993;Casey and Sposito, 1992;Anbeek, 1992;White and Peterson, 1990;Velbel, 1989Velbel, , 1993; and references cited therein). In addition, mineral reactive surface areas cannot be directly measured in nature, but only roughly estimated through laboratory experiments or field evidences.

This uncertainty also affects the rates estimated in the laboratory because of the normalization procedures required to allow their extrapolation to natural systems (e.g. Hodson, 2006).

Total surface areas, usually estimated in the laboratory by gas adsorption techniques (Brunauer et al., 1938), are often used instead of reactive surface areas in numerical calculations. However, as mineral dissolution commonly occurs at selected sites, they likely overestimate the effective reactive area (Aagaard and Helgeson, 1982). Many authors approximate total surface areas of geologic media using simple geometrical models involving the ordered arrangement of spheres (Bear, 1972;Canals and Meunier, 1995;Kieffer et al., 1999), but this approach require the availability of information, not available in the present study, on the statistical distribution of the equivalent diameters of mineral particles. In addition, the applicability of this approach is questionable in the case of shales or other rocks containing significant clay amounts, where minerals occur in many different habits (platy crystals, interstratified lamellae, etc.) other than spheres.

Summing up, despite the large number of investigations already performed, a coherent and complete kinetic database for all pertinent minerals is then still far from being available to date. In addition, a comprehensive sensitivity analysis of all kinetic parameters is not feasible when simulating very complex geochemical systems such as that investigated in the present paper.

To circumvent all these limitations, and in consideration of the prototypical nature of our simulations, we used selected kinetic parameters from the literature (Table 4), and performed a sensitivity analysis for the most relevant minerals only (see following sections). The precipitation of secondary minerals is represented using the same kinetic expression and rates as that for dissolution, except for kaolinite (Yang and Steefel, 2008), SiO 2 -polymorphs (Rimstidt and Barnes, 1980) and dolomite (Arvidson and Mackenzie, 1999). Catalysis by H + or OH -is considered only for mineral dissolution and the related kinetic parameters have been taken from selected literature (see footnote of Table 4).

Initial values for reactive surface areas have been set to two single average reference values (RSA 1 and RSA 2 ) taken from the literature (e.g. Köhler et al., 2005;Bickmore et al., 2002;Brantley and Mellott, 2000;Metz et al., 2005;White and Brantley, 1995), in order to account for the larger intrinsic surface areas of phyllosilicates and zeolites (RSA 1 =5 m 2 g -1 ) with respect to other minerals (RSA 2 =10 cm 2 g -1 ). The two average values have been selected at the lower end of the range of values reported in the literature to account for the concentrated reactivity along the edges of the phyllosilicates (e.g. Metz et al., 2005;Tournassat et al., 2003;Bickmore et al., 2002Bickmore et al., , 2001)), the ageing effect which allows in situ phases to be progressively less reactive in time due to adsorption of ions, the smoothing of surface defects by recrystallization, and the secondary minerals coating on primary minerals surface (e.g. Lasaga, 1998).

The dissolution and precipitation reactions involving fast reacting minerals, such as calcite, anhydrite and cement minerals are assumed to proceed at local equilibrium (see chapter 3.3). For these minerals, no kinetic parameters were provided in Table 4.

Mass transport mechanisms

Molecular diffusion in a single aqueous phase has been considered the dominant mass transport process in our simulations. This corresponds to assuming almost stagnant conditions all over the modeled system. This assumption could be not strictly valid for the reservoir domain, where, under field conditions, relatively high porosity and permeability conditions would allow mass transfer to occur also by advection in response to local hydraulic gradients (e.g. Bear, 1972). The treatment of the mass transfer process is based on Fick's First Law. The diffusive flux is assumed to be proportional to the concentration gradient:

x C d J i i eff i [eqn.1]
where i C is the concentration of the diffusing species in solution and i eff d is its effective diffusion coefficient. The latter coefficient can be thought as the product of the diffusion coefficient of the considered species in an aqueous solution at infinite dilution and a formation factor which accounts for the spatial variations in porosity and the tortuosity of the paths through the considered media.

Definitions of tortuosity vary in the literature (e.g. Lichtner et al., 1996). In TOUGHREACT tortuosity is a parameter which allows an effective diffusion coefficient in the porous medium ( eff d ) to be calculated from the diffusion coefficient in pure water ( w d ) according to the following equation:

w eff d d [eqn.2]
where is the porosity and the tortuosity coefficient. The tortuosity coefficient is always less than or equal to 1, and may vary in the different rock domains.

An alternative formulation for the coefficient for molecular diffusion in porous media is based on Archie's Law (Bear, 1972). Formation factors generated from Archie's Law have associated large uncertainties, of at least one order of magnitude, that largely depends on the anisotropic distribution of the formation factor in geological media and on its sensitivity with respect to different rock fabrics and grain alignments (Lichtner et al., 1996, and references therein).

Due to the lack of specific reliable experimental data and the intrinsic complexity of accurately describing the influence of the pore structure on macroscopic transport properties of aqueous solutions migration through natural porous rocks, the relationship between the porosity and the effective transport properties of geological media still remains difficult to model. Both the aforementioned approaches can be then equivalently applied to reproduce the behavior of the natural system. To overcome all these uncertainties in transport parameters, in our model we chose to explore the sensitivity of numerical outputs with respect to the effective diffusivity, instead of dwelling into second order refinements such as the testing of different relationships (actually non supported by any experimental data) between porosity and transport properties. Different tortuosity coefficients have been then set for the different rock domains and kept constant during the calculations.

Diffusion in the model is expressed with an average diffusion coefficient applied to all dissolved species. This assumption ignores the different tracer diffusivities of each species and electrochemical migration effects due to the different charges of various ions (e.g., Boudreau et al., 2004;Giambalvo et al., 2002).

Diffusion coefficients vary by about one order of magnitude among different species at the temperature of interest of our simulations (75°C). Based on compilations by Oelkers and Helgeson (1988), and Lasaga (1998), we calculate values ranging between 15 and 20×10 -9 m 2 /s, 8.7 and 11.5×10 -9 m 2 /s, 1.2 and 1.5×10 -9 m 2 /s, 3.1 and 4.2×10 -9 m 2 /s, 1.9 and 2.6×10 -9 m 2 /s for the following selected species H + , OH -, Mg and a sensitivity analysis on tortuosity coefficient has been performed. In our calculations we adopted the same value of for both caprock and cement domains, with values in the reservoir larger by one order of magnitude at least, to reflect the different texture and transport properties of local carbonate formations (see, Descostes et al., 2008;Melkior et al., 2004;Sammartino et al., 2003, and references cited therein). Reference values for porosity and tortuosity are given in Table 1.

Actual limitations of the modeling approach

Models presented in this work are an idealization of the physical and chemical processes expected to occur in the cement of a hypothetical abandoned well of the Paris Basin, after CO 2 geological disposal in the local Dogger saline aquifer. The accuracy of the numerical predictions should be improved when more detailed site-specific physical and chemical data will be available, particularly (iv) addressing reactivity at reservoir-caprock-cement interface under two-phase conditions and realistic, site-specific hydraulic conditions, eventually accounting also for the advective transport of chemicals from reservoir; (v) considering the presence of physical discontinuities and the spatial variability of hydraulic parameters, to account for the different geochemical behavior of reactive species expected in fractured and/or heterogeneous rocks from that in homogeneous porous media;

(vi) accounting for possible in-situ ageing effects induced by the interaction of the cement with host formation porewaters before the beginning of the confinement of CO 2 at depth; (vii) improving the analysis of redox conditions, when more accurate information will be available on thermodynamic data and the occurrence of redox-sensitive, Fe-bearing minerals in the cement, and the reactivity of dispersed organic matter in the caprock; (viii) exploring the effect of different, non-linear porosityeffective diffusion relationships, when experimental data will be available on transport properties of site-specific materials.

Despite all these approximations, the numerical simulations revealed some general trends that are important for the assessment of long-term integrity of CO 2 storage.

Results and discussion

Major uncertainties in the calculation of the alteration patterns are related to the effectiveness of diffusive mass transport mechanisms. Uncertainties related to the poor knowledge of kinetic parameters have a negligible impact because of the equilibrium approach used in simulating the behavior of cement minerals. The effective diffusion coefficient is the main parameter governing mass transfer throughout the system. A sensitivity analysis has been performed on its initial value by varying the tortuosity factor (Table 1) for the different rock domains considered in the calculations. Different diffusivity conditions have been explored, and two limiting reference cases have been defined, based on the values of the effective diffusion coefficient, eff d , calculated at reservoir-cement interface. We define a " low-diffusivity" (LD) and a "high-diffusivity" (HD) case simulation, for eff d values equal to 1.65×10 -11 and 1.79×10 -10 m 2 s -1 , respectively. Under these conditions, the effective diffusion coefficient within the cement domain spans the range between 3×10 -13 to 3×10 -11 m 2 s -1 . These values roughly encompass the range of effective diffusivity values reported in the literature for cement materials (e.g. Trotignon et al., 2007). Due to the feedback on porosity, mineral dissolution/precipitation reactions affect mass transfer efficiency, and the effects of this feedback on eff d , as calculated at reservoir-cement and reservoir-caprock interfaces, are plotted in Figure 3. In the following sections, the results for both the LD and HD case simulations are presented. Even larger values of eff d (up to about 10 -10 m 2 s -1 within the cement domain) have been explored by increasing tortuosity values. Under these conditions, reactive fronts propagate over wider areas, and then, larger model geometries need to be implemented. The increase in the number of cells of the domain causes the cost of the simulation to increase, and then the computational efficiency to sharply decrease. For this reason, this even higher diffusivity case is not further investigated or presented here.

Propagation of acid and alkaline fluids

To facilitate the analysis of the propagation of the acid and alkaline fluids through the system, contour plots of aqueous tracer concentration after 1000 years are displayed in Figure 4. Contour lines for tracers give a general picture of the maximum propagation distances the reactive fronts are expected to cover in our system. Numerical results show that after this simulation time, the influence of reservoir fluids in the overlying layers (marked by the presence of tracer concentrations around 1/10 of the initial value in the reservoir) extend up to less than 2 and more than 8 m from the interface, in the LD and HD case, respectively. Due to the small volume of the cement domain compared to caprock and reservoir, cement pore waters have a reduced capacity to penetrate the adjacent domains, and a strong dilution effect is observed within the cement domain, with tracer concentrations lower than 1/10 of the initial values throughout the cement slice in the HD case.

Caprock fluids are predicted to fully penetrate the cement volume in both the LD and HD cases.

Based on these considerations, which reflect the particular geometry and relative volume of different rock domains in our model (Table 1), it is expected that the alkaline pulse is markedly less effective than the acid one in controlling the geochemical evolution of the system. The contour plot of Figure 5, showing the areal variation of pH values after 1000 years, confirm this general picture.

In this figure, a rise of pH values above 7-7.5 in the caprock, due to the migration of the alkaline pulse from the cements filling, is predicted only in the upper part of the model, near cement-caprock interface. In contrast, pH sharply decreases both in the lower part of cement and caprock domains, due to the propagation of the acid fluids coming from the reservoir, and in rightmost slices of cement domain, due to the ingress of near-neutral fluids coming from the caprock.

Mineralogical transformations

In the cement domain

Due to the high buffer capacity of the cement with respect to the acid, C-and Ca-rich fluids diffusing from reservoir and the local equilibrium approach adopted in describing the behavior of cement minerals, most of the mineralogical transformations occur at the bottom of the cement slice and concentrate in the first grid block, near the interface. Mineralogical transformations induced in the cement by the propagation of acid, C-rich fluids from reservoir are shown in Figures 6A and6B with respect to the first block element CEM-01 along the slice #1. Several distinct reactive stages are recognized in this grid cell. The first stage is characterized by the dissolution of portlandite with correlated precipitation of calcite, minor precipitation of hydrotalcite, under almost pH-buffered, alkaline conditions (pH around 11). These processes can be described through the following reactions2 , whose energetic feasibility is verified by inspection of activities of aqueous species involved:

Portlandite + CO The first two stages cover a simulation time of about 250 and 25 years in the LD and HD case, respectively, and cause porosity to decrease down to 0.12±0.02 (stage 1) and 0.03±0.01 (stage 2), respectively.

A third stage occurs roughly between 210 and 370 years (between more than 25 and 38 years in the HD case), which is characterized by the dissolution of primary tobermorite and newly formed straetlingite, and by the precipitation of calcite, Ca-phillipsite, Mg-illite, and Fe(OH) 3 . The following equations can be written: (Rimmelé et al., 2008), has not been reproduced in our simulations due to the high thermodynamic stability of tobermorite (a primary OPC mineral in our calculations) compared to the easily dissolving amorphous CSH phases considered in the laboratory experiments. Quartz predominates in our model (amorphous silica and chalcedony precipitate in negligible amounts) as it is the most stable silica polymorph in our database at the reference temperature of 75°C 3 .

Additional mineralogical transformations include metastable formation of anhydrite (between 200 and 350 years) and chlorite (between 250 and 400 years), and precipitation of celestite and siderite (later than 400 years). All these minerals occur in small amounts and do not significantly affect overall mass transfer in the system.

Summarizing, the most important geochemical reactions resulting in porosity changes in the cement, close to the reservoir interface, are, in a chronological order: i) the early replacement of portlandite and other dissolving primary OPC minerals with calcite, which leads to a significant decrease in porosity, down to -15%; ii) the transition katoite → straetlingite, which causes a slight increase in porosity (up to +2%); iii) the concomitant disappearance of straetlingite and tobermorite and precipitation of Ca-phillipsite and calcite, which determine a further, slight increase in porosity (up to 2%); iv) the late disappearance of Ca-phillipsite and the concomitant precipitation of additional carbonate minerals (magnesite, dawsonite, siderite), which cause, at the end of the reaction path, an additional porosity increase of up to +4%.

Overall, after 1000 years simulation time, cement porosity is lower than in the initial state. This is due to the early carbonation of portlandite, which causes an average increase of the volume occupied by solids in the order of 11.5%. In contrast to some laboratory experiments (e.g. Kutchko et al., 2007;Duguid et al.,2005), the mechanism of late porosity reopening predicted by the code is not related to the dissolution of calcite. This different behavior mirrors, in part at least, the large 3 Despite a number of field and experimental evidences supporting the low temperature stability of quartz in natural environments (e.g. Rimstdit, 1997; White and Brantley, 1995, and references therein), the predominance of different silica polymorphs in natural and engineered systems appears to be a function of not only the temperature, but also of many additional, often ill-defined factors, such as the chemical composition of circulating fluids, the mineralogical composition of the hosting rocks, fluid flow conditions and reaction times (i.e. amorphous and other unstable silica polymorphs tend to transform into cristobalite, opal, chalcedony or quartz over time).

In order to circumvent any possible uncertainty about the predominance of different silica polymorphs at 75°C, in our simulations we then considered the possibility that chalcedony, a silica microcrystalline mineral that can be considered a particular habit of quartz rather than a different phase, could be the only secondary silica polymorph allowed to precipitate over the time frame of interest of our model. Numerical results indicate that under these conditions, no appreciable variations are observed in the mineralogical patterns compared to the reference case simulations (with quartz allowed to precipitate) discussed in this paper.

differences in the initial geochemical conditions (i.e. porewater chemistry and initial mineralogy), and both time and spatial observation scales (months and millimeters vs. hundreds of years and centimeters to meters), between laboratory experiments and our numerical exercise. Laboratory experiments set a severe, upper bound on the degradation rate of cement materials (pH values of the reacting brines as low as 2 to 4; e.g. Duguid and Scherer, 2010), whereas our numerical simulations realistically accounts for the capacity of reservoir rocks to buffer the pH of CO 2 -rich reservoir porewaters to values around 5.

Inspection of the mineral reactivity in the second grid block CEM-02, reveals that minor transformations of the primary mineralogical assemblage are predicted moving away from reservoir interface. In this grid cell, the main processes controlling the geochemical evolution of the medium are still represented by the dissolution of portlandite, the precipitation of calcite and the decrease in pH values, but the computed variations (about -1300 mol m -3 ; +1000 mol m -3 ; -0.5 pH units, respectively) are now lower than in the first grid block (CEM-01). In addition, tobermorite, C 3 FH 6(s) , hydrotalcite and ettringite tend to precipitate instead of dissolving. No additional secondary minerals are predicted to precipitate over 1000 years. In the HD case, the thickness of the heavily altered zone (along the vertical slice #1) extends up to 0.175 m from reservoir interface after 1000 years.

In the reservoir domain

The mutual interaction between cement and reservoir fluids induces very minor variations in the reservoir mineralogical assemblage. The largest variations are observed in the grid cells adjacent to cement interface. Here numerical outputs indicate that the disappearance of accessory magnesite and Mg-illite, and a small increase in calcite, kaolinite and quartz abundances are possible.

Associated to these variations, a minor reduction in porosity, in the order of 1%, is also predicted.

The effects of the alkaline pulse migration from cement into reservoir domain are also revealed by the spatial and temporal variations in chemical composition of reservoir porewaters. However, due to the relatively large volume of reservoir elements compared to cement (Table 1), alkaline fluids are diluted at very early times, and their maximum impact on reservoir geochemistry occurs in the first 100 years simulation (LD case). During this stage, pH rises up to about 5.15, CO 2(aq) concentration decreases to 0.2 mol kg H2O -1 , and aqueous cement tracer rises up to about the 30% of its initial value in the cement, in the RES-01 element, near cement interface. A quite complex behavior is shown by Ca 2+ concentration profile, which reflects the interplay between diffusive migration patterns in the aqueous phase and capture/release effects associated to dissolution/precipitation mechanisms involving Ca-bearing minerals (mainly portlandite and calcite) near the interface. Under HD conditions, the migration of acid fluids from reservoir into cement is maximized, and accordingly, the impact of cement alkaline fluids on reservoir geochemistry further reduced.

At the cement-caprock interface

Water-rock interactions occurring at the interface between caprock and cement, far from reservoir

Porosity variability

The potential impact to the integrity of wellbore cement as a result of exposure to CO 2 -charged natural brines is a primary concern with respect to long-term effectiveness of geological storage of CO 2 . Whereas CO 2 migration in defect-free cement is expected to occur over geological timescales and does not present a major risk, the possibility that some leakage will occur along fractures of the cement plugs should be taken into consideration and carefully analyzed (e.g. Wigand et al., 2009;Carey et al., 2007;Kutchko et al., 2007). Porosity variations predicted by reactive transport codes can be used to evaluate if the geochemical processes expected in the simulated systems may affect isolation properties of cement, providing high permeability rising pathways for CO 2 . The basic assumption is that any computed decrease in porosity can be thought to be beneficial in terms of containment capacity of the storage reservoir, whereas any increase in porosity can be considered of some concern with respect to the long-term potential for undesired fluids leakage. This simplified approach presents a major drawback in that it does not account for the changes in petrophysical and geomechanical properties of rocks and cement, likely associated with the largest mineralogical transformations predicted by the reactive transport code. Field and laboratory evidences, and modeling results (e.g. Taylor, 1990;Barlet-Gouéadard et al., 2007;Carey et al., 2007;Duguid, 2008;Rimmelé et al., 2008;Fabbri et al., 2009;Wigand et al., 2009;among during the injection of supercritical CO 2 , and later on, during the in-situ ageing of materials, may in fact develop fractures and preferential rising pathways for reservoir fluids. As of today, the lack of fully coupled geomechanical-reactive transport codes, prevents any quantitative assessment of these contrasting effects. Two possible scenarios can be then envisaged: i) in the best possible conditions, massive cement carbonation can be considered as a beneficial effect of water-rock interactions, which ensures site containment capability; ii) in a less favorable scenario, massive cement carbonation leads to a degradation of the desirable properties of cement plugs. Reactive transport modeling can be used to recognize as zones of potential concern all the portions of the modeled system where large variations in mineralogical composition and/or porosity are predicted.

In our model, the largest mineralogical transformations occur at the caprock-cement interface, in close proximity to the reservoir, where the interaction between the three chemical systems is maximum. For instance, portlandite preferentially dissolves along the vertical slice #4, showing an appreciable decrease in the volume fraction over a distance of about 0.5 m from reservoir interface after 1000 years (LD case; Fig. 7). In the HD case, a complete dissolution of portlandite is predicted all over the lower portion of the cement domain, up to a distance of about 0.175 m from reservoir interface, together with a decrease in volume fraction from 0.328 to about 0.24 in the upper part of the vertical slice #4, up to the top of the modeled system. Carbonation is also predicted already after a few years of cement exposure to external fluids in this part of the system.

Under the diffusion-controlled conditions of these calculations, the thickness of the calcium carbonate layer formed in the cement at reservoir interface increases progressively from the vertical slice #1 (between less than 0.1 to about 0.25 m after 1000 years in the LD and HD case, respectively), towards the vertical slice #4 (up to more than 1 m in the HD case; Fig. 8).

Abundant crystallization of carbonate minerals is also predicted in the caprock, at the cement interface, along the vertical slice #5, as a result of the coupling between the three different zones.

Porosity variations reflect all the mineralogical transformations described above. Again, patterns differ significantly along the three vertical slices used as reference. Along the slice #1, the buffer capacity of cement minerals prevents the alteration front from propagating more than a few centimeters from reservoir interface, resulting in a very steep pattern of porosity decrease (up to -8 to -14% in the LD and HD case, respectively), whereas a more effective decrease in porosity over wider portions of the system is predicted along the slices #4 and #5. Overall, the greatest effects of the geochemical perturbation are predicted along the slice #4 (Fig. 9), where, under HD conditions, an increase in porosity (up to +8%) is observed far from reservoir interface, next to the portion of the cement domain characterized by the largest porosity decrease (up to -16% in the HD case). The decrease in porosity in the lower part of the slice #4 could be put in relation with the migration of acid fluids from reservoir, whereas the increase in porosity along the upper part of the same slice reflects the penetration of fluids from the adjacent caprock and the concomitant leaching of species from cement into caprock. The alteration of the cement with the surrounding host shale formation is relevant for the assessment of the sealing integrity of the repository, because it might lead to the formation of a vertical preferential ascent route towards the surface, and allow reservoir fluids to leak over larger simulation times than considered here. Of course, this pessimistic scenario is strictly related to the possibility that the enhanced carbonation of OPC predicted at the bottom of the cement domain, may develop a set of interconnected fractures connecting the reservoir with the zone of enhanced porosity along the cement-caprock interface.

In the lowermost part of the caprock domain, near the reservoir interface (distance < 1m), minor variations in porosity are predicted, but significant mineralogical transformations are still expected, consisting of calcite, illite and siderite precipitation, and chlorite dissolution.

A sensitivity analysis on calcite dissolution/precipitation behavior indicates that the numerical model is quite robust (Fig. 10). Calcite was selected because carbonation is the most relevant process predicted at cement edges, which controls the spatial and temporal evolution of water-rock geochemistry and porosity. Again, porosity may be used as a master variable to trace the overall geochemical effects induced under equilibrium and kinetically controlled conditions. As shown in Figure 10, almost negligible differences can by recognized in terms of porosity variation when switching from kinetic to equilibrium dissolution/precipitation conditions for calcite, compared with variations related to low-(LD) and high-diffusivity (HD) end-member conditions.

Additional sensitivity analysis

On caprock porewater composition

The direct measurements of Michard and Bastide, 1988;Coudrain-Ribstein and Gouze, 1993) for a temperature close to our simulation temperature. Consequently, the chemical composition of a second, relatively "high-P CO2 " COx pore water ("Caprock #2" water of Table 5), compatible with the mineralogical assemblage shown in Table 3, has been then computed using batch numerical modeling, and used in 2-D simulations for sensitivity analysis. A reference 2 log CO f value of -0.98 has been obtained.

Despite the lower initial pH value (6.33), and the higher C TOT contents of this water (3 × 10 -3 mol kg H2O -1 ), very minor differences in the mineralogical assemblage alteration have been predicted by the code with respect to the reference case already discussed. In particular, minor porosity changes at the caprock-cement interface are computed after 1000 years (up to -0.4% and -1.2% in the cement and caprock, respectively), which are prevalently driven by calcite precipitation in the caprock (up to 318 mol m -3 -medium at the interface) and portlandite dissolution in the cement (up to -310 mol m -3 -medium). Near the caprock interface, cement primary mineralogy is also affected by katoite dissolution (up to -15 mol m -3 -medium) and precipitation of traces of ettringite, tobermorite, hydrotalcite and C 3 FH 6 (always lower than +10 mol m -3 -medium), whereas the most significant modifications in argillite mineralogy are Mg-illite and quartz precipitation (up to +25 mol m -3 -medium), and Mg-montmorillonite-Na dissolution (up to -35 mol m -3 -medium).

On cement mineralogy and porewater composition

The effects of different portlandite and tobermorite relative amounts in the initial idealized Portland-type cement have been tested by introducing into the calculations the mineralogical composition and associated pore water chemical composition ("Cement #2") listed in Table 5. The new configuration of the system (Caprock #2 + Cement #2 initial conditions) leads to more conservative results compared to the scenarios already discussed. The most remarkable differences concern the alteration products predicted in the element CEM-01 of the cement, near the reservoir interface. Here, larger amounts of precipitated silica minerals (i.e. quartz in the reference case simulation, chalcedony when quartz is not allowed to precipitate), and almost no precipitation of Mg-illite are now computed, together with slightly lower amounts of precipitated magnesite and

Conclusion

A complex 2D diffusive transport numerical model has been developed to assess the geochemical integrity of a well cement filling connected to two separate mineralogical assemblages representative of both reservoir and caprock porous rock. Cement alterations of abandoned wells may play an important role in guideline definition and selection criteriaof carbon dioxide geological storage sites. The consideration of three different mineralogical domains induces a very complex fluid rock geochemical interaction system which summarizes as the following major conclusions. (c) Geochemical interactions at the top of the reservoir generally induce some pore clogging in the lowermost part of both cement and caprock domains, primarily due to enhanced precipitation of secondary calcite. In contrast, porosity is predicted to increase in the cement, in the upper part of the system, along the caprock interface. Close to the reservoir, along the cement-caprock interface, porosity initially decreases down to about 3% (after 200 years), and then increases up to about 10% (after 750 years). Overall, the mineralogical transformations predicted at the bottom of the system, after the penetration of acid fluids from reservoir, and above, the increase in porosity caused by caprock-cement interactions, might lead to the formation of a vertical preferential ascent route, along the external edge of the cement filling, for reservoir fluids. This scenario is related to the possibility that a set of interconnected fractures may develop in the cement, as a consequence of the enhanced carbonation of the OPC plug near the reservoir interface.

FIGURES

(

  >50°C) temperature reported in the literature(Paul and Glasser, 2000), namely: (i) CSH amorphous phases are metastable with respect to crystalline CSH phases such as tobermorite, and even if they had been found to be only partially crystallized in laboratory experiments and field conditions(Glasser et al., 2005), they are expected to convert to crystalline phases over the reference time frame of our calculations; (ii) the prolonged hydration of cements favors the formation of siliceous hydrogarnets and (iii) hydrotalcite, whereas (iv) monosulphate or ettringite are expected to persist, depending on sulfate and aluminium concentrations in the aqueous phase;(v) katoite-Si and (vi) C 3 FH 6 (s) tend to replace AFm (monocarboaluminate and C 4 FH 13 (s)) and AFt (ettringite) phases. Based on recent data fromBlanc et al. (2010a and 2010b), afwillite emerges as an additional crystalline CSH phase to be considered in the numerical reconstruction of an hypothetical equilibrium OPC mineralogical assemblage at 75°C. However, the limits of the stability field of afwillite are still rather uncertain, as suggested by the discrepancies between the data recently gathered byBlanc et al. (2010a and 2010b) and the thermodynamic properties determined byDickson et al. (2004) and/or estimated byBabushkin et al. (1985). Moreover,

  Calcium carbonates are represented by calcite only. Aragonite and vaterite polymorphs are not taken into account being thermodynamically unstable at the reference temperature of 75°C.Additional carbonates are dolomite, ankerite, magnesite, siderite and dawsonite. The clay fraction of the Callovian-Oxfordian (COx) argillite is dominated by mixed-layer illite/smectite, Mg-illite, and includes minor amounts of Mg-Na-montmorillonite, but several additional types of illites, montmorillonites and saponites have been inserted in the list of potential secondary minerals. The Mg-Na-bearing montmorillonite has been also used as thermodynamic proxy for the mixed-layered fraction of illite-smectite minerals. Quartz is a major phase of the COx formation and has been considered as the primary silica polymorph in the calculations, but chalcedony and amorphous silica are allowed to precipitate as possible secondary phases. Chlorite group minerals are represented by means of two single minerals corresponding to the Mg-and Fe-enriched end-members, clinochlore and daphnite. Sulfur-bearing minerals consist of anhydrite, pyrite and celestine, whereas analcime, and three different phillipsites have been used to represent zeolite group minerals. Non-clay silicates such as muscovite, K-feldspar and low-albite have been considered in the calculations.

  +2 , CO 2(aq) , HCO 3 -, respectively, depending on the value of the activation energy (between 10 to 15 kJ/mol) used in the Arrhenius-type equation giving the dependency of the diffusion coefficient on temperature (Lasaga, 1998): for the effective diffusion coefficient, eff d , reasonably comparable with values from field and laboratory experiments (De Vera et al., 2007), the largest value of this range, corresponding to proton bulk diffusion coefficient in water at infinite dilution, H T d , has been used,

  on mineralogy and porewater composition of cement materials. The reliability of numerical predictions is also affected by the severe limitations in the current knowledge of kinetic parameters of minerals, especially of both crystalline and amorphous solid constituents of cement. In this study, all uncertainties have been tackled through simplification of conceptual and geometric models, and by carrying out a sensitivity analysis on relevant features. Future work could improve on all simplifications made by (i) considering species-dependent variations in the aqueous phase diffusion coefficients; (ii) considering the presence of amorphous and/or not completely hydrated components in the starting cementitious material and addressing the simulation of their alteration mechanisms by means of a fully kinetic treatment of precipitation/dissolution reactions; (iii) following a different numerical approach, incorporating solid solutions into equations describing reactive transport, to describe the CSH phase evolution, and to reproduce the experimental observations of progressive CSH decalcification with decreasing pH(Carey and Lichtner, 2007);

Fe

  2+ + 2 OH -+ H 2 O → Fe(OH) 3(s) + 0.5 H 2(aq) (eqn. 11), characteristic of a very short time domain around 220 years simulation time, and: Straetlingite + CO 2(aq) + 4 SiO 2(aq) → Calcite + Ca-phillipsite + 2.5 H 2 O (eqn. 12). After about 370 years (about 38 years in the HD case), the code predicts the complete disappearance of tobermorite, of secondary straetlingite and clinochlore, the precipitation of magnesite, and a drop of pH to 4.8. The Mg released in aqueous solution by Mg-chlorite dissolution, is captured by precipitating Mg-illite and magnesite. After the complete alteration of the cement initial mineralogy, the fluids coming from reservoir exert a total control on the geochemistry of the medium. This trend is further enhanced in the successive 600 years (LD case), as shown by the two additional zones predicted by the code. The first is characterized by the dissolution of Ca-phillipsite and the precipitation of Mg-illite, whereas the second one, roughly occurring between 600 and 800 years (LD case), is characterized by the precipitation of kaolinite, dawsonite and quartz. Porosity slightly increases with time up to about 0.10 after 1000 years. The late precipitation of quartz predicted by the numerical model is qualitatively in agreement with the end result of the degradation that hardened cement underwent upon carbonic acid exposure in the laboratory under similar (50°C) conditions (Kutchko et al., 2007). The early precipitation of amorphous silica observed during some laboratory experiments on Portland cement exposed to CO 2 -rich fluids for several weeks (up to six months) and at temperatures between 50 to 90°C

  interface, cause porosity to slightly increase in the cement and to decrease (up to -2%) in the caprock. The diffusion of H + component from caprock to cement allows pH to decrease in the cement (buffered at about 10.8 after few tens of years), and to progressively increase in the caprock (up to about 9.4 after 1000 years). Apart from significant portlandite dissolution, minor mineralogical transformations are predicted in the cement. Ettringite, tobermorite, hydrotalcite and calcite precipitate, and katoite dissolve in very low amounts (delta values up to a few moles m -3 -MEDIUM ). Apart from calcite, no secondary minerals form in the cement.The mineralogical transformations observed in the caprock, far from reservoir interface (CAP-03 grid block in our monitoring system; Fig.2), are coherent with the alteration patterns already described in the literature in association with high-alkaline fluids circulating through argillaceous matrix (Gaucher and Blanc, 2006; and references cited therein). The most relevant processes predicted in this part of the system are the dissolution of montmorillonite, quartz and dolomite, and the precipitation of Ca-zeolites (Ca-phillipsite), illite, chlorite and calcite. The occurrence of zeolites at the cement/clay interface is well documented over a large temperature range, up to about 150°C. Precipitation of phillipsite and analcime is documented by Bauer and Berger (1998), Chermak (1993), de la Villa et al. (2001), Mosser-Ruck and Cathelineau, (2004), Ramirez et al. (2002, 2005), Sanchez et al. (2006). Chabazite has been observed by Ramirez et al. (2005) duringalkaline alteration of COx. Field data and experimental evidences then suggest that zeolites could be thermodynamically stable at the temperature of interest of our simulations (i.e. 75°C), and that their precipitation may exert some control on aqueous Ca concentration of caprock and cement porewaters. No CSH (amorphous and crystalline) phase precipitation is predicted, because of the limited increase of Ca/Si ratio in aqueous solution which does not allow these minerals to reach the saturation state over the time frame of 1000 years. Montmorillonite dissolution is the most significant mineralogical transformation predicted in the caprock, close to cement interface. Similar to many other clays, montmorillonite swells with the addition of water, acting as a natural plug against fluid leakage. Its decrease in concentration may then negatively affect the sealing integrity of the repository all along the annulus between the cement and the caprock. In the next forthcoming years, the advancement of focused fundamental research and the development of fully coupled geomechanical-reactive transport codes will likely help better clarify on this issue currently not addressed by conventional reactive transport codes. Additional mineralogical transformations include the transitory precipitation in traces of saponites and phillipsites (in the first 600 years), and minor dissolution of celestite and pyrite.

  several others) indicate in fact that, whereas moderate carbonation can be beneficial to cement strength, the formation of cracks is usually recognized in association with extensive carbonation processes. For instance, Barlet-Gouéadard et al. (2007) reported a strength loss of up to 65% for heavily carbonated Portland cements exposed to CO 2 -saturated fluids for six weeks, and the occurrence of cracks in the carbonation layer during compressive strength measurements. Fabbri et al. (2009) observed a significant degradation of mechanical properties of carbonated cement under wet supercritical CO 2 conditions at 90°C, associated with the presence of micro-cracks in the carbonated layers of the Portland cement. Although the deposition of calcium carbonate in the pores may locally allow the porosity to decrease, massive carbonation episodes within the OPC should be then considered not only as mechanisms of possible further sealing of the cement, but also as possible sources of late mechanical degradation of the material. Stresses and potential differential pressures likely occurring

  this formation. However, as different mineralogical assemblages (with different thermodynamic values) are possible for this buffering (Gaucher et al, 2009), it is difficult to reduce the range of uncertainties for this parameter under a log unit. Moreover, we can recall the 2 log CO f values obtained for the underlying Dogger aquifer ( 2 log CO f roughly between -1.5 to -0.4,

  dawsonite. MonoSA-Fe and anhydrite appears among the secondary alteration products during a transient phase comprised between about 200 and 280 years and between about 280 and 780 years, respectively, in the low diffusivity case. Relevant amounts (up to about 1075 mol m -3 -medium) of newly precipitated kaolinite are also predicted after about 800 years simulation time. Again, porosity variations are used to trace the overall effect of all the occurring mineralogical transformations: the larger variations are predicted to occur along the slice #4, near caprock interface, with values between -11.5 to -4.5%, after about 200 and 1000 years simulation time, respectively (Figure 11 EA ).

  (a) Numerical simulations predict a two-stage evolution of the cement porous matrix, after interaction with acid fluids from reservoir: (i) a first, "clogging" stage, characterized by a decrease in porosity after precipitation of calcium carbonate (calcite in our model), and (ii) a second stage of pore volume reopening, related to a decrease in calcite precipitation, the complete disappearance of primary cement phases and the re-dissolution of large-molar volume secondary minerals, such as Ca-zeolites. The penetration of these alteration fronts is limited to the first few to few tens of centimeters from reservoir interface (up to a maximum of about 25 cm after 1000 yrs), depending on the efficiency of mass transfer mechanisms, i.e. on variations of aqueous diffusion coefficient, porosity and tortuosity. (b) This mechanism of porosity reopening differs from experimental observations made by Kutchko et al. (2007) and Duguid et al. (2005), who related the porosity reopening to calcite dissolution. In this work, calcite is not predicted to dissolve because the numerical model realistically accounts for the capacity of reservoir rocks to buffer the pH of CO 2 -rich reservoir porewaters to values around 5. In order to reproduce the short-term re-dissolution of calcite observed during laboratory experiments, less realistic geochemical boundary conditions, i.e. lower pH fluids with timeinvariant chemical composition in the reservoir, should be adopted. This has already been shown for a Class G high Sulfate Resistent commercial oilwell cement exposed to CO 2 -rich aqueous solutions over a period of 100 years and at a constant temperature of 45°C (Biagi and Gherardi, 2010). Initial and boundary conditions of these calculations set an upper limit to the degradation rate of cement materials, and closer represent laboratory experiments already published in the literature (e.g. Duguid and Scherer, 2010).

  Table 1 of the Electronic Annex, thereinafter named Table 1 EA . The

	Equilibrium constants for portlandite, brucite, magnesite, gibbsite, and Friedel's salt are based on
	data from Robie and Hemingway (1995), Chase (1998) and Akinfiev and Zotov (1999), whereas
	those for Fe-bearing minerals (C 4 FH 13 , C 3 FH 6 , ettringite-Fe, Monosulfate-Fe, ferrihydrite) have
	been calculated from Babuskin et al. (1985) and Chivot (2004). Thermodynamic properties of
	minerals belonging to CaO-Al 2 O 3 -SiO 2 -H 2 O system (katoite-high-SiO 2 , straetlingite) were
	calculated on the basis of data from Atkins et al.

equilibrium constants are taken from the Thermoddem database (February 2008 release), a thermodynamic database compiled at BRGM (http://thermoddem.brgm.fr/). Thermoddem uses internally consistent and critically evaluated thermodynamic data derived from a comprehensive review of the published literature. The primary source for equilibrium constants for aqueous species and minerals is the Slop98.dat database

(Shock et al., 1998)

, which is in turn a substantial revision of the SUPCRT92 database

(Johnson et al., 1992)

. A major modification with respect to the original Slop98 data set concerns data for Al aqueous species, which have been taken from the experimentally-derived compilation of

Tagirov and Schott (2001)

. Data for minerals have been thoroughly revised and major changes concern cement minerals

(Blanc et al., 2006)

. Solubility constants for CSH (amorphous and crystalline) phases have been calculated on the basis of thermodynamic data and experimental evidences from

Babuskin et al. (1985)

,

Hong and Glasser (2004), Shaw et al. (2000), and Glasser et al. (2005)

. Thermodynamic properties of minerals belonging to CaO-Al 2 O 3 -SO 3 -H 2 O system (ettringite, monosulfate, hydrogarnet) have been extracted from

Warren and Reardon (1994), Perkins and Palmer (1999), and Babuskin et al. (1985)

.

  This stage is characterized by the dissolution of all primary cement minerals (tobermorite, katoite-Si, C 3 FH 6(s) , ettringite, hydrotalcite), and the precipitation of secondary straetlingite, calcite, and clinochlore. During this phase, alkaline conditions still persist, but some decrease in pH is observed (pH around 10.2), because of the partial buffering of H + ions coming from the adjacent reservoir. The overall process can be tentatively represented in terms of the following elementary processes:

	Tobermorite + 5 CO 2(aq) → 5 Calcite + 10.5 H 2 O + 6 SiO 2(aq)		(eqn. 6)
	Katoite + CO 2(aq) + 3.5 H 2 O → Calcite + Straetlingite		(eqn. 7)
	Hydrotalcite + 0.4 OH -+ 2.4 SiO 2(aq) → 0.8 Clinochlore + 0.4 AlO 2	-+ 7 H 2 O	(eqn. 8)
	Ettringite + 6 CO 2(aq) + 8 OH -→ 6 Calcite + 2 AlO 2	-+ 36 H 2 O + 3 SO 4	2-	(eqn. 9)
	C 3 FH 6(s) + 4 H + + 3 CO 2(aq) + H 2(aq) → 3 Calcite + 9 H 2 O + 2 Fe 2+		(eqn. 10)
	2(aq) → Calcite + H 2 O			(eqn. 4)
	4 Mg 2+ + 6 OH -+ 2 AlO 2	-+ 7 H 2 O → Hydrotalcite			(eqn. 5)
	and with Al TOT in the aqueous solution supplied by dissolving ettringite.

The onset of the second reactive stage, which occurs after about 180 and 20 years in the LD and HD cases, respectively, is concomitant with the disappearance of portlandite.

Reactions are written in terms of predominant species and not basis species, i.e. CO 2(aq) and OH -are the predominant C-and H-species in the pH-acid and pH-basic fluids of reservoir and cement domains, respectively

Acknowledgements

As part of Carnot Labeled Research Institutes, BRGM received fundings to elaborate this work and to collaborate with the Istituto di Geoscienze e Georisorse (Pisa, Italy). The French National Carnot Program aims at financing exchange expertise and research development between public and private research entities. The authors thank Philippe Blanc, Francis Claret, André Burnol, Nicolas Jacquemet, Nicolas Marty (all at BRGM) for fruitful discussions. Comments by the Editor L.Charlet, the Associate Editor T.Kozaki, two anonymous reviewers and Ian C.Bourg largely contributed to improve the quality of the manuscript.