Numerical modeling of carbon dioxide sequestration on the rate of pressure solution creep in limestone: Preliminary results

When carbon dioxide (CO2) is injected into an aquifer or a depleted geological reservoir, its dissolution into solution results in acidification of the pore waters. As a consequence, the pore waters become more reactive, which leads to enhanced dissolution-precipitation processes and a modification of the mechanical and hydrological properties of the rock. This effect is especially important for limestones given that the solubility and reactivity of carbonates is strongly dependent on pH and the partial pressure of CO2. The main mechanism that couples dissolution, precipitation and rock matrix deformation is commonly referred to as intergranular pressure solution creep (IPS) or pervasive pressure solution creep (PSC). This process involves dissolution at intergranular grain contacts subject to elevated stress, diffusion of dissolved material in an intergranular fluid, and precipitation in pore spaces subject to lower stress. This leads to an overall and pervasive reduction in porosity due to both grain indentation and precipitation in pore spaces. The percolation of CO2-rich fluids may influence on-going compaction due to pressure solution and can therefore potentially affect the reservoir and its long-term CO2 storage capacity. We aim at quantifying this effect by using a 2D numerical model to study the coupling between dissolution-precipitation processes, local mass transfer, and deformation of the rock over long time scales. We show that high partial pressures of dissolved CO2 (up to 30 MPa) significantly increase the rates of compaction by a factor of ~ 50 to ~ 75, and also result in a concomitant decrease in the viscosity of the rock matrix.


Introduction
The subsurface sequestration of CO 2 in geological repositories is frequently cited as a promising solution for reducing the amount of anthropogenically-produced CO 2 in the atmosphere. Some of the important issues involved in the long-term sequestration of CO 2 in such sites are discussed in an overview by Wawersik et al. (2001), for example.
Therein it has been judged essential that models need to predict CO 2 sequestration behavior over time periods of several thousand years, which is the same order of magnitude as some climatic cycles. Therefore, in order to advance our knowledge of processes involved in the geological sequestration of CO 2 , one of the most important objectives confronting geoscientists is understanding and quantifying all of the mechanochemical processes, at both short and long time scales, that are relevant to CO 2 storage in geological formations.
When considering O H CO 2 2 − injection into a geological repository, the following general chemical reactions (Stumm and Morgan, 1996)  (2) Since we examine CO 2 sequestration within the context of pervasive pressure solution creep (PSC) in limestone, Eq. (2) is thus of particular relevance to this study.
PSC is a mechano-chemical process characterized by ductile deformation and local mass transfer affecting water saturated porous rocks (e.g. Weyl, 1959;Rutter, 1976;Gratier and Guiguet, 1986;Dewers and Ortoleva, 1990;Spiers and Brzesowsky, 1993;Gundersen et al., 2002;Yasuhara et al., 2003). This ductile deformation mechanism occurs in the upper crust and plays an important role in the compaction of sedimentary rocks during diagenesis (Ortoleva, 1994;Tuncay et al., 2000). PSC is driven by differences in chemical potential induced by differential stress along grain surfaces in the rock matrix. PSC can be modeled as a serial process involving four successive steps: -stress-enhanced dissolution at grain-grain interfaces subject to elevated normal stress; -diffusion of dissolved material (solutes) through intergranular fluid films; -precipitation of dissolved material in adjacent pore spaces (grain surfaces in contact with pore fluid); -transport of dissolved material to distant pores, which can induce local mass transfer (Gundersen et al., 2002).
Since it is assumed that PSC operates as a serial process, the slowest step imposes the overall rate for deformation (Rutter, 1976;Gratz, 1991). The first three individual steps of PSC (dissolution, diffusion, precipitation) are in turn influenced by local parameters such as temperature, stress state of the rock matrix, fluid pressure, and fluid chemistry (Rutter, 1976). The PSC mechanism of deformation is slow and operates over long geological time scales. Because of this, it is even possible that the slow step can change from one process to another over time during compaction of the rock matrix.
Because PSC operates over geological time scales, correctly predicting the longterm stability of a CO 2 repository requires accurate modeling. PSC models are based for the most part on kinetic and equilibrium parameters derived from laboratory dissolution and precipitation experiments, as well as from pressure solution experiments that typically run for only a small fraction of the time scales associated with natural PSC deformation.
Injection of CO 2 causes chemical and flow regime perturbations that affect the PSC process, causing the porosity, permeability, and mechanical stability of the porous rock matrix to evolve over time. There are multiple reasons that are responsible for this.
The acidification of pore fluids due to the dissolution of CO 2 generally increases rates of fluid-rock interactions. This is particularly important in limestones where a decrease in pH increases both the rate of calcite dissolution and calcite solubility. In addition, higher concentrations of dissolved calcium carbonate can result in increased rates of precipitation. The porosity of the rock matrix is reduced during PSC, due primarily to grain indentation and precipitation in pore spaces. Taken together, the dissolution, diffusion, and precipitation processes have the potential for modifying the long-term porosity and permeability of the repository rock, as well as its mechanical stability.
The deformation of chalk, both dry and in the presence of fluids, has been widely investigated experimentally. Most of the published work is based on chalk deformation data that have been interpreted solely in terms of mechanical processes (e.g., Botter, 1985;Da Silva et al., 1985;Jones and Leddra, 1989;Monjoie et al., 1990;Shao et al., 1994;Piau and Maury, 1995;Schroeder and Shao, 1996;Homand et al., 1998;Risnes and Flaageng, 1999). Recently, however, a few studies have examined chalk-fluid deformation within the context of chemical processes associated with PSC (Hellmann et al., 2002a, b;Heggheim et al., in press). In addition, PSC rates have also been measured for calcite aggregates (Zhang et al., 2002) and single calcite grains (Zubtsov et al., 2004).
In this study, we examine the effect of elevated concentrations of dissolved CO 2 on the overall PSC rate of limestone dominated aquifers or reservoirs at burial depths relevant to CO 2 storage (1000-3000 m). The model treats the post-injection phase of sequestration, where the aqueous pore fluids have been homogeneously acidified by the presence of CO 2 . The partial pressure of CO 2 is fixed and remains constant for each simulation. The model also makes the approximation that the 2 CO p equals the pore pressure in the reservoir. Using a 2D numerical model, we examine how various parameters such as grain size, burial depth, rock texture, and the partial pressure of CO 2 modify the rate of matrix compaction by pervasive PSC. Only the effect of CO 2 dissolved in water is considered (i.e. a single subcritical aqueous phase); we do not consider the injection of a supercritical CO 2 phase since the reactivity (i.e. solubility) of calcite therein is predicted to be minimal.
Below, we first present a brief review of the relevant thermodynamics and kinetics of the calcite-H 2 O-CO 2 system that is used to model the dissolution-precipitation processes. This is followed by a description of the basics of our PSC model. Lastly, we present several results from the 2D simulations.

Conventions
Before we address the thermodynamics of the calcite-H 2 O-CO 2 system, it is important to note the conventions we use. First, even though equilibrium constants are by definition based on activities of product and reactant species and phases, we make a simplifying assumption of equating activities of aqueous species with concentrations (in mol m -3 ); that is we assume that their activity coefficients are equal to 1. This approximation is based on the relatively low ionic strength of the solutions (the maximum ionic strength (I) = 0.115 molal, calculated using EQ3NR, Wolery, 1992, and the activity coefficient is greater than 0.85, see Fig. 1 of Kervevan et al., 2005). Solid phases are always assigned an activity equal to 1. We also assume that the molarities are equal to the molalities and the density correction is neglected. Lastly, the use of the term 'fluid' has the same meaning as 'solution', no supercriticality is implied.

Overview of equilibrium thermodynamic relations at 25 °C and 0.1MPa
In this study we consider the system CaCO 3 -H 2 O-CO 2 with the following equilibrium constants (Nordstrom et al., 1990) 6.5, Langmuir, 1997). The only exception to this is the species − 2 3 CO whose concentration, albeit low, is necessary for the calculation of the equilibrium constant sp K .
In addition to the above equilibria, the following aqueous charge CO , and OHare of minor importance and can be neglected in this charge balance.
Using the chemical equilibria and charge balance relations given in Eqs. 3a-i, the concentrations of all chemical species pertinent to the calcite-H 2 O-CO 2 system at STP can be calculated. However, in order to be useful for PSC, the code recalculates these equilibrium concentrations to conform to the relevant temperature and pressure conditions of the fluids present either in the contact zone or in the pore space (as discussed in the following sections). Note, however, that the charge balance relation in Eq. 3i is not pressure or temperature dependent.
In order to simplify the numerical code, only the + 2 Ca concentration is allowed to vary as a function of time. All of the other species' concentrations remain fixed at their initial equilibrium values calculated for t = 0. Thus, at this stage in the development of the code, we have made the simplifying assumption that PSC depends only on one chemical species, such that only + 2 Ca dissolution, precipitation, and local transport are rate-determining.

1.3.
Equilibrium thermodynamic relations at elevated temperatures and pressures 1.3.1. Dependence on temperature An empirical expression for the solubility of calcite (Eq. 3f) as a function of temperature T and salinity S has been derived by Mucci (1983) and is given by where T is in Kelvin. As defined by Mucci (1983), this equation is valid for 0-40°C and S = 5-35 g/kg. For the purposes of this study, we restrict ourselves to the first four terms of Eq. 4, which is the expression originally determined by Plummer and Busenberg (1982) and is applicable at 0-90°C.
The values of a, b, c, d, and e in Eq. (5) are given in Table 2. 1.3.2. Dependence on pressure of the equilibrium constant of calcite dissociation The expression for the pressure dependence of K sp (Eq. 3f) that we use is the following (Lown et al., 1968;Millero, 1982) ln K sp Here, ΔV r 0 and Δκ r are the molal volume and compressibility changes for the dissolution reaction, P is the applied pressure in MPa, and R = 8.32 cm 3 MPa mol -1 K -1 . Thus, using 0 Ksp V Δ = -58.3 cm 3 mol -1 at 25°C (Langmuir, 1997) and Δκ r = -1.5·10 -3 cm 3 MPa -1 mol -1 at 25°C for the solubility of calcite in pure water (Owen and Brinkley, 1941), we can then calculate K sp at elevated pressures. As an example, a change in system pressure from P 0 =0.1 to P=10 MPa results in the following increase in the solubility product: 26 .
In the present version of the code, the pressure dependencies of the other equilbrium constants (Eqs. 3a-e, g-h) have not been considered.

Overview of kinetic rate laws for calcite dissolution
The rates of calcite dissolution and precipitation in aqueous solutions have been actively investigated over many decades (e.g. see a thorough review by Morse and Arvidson, 2002; and references therein). One of the most widely used empiricallyderived kinetic rate laws is based on the work of Busenberg, Plummer, and co-workers (Plummer et al., 1978;Busenberg and Plummer, 1986), in which the rate of reaction (dissolution or precipitation) of calcite can be expressed as a function of three parallel forward rates and one backward rate where R is the overall reaction rate, k 1 , k 2 , k 3 , k 4 are kinetic rate constants, and [i] represents the bulk solution concentrations (in mol/m 3 of water) of aqueous species i.
When considering only the three forward reactions (k 1 , k 2 , k 3 ) in Eq. (7) space there is also a region where all three terms must be considered (see Fig. 12, Plummer et al., 1978).
In circum-neutral pH solutions, and close to calcite equilibrium ( Ω > 0.6), the dissolution reaction can be represented by a simpler alternative kinetic rate law for dissolution that has the following form (Berner and Morse, 1974;Busenberg and Plummer, 1986;Wollast, 1990;Hales and Emerson, 1997) n sp diss n diss Here k diss is the rate constant for calcite dissolution (denoted by k 5 in Busenberg and Plummer, 1986), is the chemical reaction affinity term (i.e. degree of solution undersaturation or oversaturation), n is the reaction order, and Ω is defined as the ratio of the ion activity product (see Section 1.1) to the equilibrium constant K sp (Eq. 3f). For calcite, n = 1 and log k diss = -9.93 at 25 °C and 0.1 MPa (Plummer et al., 1978, units of k diss in mol cm -2 s -1 ). A linear relation between the rate of dissolution and chemical affinity (i.e. n ≈ 1) at conditions close to equilibrium has been reported in other studies, as well (e.g. Cubillas et al., 2004). The rate law given in Eq. (8) has in fact been described as, "the most commonly used equation in geosciences to describe the rate of carbonate mineral dissolution" (Morse and Arvidson, 2002). The kinetic expression in Eq. (8), which describes the net rate of reaction (dissolution and precipitation), has both a mechanistic and empirical origin (Wollast, 1990), based on the However, the simple, quasi-linear relationship between rate and chemical affinity shown in Eq. 8 is perhaps still debatable. As an example, Svensson and Dreybrodt (1992) report that for natural (i.e. impure) calcite, the dissolution rate-affinity relation deviates from a linear relation chose to equilibrium, yielding a value of n ≈ 3-4 in Eq. 8. In the following treatment, however, we neglect the uncertainty associated with this potential non-linearity, as it is much smaller than the uncertainty associated with the measurement of the value of k diss that we adopt.

Effect of temperature and pressure on the rate constant k diss
Elevated T and P conditions have two effects on the overall rate of calcite dissolution: they modify the rate constant k diss , and also the chemical affinity term . The dependence of k diss on T can be expressed using the classical Arrhenius where A is a pre-exponential frequency factor, and E a is the overall thermal activation energy. Plummer et al. (1978) determined an E a = 33.05 kJ mol -1 for the rate constant k 3 (Eq. 7); this corresponds to A = 11.7. An E a value of 35 kJ mol -1 was also reported by Sjöberg (1978). Morse and Arvidson (2002) proposed that activation energy values in this range are consistent with surface reaction control kinetics, but the value of 35 kJ/mol suggests a mixed control: activation energies in the range 0-25 kJ/mol are typical for diffusion in a fluid media, whereas 50 kJ/mol and greater is typical for surface reaction control (Lasaga, 1998). In the study of Busenberg and Plummer (1986), the values of k 3 (Eq. 7 above) and k diss (Eq. 8 above, equiv. to k 5 ) are very close in value, and therefore we base the temperature dependence of k diss on the following empirical expression for k 3 (25-48°C) given in Plummer et al. (1978) ( ) The effect of fluid pressure (neglecting the effect of 2 CO p -see next section) on the rate constant k diss is not considered, given the lack of pertinent data in the literature.  (Plummer et al., 1978;Busenberg and Plummer,1986;Arakaki and Mucci, 1995 In addition, these authors have also shown that the ionic strength has no effect on the dissolution rate constant. In our model, we use the results of Pokrovsky et al. (2004) for defining the change in k diss as a function of 2 CO p .

Coupling between dissolution-diffusion-precipitation processes and matrix deformation
Using our model we compute how the petrophysical properties of the rock matrix (permeability, porosity) and the volumetric strain (compaction) evolve with time due to PSC. The rock aggregate is a monomineralic limestone (i.e., pure calcite), with either a homogeneous or variable grain size (i.e. random spatial variation based on a uniform grain size distribution). The model treats the coupling between chemical reactions of calcite dissolution and precipitation, diffusion, and PSC deformation in two steps: textural equations at the grain scale allow for grain deformation due to chemical reactions occurring at contact and pore surfaces (including diffusion in the contact fluid); -a mass conservation relationship at a global scale which takes into account the transfer of dissolved material by diffusion in the interconnected pore spaces.
Therefore, microscopic grain scale processes and macroscopic mass transport within the model volume are fully coupled (for a complete review, see Gundersen et al., 2002).
The basis of our model for treating the effect of elevated normal stress on individual grains comprising a rock matrix is the thermodynamic relationship which relates normal stress to chemical potential (Gibbs, 1878;Kamb, 1961;Paterson, 1973;Lehner, 1995;Renard et al., 1999). The higher chemical potential of calcite surfaces in the contact zone results in a higher solubility with respect to calcite surfaces of the pore space. Ultimately, this difference in chemical potentials drives two processes: stressenhanced dissolution of the contact zone surfaces and diffusion of dissolved material from the contact zone out to the pore fluid. The diffusion of dissolved material occurs within a trapped water film separating the grains (Rutter, 1976;de Meer et al., 2002;Dysthe et al., 2002b).
Finally, grain-scale dissolution-precipitation processes in pore spaces are coupled to bulk diffusion of solute within the entire porous medium. The conservation relationships are derived by a global mass balance of the solute phase in the pore volume and a local mass balance at each grain contact (Gundersen et al., 2002). These relationships are then coupled to equations which express the deformation of the grains and the evolution of the rock texture.

Deformation of individual grains and rock matrix
In our model the rock matrix is a 20x20 cm domain, located at a specified depth.
The geological conditions (stress at the boundaries, pore pressure, temperature) are assumed to stay constant and are chosen to represent relevant conditions for CO 2 sequestration (depths between 1 and 3 km, see Table 3).
The rock matrix is modeled as a network of solid calcite grains with well-defined grain-grain contacts and interconnected pore spaces between grains (Fig. 1). The aggregate grain framework is represented by a dense cubic packing of truncated spheres that represent the individual grains. In Fig. 1 it is also important to distinguish the two different types of grain-fluid interfaces that are treated in the model: -grain surfaces at intergranular (i.e. grain-grain) contacts that are separated by a trapped fluid film, -'free' grain surfaces that are in contact with pore fluids.
The intergranular grain boundary is considered as a flat interface (Hickman and Evans, 1995) with a mean roughness of several nanometers averaged over the contact.
This interfacial contact zone contains a trapped fluid film that is postulated to have the unusual characteristic of being able to support a shear stress transmitted by the normal stress imposed on the grains. In addition to nanoscale topography, there is also ample evidence that grain-grain contact zones have a structure of channels and islands at a micrometer scale (Hickman and Evans, 1992;Schutjens and Spiers, 1999;Renard et al., 1999;Dysthe et al., 2002a). The exact relation between these nano-and microscale features is not yet well understood.
The initial radius of the spherical grains, L f , can vary between the different elements. The other lengths, describing the truncations of the spheres, L x and L z , are given . The choice of this value allows initial porosities close to 30%. Each simulation lasts until the porosity in the whole simulation domain is less than a threshold value equal to 5%. Below this value, transport by diffusion between the pores ceases and pressure solution only continues locally inside this element, as a closed system. The decrease in porosity from an initial 30% to a final value of 5% represents a finite volumetric strain of approximately 20%. The time associated with this porosity reduction defines the average rate of deformation for PSC and forms the basis for comparing the various cases examined (variable depth, 2 CO p ) The entire matrix of grains with cubic packing is subjected to a constant normal stress component (σ n ). We then define a contact normal stress component (σ c ) as the mean stress at each grain contact surface. This stress depends on the relationship between the surface area and the diameter of the sphere. This relationship is independent of the initial size of the individual grains (Dewers and Ortoleva, 1990), and the stress at the contact is proportional to some (positive) power of the porosity.
The evolution of the rock texture, which leads to a compaction of the aggregate, is computed as a result of the coupling between dissolution, diffusion, precipitation and mass transport in the fluid. This textural model allows us to take into account the sequential coupling between chemical processes (dissolution at grain contacts, diffusion, and precipitation on free pore surfaces) and the mechanical evolution (deformation of the grains) of the rock matrix. The resulting model is then a set of highly coupled non-linear equations that can only be solved using numerical methods. All of the parameters used in the following equations and their units are given in Table 1.

Matrix-fluid chemical equilibrium at t = 0
The initial state of the model is based on two separate equilibrium chemical states, such that at t = 0, the rock matrix is in equilibrium with, respectively, the contact and Thus, an initial chemical imbalance exists between the contact zones and the pore spaces.
The evolution of the chemistry of the contact fluid and the pore fluid is treated in a different manner once deformation starts, that is for t > 0. The initial chemical imbalance between the contact zone fluid and the pore space fluid drives diffusion that initiates the PSC deformation of the rock matrix. This diffusion process, which is constrained by the code to depend only on the negative [

Conservation equations in the pore fluid
Defining c p, Ca to represent the concentration of where the lengths L x , L y , L z are geometric grain parameters defined in Fig. 1, D  Ca on the pore walls (i.e. pore free faces) and the solute moles derived from diffusion from the grain contacts to the pore fluid, respectively (see Table 1 for symbols). The first term on the right side in Eq. (11) describes the diffusive transport of solutes in the pore space by a standard Fickian diffusion relation.
The second term on the right hand side in Eq. (11) represents mass loss due to precipitation on the grain surfaces and is given as: where k p is the rate constant for calcite precipitation, A p is the grain surface area in contact with the pore fluid (reactive area), is the ion activity product in the pore fluid, and K sp,p is the calcite solubility product in the pore fluid. Notice that as the geometry of the grains changes during compaction, A p varies as well. At conditions close to equilibrium, the rate expression used in the model to calculate the overall reactivity of calcite does not differentiate between dissolution and precipitation (i.e. k diss = k prec ). This appears to be a reasonable approximation based on available experimental evidence at conditions very close to equilibrium (e.g. Fig. 8, Busenberg and Plummer, 1986).

Conservation equations at the grain contacts and grain shape evolution
The last term in Eq. (11) describes mass production due to diffusional mass transport out of the intergranular contact zone and reads: where c c,i,Ca is the concentration of + 2 Ca in the grain contact perpendicular to axis i (i=x, y, z). The driving force for this transport is a concentration gradient related to the difference in stress concentration at intergranular and free grain contacts. In the model, the length of the intergranular diffusion path relates directly to the rate of compaction (i.e. smaller grain matrices compact faster than larger grains). The transfer of dissolved + 2 Ca within the interfacial fluid from grain-grain contacts to the pore fluid occurs via a diffusion process. Thus, the chemistry of the interfacial fluid evolves not only by the addition of material via dissolution, but also by the removal of dissolved material by diffusion.
In general, the rate of ion diffusion (D c ) in thin, confined interfacial fluid films is not directly measurable, but recent studies (e.g. Dysthe et al., 2002b;Alcantar et al., 2003) conclude that the rate is no more than two orders of magnitude slower than the corresponding bulk fluid value (at a given T and P). The adopted value for the diffusion coefficient D c of + 2 Ca along the grain contact used in the model is 0.01·D p (D p is the diffusion coefficient in the pore fluid) at the T and P of the contact fluid. The coefficient of diffusion within the water film, D c , follows an Arrhenius law with an activation energy of 15 kJ/mol (Dysthe et al., 2002b). Considering the uncertainty in D c , we can estimate that the true value probably differs by at least one order of magnitude from the value we have chosen. In Eq. (13) the thickness of the water film at the grain contacts, Δ, is divided by two since it is shared between two contact surfaces. We assume that this water film has a constant thickness of 3 nm (Dysthe et al., 2002b).
The global mass balance equation (Eq. 11) is then coupled to Eq. (14) which represents the local mass balance for each contact on a truncated spherical grain. The concentration of dissolved material within the intergranular fluid phase is given as the difference between the mass produced by dissolution and the mass lost by diffusion: In the expression above, R c,i is the radius of the contact surface in the i-direction, which is a function of the truncations of the spherical grains. The first term in Eq. (14) represents the local production of dissolved material by dissolution at the grain contacts. The model defines the chemical reaction (i.e. dissolution) of grain-grain contacts in terms of the rate law already presented in Eq. 8, where R is the overall rate (mol m -2 s -1 ), k diss is the rate constant, Q c,i is the ion activity product in the contact zone, and K sp,c,i is the calcite solubility product in the intergranular contact zone perpendicular to the i axis (i=x, y, z): Using this law, one can write the flux corresponding to the dissolution process as: is not limited to an increase in the rate of dissolution (via k diss ), however, since this also results in a concomitant increase in calcite solubility (i.e. K sp,c and K sp,p ), which will also affect the rate of PSC.
The grain shape evolution is given by the equations where V is the molar volume of calcite (3.69·10 -5 m 3 mol -1 ) and the calcite precipitation rate k prec is equal to the dissolution rate constant k diss . With these equations the code continually updates the texture of the grains and thus couples chemical reactions to grain deformation.
It is important to note that even though the two kinetic processes, dissolution (Eq. 17) and precipitation (Eq. 18), respectively control the grain shape evolution in the contact zones and the pore spaces, these processes are coupled to diffusion in the contact zones, and therefore are not independent. This coupling results in a modification of the chemical affinity terms in both rate equations above (bracketed terms), thereby changing the overall rate of grain shape evolution. As an example, taking the rate of dissolution in The pore space, however, is a special case since dissolved material enters by diffusion from adjacent contact zones, but it can also diffuse to other pore spaces further away.

Stress at grain contacts
The simulation domain is a square rock matrix located at depths of 1, 2 or 3 km and is exposed to constant lithostatic stress and temperature. The normal stress component on an i-contact, σ c, i , is related to the lithostatic stress σ n on the boundaries Γ i by the three relationships which take into account the stress concentration effects due to the porosity of the rock (Dewers and Ortoleva, 1990). These three equations are obtained by circular permutation of the indices in the following relationship: where i, j, k represent the three space directions x, y, and z.
When the lithostatic stress is kept constant, the contact surface area A i increases as the grains become more and more truncated (Fig. 1). As a consequence the stress normal to this contact decreases with time, and therefore, the driving force for pressure solution also decreases with time. All of these phenomena are due to the continual removal of material at grain-grain contacts by dissolution and the outward diffusion of this dissolved material into the pore fluid.
The final equation updates the porosity, which is defined as the difference between the volume of a cubic box, with lengths L x , L y , L z , and the volume of the truncated grains (Dewers and Ortoleva, 1990). This relationship arises from the textural model of truncated spheres ordered in a centered cubic network (Fig. 1).

Numerical methods and boundary conditions
The equations presented in sections 2.4-2.6 are highly coupled. The deformation of the grains is coupled to transport in the pore fluid through the diffusion step at the contacts and the Fickian diffusion in the pores. Therefore we choose two different numerical modelling techniques. We couple a discrete treatment of the grains (Eqs. 17-

18) with a continuous determination of the dissolved
Ca in the pore fluid (Eq. 11). The space domain is a 2D rock sample in the x-z plane and is modeled as a cubic packing of truncated spheres. It has a thickness of one layer of grains in the y direction (Fig. 1) The geometrical evolution of the grains is modeled as a discrete process where the grains in each grid element are treated separately (see . In order to visualize the deformation of the rock matrix as the grains deform during PSC, an adaptive grid is developed. The nodes change position after each time step as the grains deform. In addition, the contact surface area of each element is updated according to the lengths and radii of the truncated spheres, and is subsequently used to calculate the porosity loss in the xz-plane. The total change in volume of each element is then obtained by multiplying the volume variation of each grain by the number of grains in the element. The total deformation of the rock sample is determined by a summation over all the elements. The mass conservation equation (Eqn. 11) is derived from a continuum model, assuming a continuous solute phase. This equation is solved using a Galerkin finite element method, with linear basis functions on the spatial domain and a Crank-Nicolson scheme in the time domain. All the numerical schemes were programmed in C++ and make use of the numerical library Diffpack (Langtangen, 1999;Gundersen et al., 2002).
The program was extensively tested for numerical stability. In all of the simulations, we verified that the total mass of solid was conserved as it should be.
The time step is adjusted such that the maximum deformation is less than 0.1% in all the elements between each time increment. Therefore, the time step automatically adapts in the regions with the most rapid deformation.

Results
In our model, three different geometries of rock matrix are studied: -a single homogeneous layer with an initially constant grain size; -a sedimentary layered rock with an initial local variation of grain size; -a gouge filling a fracture.
In the results presented here, we aim to compare the evolution of the rock with and without CO 2 injection. In order to compare the various simulations, we have chosen a reference case that provides a normalized time scale for the measured PSC rates: the reference (see Figs. 2a and 6) is based on the time needed to achieve a porosity of 5% for a rock at 1 km depth and 40 °C, with an initial grain size of 2 mm (L f = 1 mm), and a pore fluid that is CO 2 -free (i.e. 2 CO p = 10 -4.5 MPa). In all the simulations, the times scales of deformation are normalized to this reference case (total time for reference deformation equals 740,000 years-see Table 3). In this way we measure the degree to which the sequestration of CO 2 perturbs the system and enhances deformation and compaction by PSC. The absolute time scales for PSC deformation for all of the simulations are given in Table 3. Because of uncertainties with respect to several parameters (k diss , k prec , D c ), the absolute time scales related to PSC deformation are difficult to quantify and should be used with caution (see Discussion).
In all cases, both CO 2 -free and at elevated 2 CO p , the PSC process is controlled by the diffusion process in the contact zones. This indicates that the intergranular diffusion process is slower than either dissolution in the contact zones or precipitation in the pores.
Diffusion in the contact zones, controlled by the + 2 Ca concentration gradient between the intergranular fluid and the pore fluid, therefore, plays the dominant role in defining the rate of diffusion, the rate of grain indentation, and most importantly, the overall rate of PSC deformation of the rock matrix. This result is predictable, given the rapid kinetics of calcite dissolution and precipitation. Figure 2c, which shows that the + 2 Ca concentration in the contacts is always greater than that in the pore fluid, graphically demonstrates that the rate of diffusion is the limiting process of PSC for these simulation conditions.

CO 2 -free simulations
The rock matrix, located at 1 km depth, consists of grains with a homogeneous size of 2 mm, which corresponds to a coarse oolitic limestone. We consider here a CO 2free pore fluid (i.e. atmospheric 2 CO p = 10 -4.5 MPa). In our model, a typical compaction process without CO 2 can last between 10 4 to 10 6 years, depending on the choice of the parameters D c , k prec, k diss . As already stated, this configuration represents a reference case for comparison with other simulations at elevated 2 CO p .
As PSC proceeds, the porosity decreases progressively with time both by grain indentation and by precipitation in the pores (Fig. 2a). As shown in Fig. 2b, L f increases whereas L z decreases. With time the grains become truncated and the stress on the contacts decreases. As a consequence, the gradient in [ Ca ] between the grain contacts and the pore fluid also decreases (Fig. 2c), leading to a decrease in the rate of overall PSC deformation For a heterogeneous medium, such as a layered rock, the grain size may vary spatially (Fig. 3). In this case, the PSC rate is faster in the layers with small grains than in layers with coarse grains. Because of this initial grain size difference, the different domains compact at different rates. Layers with smaller grains reach the limit of 5% porosity faster than coarser grain layers. This result is an example of the well-known inverse dependence of pressure solution creep rate on grain size (Rutter, 1976;Gratz 1991). Similar behavior can be observed for a fracture filled with a gouge containing smaller grains (Fig. 4). Both examples show that compaction is homogeneous in the domains with a constant grain size. We do not observe mass transfer between layers as we did for the simulation of PSC in sandstones (Gundersen et al., 2002). This is related to the faster compaction rate in limestones, and the smaller time scales involved for long distance diffusion in the pores.

Effect of elevated partial pressures of CO 2
When CO 2 is present at hydrostatic pressure (i.e. 2 CO p = P p ), PSC is significantly faster compared to the CO 2 -free case. This situation is clearly visualized in Figure 5 which shows the compaction process for two fractured rocks at 2 km depth, one with 2 CO p =10 -4.5 MPa (i.e. CO 2 -free), the other where 2 CO p = P p (see Table 3). The pattern of the porosity reduction is similar in both cases, however the characteristic time scales are very different; at 2 CO p = P p , the process is roughly 65 times faster.
We also compare the compaction process at different depths with and without injection of CO 2 . With no injection of CO 2 the rock matrix compacts naturally as PSC progresses. As depth increases, the temperature and the stress increase as well. On the one hand, the effect of stress enhances PSC, mainly by increasing the solubility and the rate of dissolution. On the other hand, increasing the temperature decreases the solubility of calcite (i.e. calcite has retrograde solubility), and this therefore slows down the PSC rate. These two competing factors explain why the PSC rates increase by less than a factor of 3 between 1 and 3 km burial depth (Figure 6a). In a similar manner, for elevated 2 CO p , the relative increases in the PSC rates with increasing depth are modest (Table 3), again reflecting the competing effects of stress, solubility and kinetics.
When CO 2 is present at elevated concentrations, where 2 CO p equals the pore fluid pressure, the rates of PSC are greater by factors of ~50-75 (Fig. 6b), as compared to the CO 2 -free simulations (Fig. 6a), the exact amount depending on depth (Table 3). This translates to a reduction in porosity significantly more rapid than compaction at the same depth in the absence of CO 2 . This can be explained by two factors. First, at elevated levels of CO 2 , the pH decreases dramatically (  (Pokrovsky et al., 2004) and precipitation.
However, an increase in the rate of calcite dissolution in the intergranular zone should not affect the rate of deformation since it appears that the diffusion process is the rate-limiting step in the overall PSC process. On the other hand, an increase in the rate of precipitation in the pores has a positive effect on the rate of diffusion from the contact zone to the pores. What, if any, effect increased levels of CO 2 has on the rates of diffusion is unknown.

Discussion and conclusions
The injection of CO 2 into a geological reservoir and the resulting increase in the 2 CO p cause acidification of the pore fluids. As the pH of the pore fluids decreases, they become more reactive with respect to the rock matrix; this effect is especially pronounced for carbonates since both the solubility and the reaction kinetics increase dramatically with decreasing pH. This is in accord with the general finding that, to a first approximation, PSC rates depend linearly on the solubility of minerals (Rutter, 1976).
Our model shows that for a given depth and temperature, elevated concentrations of dissolved CO 2 in the pore fluids lead to rates of PSC deformation that are up to 75 times faster than reference simulations based on CO 2 -free pore fluids. Thus, our results convincingly demonstrate the direct relation between elevated 2 CO p levels, calcite solubility, pore fluid pH, and increased rates of PSC deformation.
Of particular importance to predicting the behavior of future CO 2 geological repositories is the accuracy of theoretical PSC deformation rates. Taking just one example from our model, the absolute time scale associated with porosity reduction from 30 to 5% for a carbonate rock matrix with CO 2 -free pore fluids at 1 km depth (40 °C, 2mm grain size) is on the order of 700,000 years (Table 3). This time frame is most probably unrealistically rapid, given that in nature PSC occurs over much longer time scales. As an example, it is observed that Mesozoic chalks and limestones retain some porosity after several tens of million years, which represents a time period several orders of magnitude greater. There are undoubtedly a myriad number of reasons for this, for which a detailed analysis is beyond the scope of this article. Nonetheless, we briefly discuss how some important parameters influence theoretical PSC deformation rates derived from our model, and how they may contribute to discrepancies with natural PSC rates.
The present study reports 'preliminary results' that are derived from a model that incorporates many simplifications that are both chemical and physical in nature. H into a consistent model. Even though a fully coupled model should produce more realistic PSC rates, we estimate that the deformation patterns would not differ too much from those generated in the present study.
The thermodynamic, kinetic, and diffusion parameters that are incorporated in the code play perhaps the most important role in determining the accuracy of theoretical PSC rates of deformation. The rates that we have chosen for the kinetic rate parameters k diss and k prec , the value of n (see Eq. 8), and the overall form of the kinetic rate equation(s), determine the rates of dissolution and precipitation in the contact zone and pore space, respectively, and therefore potentially influence the overall rate of PSC deformation. The kinetic parameters we use are not unique and therefore future modeling will include sensitivity analyses based on a variation of these kinetic parameters. In addition, the accuracy of theoretical PSC rates depends on the availability and quality of experimental data, both kinetic and thermodynamic, for conditions relevant to CO 2 sequestration.
Further refinements in PSC models will depend on additional experimental work to unravel more precisely the effect of elevated 2 CO p on k diss and k prec of calcite (and other minerals, as well) at conditions close to equilibrium at elevated T and P.
If we consider CO 2 sequestration and PSC deformation within the context of conditions in a geological repository, the chemical complexity of natural pore fluids may dictate the need for additional kinetic parameters to be used in the kinetic rate laws. As an example, the presence of certains ions (e.g. PO ; Svensson and Dreybrodt, 1992;Zuddas and Mucci, 1998;Zhang et al., 2002) or humic acids (e.g. Zuddas et al., 2003) have been shown to inhibit calcite dissolution/precipitation rates.
The inhibition of calcite kinetic rates has an important implication for carbonate PSC deformation since the limiting process may switch from diffusion limited to reaction limited. This points out that more experimental work is needed to investigate the effect of fluid chemistry on rates of PSC deformation. The mineralogy of the solid phase is also important. At present, our code only considers a monomineralic (pure calcite) limestone.
The effect of impurities, such as clays, potentially has a large effect on PSC deformation.
For example, it has been experimentally shown that the presence of clays can dramatically increase rates of PSC (Hickman and Evans, 1995;Renard et al., 2001).
Even though the thermodynamics of the calcite-H 2 O-CO 2 system at STP are well known and quite accurate, future versions of our code need to further develop the temperature and pressure dependencies of most of the equilibria given in Eqs. 3a-h.
Accounting for non-ideal behavior at high 2 CO p is also necessary with respect to Henry's law constant in Eq. 3b. Moreover, future versions of the code will need to consider equilibria relations in terms of activities and not concentrations; this will be especially important for solutions with elevated ionic strengths. Two other parameters of considerable importance, especially for the case of diffusion-limited PSC deformation, are the rate of diffusion and the thickness (and continuity) of the fluid film in intergranular contacts. Measurement of diffusion coefficients in thin films is experimentally challenging, and therefore the diffusion coefficient and film thickness will remain uncertain quantities that probably contribute significantly to the discrepancy between natural and theoretical PSC deformation rates.
At this stage in the development of our model, many uncertainties exist with respect to the thermodynamic, kinetic, and physical parameters that have been incorporated in the code. Future versions of the code will hopefully be able to address many of these issues. We believe that the main value of the results generated in the present study lies not in the absolute rates, but rather in the relative rates of PSC deformation. Intercomparison of these relative PSC rates will help in better understanding and evaluating the long-term effect of increased rates of PSC on the porosity and viscosity of rock matrices in contact with high 2 CO p fluids. Another important benefit of PSC models such as ours is that they may help in focusing laboratory-based PSC experiments. The relation between solubility, pH and PSC rates has important ramifications with respect to future experimental laboratory protocols for investigating the sequestration of CO 2 . As an example, the injection of subcritical CO 2 -H 2 O mixtures into a rock sample could be simplified by the simple injection of an acidic aqueous fluid, where the pH is adjusted to be equivalent to the pH due to CO 2 acidification alone.
Nonetheless, such an approach has its limitations and would not be applicable under all circumstances, especially in the case of injection of supercritical CO 2 fluids.