Reactive Transport Modeling: Applications in Subsurface Energy and Environmental Problems

Teaches the application of Reactive Transport Modeling (RTM) for subsurface systems in order to expedite the understanding of the behavior of complex geological systems 

This book lays out the basic principles and approaches of Reactive Transport Modeling (RTM) for surface and subsurface environments, presenting specific workflows and applications. The techniques discussed are being increasingly commonly used in a wide range of research fields, and the information provided covers fundamental theory, practical issues in running reactive transport models, and how to apply techniques in specific areas. The need for RTM in engineered facilities, such as nuclear waste repositories or CO2 storage sites, is ever increasing, because the prediction of the future evolution of these systems has become a legal obligation. With increasing recognition of the power of these approaches, and their widening adoption, comes responsibility to ensure appropriate application of available tools. This book aims to provide the requisite understanding of key aspects of RTM, and in doing so help identify and thus avoid potential pitfalls.

Reactive Transport Modeling covers: the application of RTM for CO2 sequestration and geothermal energy development; reservoir quality prediction; modeling diagenesis; modeling geochemical processes in oil & gas production; modeling gas hydrate production; reactive transport in fractured and porous media; reactive transport studies for nuclear waste disposal; reactive flow modeling in hydrothermal systems; and modeling biogeochemical processes. Key features include:

• A comprehensive reference for scientists and practitioners entering the area of reactive transport modeling (RTM)
• Presented by internationally known experts in the field
• Covers fundamental theory, practical issues in running reactive transport models, and hands-on examples for applying techniques in specific areas
• Teaches readers to appreciate the power of RTM and to stimulate usage and application

Reactive Transport Modeling is written for graduate students and researchers in academia, government laboratories, and industry who are interested in applying reactive transport modeling to the topic of their research. The book will also appeal to geochemists, hydrogeologists, geophysicists, earth scientists, environmental engineers, and environmental chemists. 

• Language: English
• Publisher: Wiley
• Release date: Mar 14, 2018
• ISBN: 9781119060024

Book preview

1

Application of Reactive Transport Modeling to CO2 Geological Sequestration and Chemical Stimulation of an Enhanced Geothermal Reservoir

Tianfu Xu, Hailong Tian Jin Na

Key Laboratory of Groundwater Resources and Environment, Ministry of Education, Jilin University, Changchun, 130021, China

This chapter is devoted to giving a deep insight into the applications of reactive transport modeling to two problems where geochemistry plays an important role. The two problems involved here are CO2 geological sequestration and the development of geothermal energy. These applications elucidate the use of reactive transport modeling in understanding the reactions among water, solute and minerals occurring in the areas of CO2 geological sequestration and geothermal energy development.

1.1 Introduction

Coupled modeling of subsurface multiphase fluid and heat flow, solute transport, and chemical reactions can be applied to many geological systems and environmental problems, including geothermal systems, diagenetic and weathering processes, subsurface waste disposal, acid mine drainage remediation, contaminant transport, and groundwater quality.

The processes involved in coupling include geochemical, thermal and hydrological interactions. These processes are strongly coupled through linear and non‐linear relationships, and make up a significant part of the subject area of ‘reactive transport’, in which we focus on the development of models that quantitatively represent the evolution of subsurface geological environments.

Reactive transport models that consider conductive and convective heat transport, dissolution and precipitation kinetics, aqueous speciation and coupled porosity evolution within saturated and unsaturated geological materials are very demanding in terms of computing time and memory, and results might either be difficult to obtain or limited in terms of precision. Very recently, such modeling work has achieved a high degree of sophistication as access to high‐performance computer increases.

Recent numerical modeling for subsurface environment applications can be split into three categories, including hydrodynamic modeling, geochemical modeling (batch reaction) and reactive transport modeling. Reactive transport modeling is the most realistic modeling technique to quantify the complex of coupling processes over a long period of time.

For each category there are special modeling packages for performing specific simulations. For example, ECLIPSE, TOUGH2, SIMED, and so on, were developed for hydrodynamic modeling and associated issues. For batch geochemical modeling, usages of PHREEQC and TOUGHREACT are reported in previous works (Gaus et al., 2004; Xu et al., 2004a). Batch geochemical modeling is ideal when modeling laboratory batch experiments. However, it becomes limited when the hydrodynamics plays a dominant role in the whole process, and then TOUGHREACT, SIMUSCOPP, STOMP, HYTEC, UNFT, SUPCRT92 and GEMBOCHS were developed for reactive transport modeling.

Among these codes, TOUGHREACT was developed by introducing geochemistry into the existing TOUGH2, a framework of a non‐isothermal, multicomponent fluid and heat in a porous and fractured media simulator. This program has been widely used for many geological systems and environmental problems. For the worldwide recognition of TOUGHREACT, this chapter will briefly introduce the fundamental theories of reactive transport modeling, and illustrate its basic applications to two problems: CO2 geological sequestration and geothermal energy development.

1.2 Fundamental Theories

1.2.1 Governing Equations for Flow and Transport

The primary governing equations for multiphase fluid and heat flow and chemical transport have the same structure, derived from the principle of mass (or energy) conservation. These equations are presented in Table 1.1.

Table 1.1 Primary governing equations for fluid and heat flow, and chemical transport. Symbol meanings are given in Table 1.2.

Table 1.2 Symbols used in Table 1.1.

Aqueous (dissolved) species are subject to transport in the liquid phase as well as to local chemical interactions with the solid and gas phases. The transport equations are written in terms of total dissolved concentrations of chemical components, which are concentrations of the basis species plus their associated aqueous secondary species (Yeh and Tripathi, 1991; Steefel and Lasaga, 1994; Xu and Pruess, 2001). If kinetically controlled reactions occur between aqueous species, then additional ordinary differential equations need to be solved to link the total concentrations of the primary species with the evolving concentrations of the secondary species (Steefel and MacQuarrie, 1996). Advection and diffusion processes are considered for both the aqueous and gaseous species. Aqueous species diffusion coefficients are assumed to be the same. Gaseous species, having a neutral valence, can have differing diffusion coefficients calculated as a function of T, P, molecular weight, and molecular diameter.

The local chemical interactions in the transport equations are represented by reaction source/sink terms.

1.2.2 Equations for Chemical Reactions

To represent a geochemical system, it is convenient to select a subset of NC aqueous species as basis species (or component or primary species). All other species are called secondary species that include aqueous complexes, precipitated (mineral) and gaseous species (Reed, 1982; Yeh and Tripathi, 1991; Steefel and Lasaga, 1994). The number of secondary species must be equal to the number of independent reactions. Any of these secondary species can be represented as a linear combination of the set of basis species, such as

(1.1)

where Si represents chemical species, j is the basis species index, i is the secondary species index, NR is the number of reactions (or secondary species), and vij is the stoichiometric coefficient of the jth basis species in the ith reaction.

Dissolution/precipitation for equilibrium minerals

The equations for chemical equilibrium are similar to those by Parkhurst et al. (1980), Reed (1982), Yeh and Tripathi (1991), Wolery (1992), and Steefel and Lasaga (1994), and are presented in Table 1.3.

Aqueous species activity coefficients are calculated from the extended Debye‐Hückel equation (Helgeson et al., 1981). The calculations of gas fugacity coefficients are given in Spycher and Pruess (2005).

The mineral saturation ratio can be expressed as:

(1.2)

where m is the equilibrium mineral index and Km is the corresponding equilibrium constant. At equilibrium, we have:

(1.3)

where SIm is the mineral saturation index.

Kinetic mineral dissolution/precipitation

Kinetic rates can also be functions of non‐basis species. Usually the species appearing in rate laws happen to be basis species. In this model, we use a rate expression given by Lasaga et al. (1994):

(1.4)

where positive values of rn indicate dissolution, and negative values precipitation, kn is the rate constant (moles per unit mineral surface area and unit time) which is temperature dependent, An is the specific reactive surface area per kg H2O, and n is the kinetic mineral saturation ratio. The parameters must be determined from experiments; usually, but not always, they are taken equal to one.

The temperature dependence of the reaction rate constant can be expressed reasonably well via an Arrhenius equation (Lasaga, 1984; Steefel and Lasaga, 1994). Because many rate constants are reported at 25 °C, it is convenient to approximate rate constant dependency as a function of temperature thus:

(1.5)

where Ea is the activation energy, k25 is the rate constant at 25 °C, R is the gas constant, and T is the absolute temperature.

For many minerals, the kinetic rate constant k can be summed from three mechanisms (Palandri and Kharaka, 2004):

(1.6)

where superscripts or subscripts nu, H, and OH indicate neutral, acid, and base mechanisms, respectively, Ea is the activation energy, k25 is the rate constant at 25 °C, R is the gas constant, T is absolute temperature, a is the activity of the species, and n is an exponent (constant). The rate constant k can be also dependent on other species such as Al³+ and Fe³+. Two or more species may be involved in one mechanism by:

(1.7)

where superscript or subscript i is the additional mechanism index, and j is the species index involved in one mechanism that can be primary or secondary species.

Table 1.3 List of equations governing chemical equilibrium (illustrated by specific examples; in fact the model is valid for general geochemistry). Symbol meanings are given in Table 1.4.

Table 1.4 Symbols used in Table 1.3.

1.2.3 Solution Method for Transport Equations

Most chemical species are only subject to transport in the liquid phase. A few species can be transported in both liquid and gas phases, such as O2 and CO2. We first derive the numerical formulation of reactive transport in the liquid phase. This will then be extended to transport in the gas phases for some gaseous species.

Transport in the liquid phase

In the sequential iteration approach (SIA), the mass transport equations and chemical reaction equations are considered as two relatively independent subsystems. They are solved separately in a sequential manner following an iterative procedure. If reactions taking place in the liquid phase are assumed to be at local equilibrium, mass transport equations can be written in terms of total dissolved component concentrations. By lumping all mass accumulation terms due to mass transfer between aqueous and solid phases including precipitated (kinetics and equilibrium), exchanged and sorbed species, we can write equations for multicomponent chemical transport in the liquid phase as:

(1.8)

where n labels the grid block, m labels the adjacent grid blocks connected to n, j labels the chemical component, NC is the total number of chemical components, l labels the liquid phase, k labels the number of the time step, s labels the number of the transport‐chemistry iteration, unm is the liquid volumetric flux or the Darcy velocity (m/s), Dnm is the effective diffusion coefficient (including effects of porosity, phase saturation, tortuosity and weighting factors between the two grid blocks), dnm is the nodal distance, and are the overall chemical reaction source/sink terms.

Transport in the gas phase

Gaseous species concentrations can be related to partial pressures by

(1.9)

where Cg are gaseous species concentrations (in mol/l), Pg is the gaseous species partial pressure (in bar), R is the gas constant (8.314 J mol−1 K−1) and T is the absolute temperature. By following the same principle as used for transport in the liquid phase and by considering Equation (1.9), the numerical formulation of gaseous transport in the gas phases can be expressed as:

(1.10)

where Ng is the number of gaseous species.

1.2.4 Solution Method for Mixed Equilibrium‐Kinetics Chemical System

Aqueous complexation and gas dissolution/exsolution proceed according to local equilibrium, while mineral dissolution/precipitation is subject to equilibrium and/or kinetic conditions. Gas dissolution/exsolution is included in the model and treated in a similar way as equilibrium mineral dissolution/precipitation, but with fugacity correction. The formulation is based on mass balances in terms of basis species as used by Parkhurst et al. (1980) and Reed (1982) for the equilibrium chemical system. The kinetic rate expressions for mineral dissolution/precipitation are included in the equations along with the mass balances of basis species. At time zero (initial), the total concentrations of basis species j in the system are assumed to be known, and are given by

(1.11)

where superscript 0 represents time zero; c denotes concentration; subscripts j, k, m, n, z and s are the indices of the primary species, aqueous complexes, minerals at equilibrium, minerals under kinetic constraints, and exchanged and surface complexes, respectively; Nc, Nx, Np, Nq, Nz and Ns are the number of corresponding species and minerals; and vkj, vmj, vnj, vzj and vsj are stoichiometric coefficients of the primary species in the aqueous complexes, equilibrium, kinetic minerals, and exchanged and surface complexes, respectively.

After a time step ∆t, the total concentration of primary species j (Tj) is given by

(1.12)

where rn is the kinetic rate of mineral dissolution (negative for precipitation, units used here are moles per kilogram of water per time), for which a general multi‐mechanism rate law can be used. ∆t and Tj are related through generation of j among primary species as follows:

(1.13)

where l is the aqueous kinetic reaction (including biodegradation) index, Na is the total number of kinetic reactions among primary species, and rl is the kinetic rate, which is in terms of one mole of product species.

By substituting Equations (1.11) and (1.12) into Equation (1.13), and denoting residuals as F cj (which are zero in the limit of convergence), we have:

(1.14)

According to mass‐action equations, the concentrations of aqueous and exchanged complexes ck and cz can be expressed as functions of concentrations of the primary species cj. Kinetic rates rn and rl are functions of cj. Surface complexes are expressed as the product of primary species and an additional potential term Ψ. Additional equilibrium equations depending on the surface complexation model have to be solved together with Equation (1.14). No explicit expressions relate equilibrium mineral concentrations cm to cj. Therefore, NP additional mass action equations (one per mineral) are needed. At equilibrium we have the residual functions

(1.15)

where SIm is the mth equilibrium mineral saturation index, which is a function of cj. NC equations in (1.14) and Np in (1.15) constitute a complete set of equations needed for NC + NP primary variables (c1, c2,…, cNc; c1, c2,…, cNp).

Equations (1.14) and (1.15) can be solved together with the robust Newton‐Raphson iterative method.

1.3 Application to CO2 Geological Storage (CGS)

1.3.1 Overview of Applications in CGS

Carbon capture and storage (CCS) has been validated by the IPCC (2005) as an important option to mitigate the increasing atmospheric concentrations of anthropogenic carbon dioxide (CO2). In order for the technology to be deployed on a large scale, its viability in terms of injectivity, the containment of the injected CO2 and the long‐term safety with respect to humans and the environment needs to be guaranteed. Carbon dioxide geological storage (hereafter abbreviated to CGS) has been conducted on a commercial scale (e.g. a CO2 injection rate more than 1.0 Mt/a) at several sites including the Weyburn‐Midale field, Canada (Petroleum Technology Research Center, 2011). Due to the corrosive character of CO2, certainly once it is dissolved, geochemical reactions play an important role and might affect the chemical and physical properties of the wells, the reservoir and its surroundings.

Geochemical interactions caused by the presence of CO2 in geological sequences where CO2 occurs naturally are particularly valuable, since this illustrates the long‐term impact of CO2 on natural rocks that cannot be reproduced during experiments or field tests. Chemical equilibrium seems not to be reached in some natural analogues, even when subjected to very long (geological) contact time, suggesting that chemical equilibrium might not be the natural state of injected CO2 during the whole duration of the storage period (thousands to hundreds of thousands of years).

Numerical modeling of a storage site requires estimates for both the short‐ and long‐term fate of the injected CO2. Such modeling work can be very demanding in terms of CPU time and memory, and results might either be difficult to obtain or limited in terms of precision. Recent numerical modeling for CO2 storage applications can be split into three categories: (i) hydrodynamic modeling simulating structural, residual gas and dissolution trapping processes; (ii) geochemical modeling simulating batch geochemical reactivity (closed system without any fluid flow); and (iii) reactive transport modeling combining the two previous types of simulations, and therefore providing a complete calculation over time of the amount of CO2 trapped through structural, dissolution or mineral trapping. Reactive transport modeling is the most realistic modeling technique to quantify the long‐term fate of CO2 as well as other aspects during geological storage, but also the most challenging to perform.

Geochemical and solute transport modeling has many application domains when assessing the geochemical impact of CO2 storage, and each has its own spatial scale and timeframe of interest and therefore requires an adapted modeling approach.

One application of reactive transport modeling is to assess the ultimate fate of the injected CO2 and its impacts on physical properties. There are four distinct processes during sequestering the injected CO2, including structural trapping, residual CO2 trapping, dissolution trapping and mineral trapping. Mineral trapping kinetics is controlled by both dissolution and precipitation kinetics. While the dissolution kinetics for carbonates and sulfate reactions are generally fairly rapid, the kinetics of alumino‐silicate mineral reactions are much slower (total dissolution of these minerals can take up to thousands of years at ambient temperatures), making mineral trapping a slow process since mainly the latter are involved in CO2 trapping processes.

Another application domain of interest for researchers is to predict the risk potentially occurring in CGS through reactive transport modeling, such as caprock integrity and damage, CO2 leakage, or groundwater quality reduction. However, the caprock as the first natural barrier of the reservoir particularly concerns scientists. Caprock is generally defined as a low to very low permeability formation, and composed mainly of clay‐rich shale or mudstone that consists of aluminates or aluminosilicates. Because aluminosilicate mineral alteration is very slow under ambient deep‐formation conditions, and is not amenable to experimental studies, the numerical modeling of hydrogeochemical processes is necessary to investigate long‐term evolution of caprock sealing efficiency.

1.3.2 Long‐Term Fate of Injected CO2 in Deep Saline Aquifers

As the negative effects of global warming are widely discernible, the need for the deployment of CO2 capture and storage (CCS) incorporating CO2 geological storage (IPPC, 2005) has been increased to confront this worldwide environmental issue. Geological CO2 sequestration offers a most promising solution for reducing net emissions of greenhouse gases into the atmosphere. Injection of CO2 produced by different human activities into depleted oil and gas reservoirs has been considered, but aquifers containing non‐drinking water offer the largest storage volume. CO2 is injected into these aquifers in the supercritical state with a higher density, and therefore occupying less volume underground. The supercritical state can be obtained at pressures greater than 7.4 MPa and temperatures higher than 31.1 °C. These conditions correspond to a reservoir at a depth of 800 m beneath the Earth’s surface. There are three interdependent yet conceptually distinct processes storing CO2 in deep saline aquifers, which are distinguished as hydrodynamic trapping of the gas (or supercritical) phase, dissolution trapping in the liquid (groundwater) phase (dissolution of CO2), and mineral trapping of the solid phase (mineral alteration, leading to precipitation of secondary mineral phases).

In terms of migration, injected CO2 moves by volumetric displacement of formation waters, with which it is largely immiscible; by gravity segregation, which cause the immiscible plume to rise owing to its relatively low density; and by viscous fingering, which causes it to migrate preferentially into local high‐permeability zones owing to its relatively low viscosity compared with water. In terms of sequestration, some fraction of the rising plume will dissolve into formation waters (solubility trapping); some fraction will react with formation minerals to precipitate carbonates (mineral trapping); and the remaining fraction will reach and become isolated beneath the caprock (hydrodynamic trapping), migrate up‐dip along this interface, and accumulate in any local topographic highs (structural trapping). Numerical simulation of these interdependent migration and sequestration processes requires a computational capability that explicitly represents and couples multiphase flow and kinetically controlled geochemical processes within porous media characterized by physical and compositional heterogeneity.

The ultimate fate of CO2 injected into saline aquifers for environmental isolation is governed by the above‐mentioned three interdependent yet conceptually distinct processes. The first process is directly linked to hydrodynamic trapping, the second to solubility trapping and pH evolution, and the third to mineral trapping (and pH evolution). In this section, using the reactive transport modeling approach, we will quantify and compare the relative effectiveness of these trapping mechanisms during propagation of injected CO2 within a specified deep saline aquifer.

1.3.2.1 Brief Description of CO2 Storage Site in the Songliao Basin

In order to ensure the effectiveness and safety of any geological storage project, site selection is crucial to the long‐term containment of injected CO2 in the deep reservoir and a very important component in the life‐cycle assessment of CCS, and also plays a key role during the implementation process of CCS (Bachu, 2008; Li et al., 2012). Sites are chosen under different conditions, such as performance or risk. Integrated scales of site assessment are divided into four stages: country/state‐scale screening, basin‐scale assessment, site characterization, and site deployment (CO2CRC, 2008; CSLF, 2008). Screening criteria for CO2 sequestration have been proposed for reservoir properties and surface facilities, and different value ranges of indicators have been given (Kovscek, 2002; LBNL, 2004). Chadwick et al. (2008) defined the criteria as reservoir efficacy, reservoir properties and caprock.

The site used in this study, Songliao Basin, is located in the northeast of China, which spans Heilong Jiang, Jilin and Liaoning provinces and the Inner Mongolia Autonomous Region. The basin was formed first by rifting and then by subsidence and inversion. The Songliao Basin is one of the most important productive oil basins in China, as well as one of the most important heavy industrial bases. Every year a large amount of CO2 is generated in this basin. Therefore, it is essential to implement CGS in this basin for mitigating these CO2 emissions. The target layer used in this study is the sandstone formation located in the third and fourth Members of the Lower Cretaceous Quantou Formation (K1q3 and K1q4), which is a deep saline formation and oil‐producing unit, with a thickness of about 50 m. This formation is expected to be a potential candidate for CO2 geological storage.

In this study, based on the geological, hydrogeological and geochemical conditions of K1q3 and K1q4, the ultimate fate of injected CO2 and efficiencies of three trapping processes will be simulated with a multiphase, reactive transport simulator, the TOUGHREACT program.

1.3.2.2 Conceptual Model

Geometry and boundary conditions

A 2‐dimensional (2D) radial model was developed to investigate the transformation of injected CO2 within the target reservoir (Figure 1.1). The 2D model is a homogeneous sandstone formation of 50 m thickness with a cylindrical geometry. Ten layers were used vertically with a constant spacing of 5 m. In the horizontal direction, a distance of 100 km was modeled with a radial spacing that increases gradually away from the injection well. This distance was discretized into 45 radial grid elements. CO2 is injected through two elements at the bottom of the well (Figure 1.1) at a constant rate of 20 kg/s, and the total amount of CO2 is calculated to be 0.64 Mt/a. The injection lasts for 100 years, while the simulation time is 10,000 years.

Left: Schematic of 10 rings forming a cylinder measuring 50m. At the ring's center is a vertical bar. Right: A rectangle with a bar (at center) as injection well having shaded portion at its bottom.

Figure 1.1 Two‐dimensional radial schematic representation and discretization for CO2 injection model in the formation.

Physical parameters

Without consideration of the caprock for this reservoir, the model was treated as impermeable at the top and bottom boundaries. Considering that the model length of up to 100 km is out of pressure‐affecting range, lateral boundaries were assigned with fixed pressure. Hydrostatic pressure was imposed along the vertical layers.

The hydrogeological parameters used in the simulations are listed in Table 1.5. The whole formation was assumed to be homogeneous and isotropic. The petrophysical characteristics such as porosity and permeability were taken from the averaged experimental measurements and data collection (Zhang et al., 2009). Other parameters such as those used in the capillary pressure and relative permeability models were taken from Zhang et al. (2009) and are summarized in Table 1.5.

Table 1.5 Hydrogeological parameters used in the simulations.

Initial mineral composition

According to X‐ray diffraction (XRD) mineralogical analysis of the clastic rocks from the Songliao Basin, they are mainly composed of quartz, illite, chlorite, calcite, plagioclase and K‐feldspar. The composition of plagioclase is uncertain, because the ratio between albite and anorthite may affect the simulation results. In the work of Gaus et al. (2005), anorthite did have an impact on the mineral trapping of CO2. However, in the modeling of Audigane et al. (2007), all plagioclase is treated as pure albite because albite is much more resistant to weathering than anorthite, and Na will be the dominant cation in plagioclase. For the simplification in this study, plagioclase was represented using pure albite. The choice of secondary minerals may also affect the simulation results significantly. Almost all possible secondary carbonate and clay minerals are considered in the modeling, which refers to previous simulation studies (Xu et al., 2006; Gaus et al., 2005), laboratory research (Wolf et al., 2004), and field observations (Moore et al., 2005). The resulting mineralogical composition is presented in Table 1.6.

Table 1.6 Initial volume fractions of primary and secondary minerals used in the simulations.

Water chemistry

Prior to the simulations, a batch geochemical modeling of water–rock interaction was performed to obtain a quasi‐stable aqueous chemistry, where a 0.171 M (mol/kg H2O) solution of NaCl was used to react with the initially present minerals (listed in Table 1.6) at an initial CO2 partial pressure of 0.01 bar. After about 10 years the whole geochemical system reached an approximate steady state, resulting in a quasi‐stable solution composition (Table 1.7). This composition is used as the starting‐point water for the subsequent reactive transport modeling.

Table 1.7 Initial total dissolved component concentrations of the formation water at reservoir conditions of 50 °C and 120 bars.

Thermodynamic and kinetic parameters

The mineralogy in this study is described in terms of 16 minerals including primary minerals and secondary minerals, as shown in Table 1.6. The thermodynamic data for minerals, gases and aqueous species are mostly taken from the EQ3/6 V7.2b database of Wolery (1992). For kinetically controlled mineral dissolution and precipitation, a general form of rate law (Lasaga, 1984; Steefel and Lasaga, 1994) is used (Equation 1.4). The kinetic rate constants k(T) of concerned minerals are calculated using Equation (1.6). The values of the most relevant kinetic parameters are summarized in Appendix A to this chapter (Xu et al., 2011; Gherardi et al., 2007).

1.3.2.3 Results and Discussion

The injection of CO2 lasts 100 years with a rate of 20 kg/s (0.64 Mt/a), and the total reactive transport simulation was performed for a period of 10,000 years.

Changes in the distribution of gas saturation (Sg) with time

Once carbon dioxide is injected into the reservoir in a supercritical state (hereafter referred to as gas CO2), it will displace the ambient porewater, which will cause different regions with varying combinations of CO2 and brine. For about 100 m distance around the injection well, the pore spaces are completely filled with gas CO2 after 100 years. These zones are called the ‘dry area’, within which the CO2 gas saturations (abbreviated to Sg) are higher or at least not lower than 0.95. Beyond 100 m the Sg lies in the range from 0.4 to 0.6, where gas CO2 and porewater exist together. In the outermost zones Sg declines to almost zero.

Due to the relatively lower density of supercritical CO2 compared with brine, the buoyancy force moves supercritical CO2 upwards until it reaches the caprock. With the return of groundwater, a very small amount of the CO2 is immobilized and permanently trapped; this trapping process is referred to as residual sequestration. However, a large amount of CO2 will move towards and accumulate beneath the caprock, which gives rise to pressure under the base of the caprock. Locally higher pressure forces CO2 to migrate laterally.

At the beginning of injection, CO2 mainly distributes around the injection well, displacing groundwater. Sg of CO2 is able to reach 0.95. After injection (100 years), along with vertical and lateral propagation of CO2, the displaced groundwater returns, which causes a mutual displacement between supercritical CO2 and groundwater. At the time point of 5000 years, the space occupied by CO2 is small due to upward movement (Figure 1.2). At the end of the simulation, most of the injected CO2 accumulates beneath the caprock, and supercritical CO2 gradually dissolves into groundwater.

Elevation vs. distance depicting the spatial distributions of CO2 saturation after 5,000 (top) and 10,000 years (bottom). Each graph has various shades (light–dark) with the dark shade being more prominent.Elevation vs. distance depicting the spatial distributions of CO2 saturation after 5,000 (top) and 10,000 years (bottom). Each graph has various shades (light–dark) with the dark shade being more prominent.

Figure 1.2 Spatial distributions of CO2 saturation after (a) 100, (b) 1000, (c) 5000 and (d) 10,000 years.

Solubility trapping

Part of the injected supercritical CO2 dissolved into the formation water. This dissolved CO2 initiates the increase of aqueous solution density (density change depends on the salinity, temperature and pressure), consequently causing gravitational instability of the ‘gas’‐water interface. The dissolved CO2 migrates through molecular diffusion, diffusion, and convection. As the large amount of dissolved CO2 accumulates beneath the caprock, a high‐density solution saturated with dissolved CO2 moved downwards under gravity (Figure 1.3a, b and c), which was accompanied by upward migration of gaseous CO2. This phenomenon is so‐called ‘convective mixing’. At the end of the simulation at 10,000 years (Figure 1.3d), an obvious fingering phenomenon can be observed, and the maximum distance of dissolved CO2 migration reaches 9000 m beneath the bottom of the caprock.

Elevation vs. distance depicting spatial distributions of dissolved CO2 at 100 (top) and 1,000 years (bottom). Each graph has various shades with darker shade being more dominant. Top graph presents a very dark shade.Elevation vs. distance depicting the spatial distributions of dissolved CO2 at 5, 000 (top) and 10,000 years (bottom). Each graph has various shades with the darker shade being more dominant.

Figure 1.3 Spatial distribution of dissolved CO2 at (a) 100, (b) 1000, (c) 5000 and (d) 10,000 years.

These above‐mentioned processes triggered the solubility trapping mechanism of injected CO2. Subsequently, the dissolved CO2 could react with formation water forming carbonic acid (H2CO3), and then this acid would dissociate into proton (H+) and bicarbonate ion (HCO3−), resulting in the acidification of brine.

Variation of pH

As mentioned above, the injected CO2 acidified the formation water, and this process is reflected by the spatial and temporal distribution of pH. As depicted in Figure 1.4, during CO2 injection pH continued dropping. At the end of injection (Figure 1.4a), pH declined to about 4.5. Along with the increment of distance from injection well, owing to the consumption of dissolved CO2, pH of the formation water decreased. After injection stopped, CO2 was trapped in three forms: gaseous, aqueous and solid (trapped during mineral alteration). The zones with low pH also move upwards and laterally (Figure 1.4c,d). Spatial and temporal distribution of pH is subject to the dissolution of CO2. Therefore, the patterns of gas CO2 and dissolved CO2 are consistent.

Elevation vs. distance depicting the spatial distributions of pH at 100 (top) and 1000 years (bottom). Both graphs present various shades (light–dark), with the very dark shade being more prominent.Spatial changes in mineral volume fraction and porosity after 500 (top) and 1000 (bottom) years in the case with Sg=0.5, each displaying curves with markers for K-feldspar, oligoclase, chlorite, quartz, etc.

Figure 1.4 Spatial distribution of pH at (a) 100, (b) 1000, (c) 5000 and (d) 10,000 years.

Mineral alterations

The increased acidity (decreased pH) induces mineral dissolution/precipitation. In this study (the end point of simulation at 10,000 years was chosen to conduct discussion), spatial distributions of the volume fraction changes in minerals, which experienced dissolution/precipitation, are presented in Figures 1.5 and 1.6, respectively.

Spatial changes in mineral volume fraction and porosity after 100 (top) and 10 (bottom) years in the case with Sg=0.9, each displaying 2 curves with markers for Phi_Delta and calcite and a horizontal dashed line.Spatial changes in mineral volume fraction and porosity after 500 (top) and 1000 (bottom) years in the case with Sg=0.9, each displaying curves with markers for Phi_Delta, k-feldspar, illite, oligoclase, etc.

Figure 1.5 Spatial distribution of volume fraction change for dissolved minerals at 10,000 years (negative values indicate dissolution, positive precipitation).

Schematic representation of the 1D model displaying a shaded square labeled Matrix with solid and dashed lines depicted by arrows labeled Injection, Production, Fractured vein, 1m, and 600 m.Schematic of MINC method displaying a box containing intersecting lines with 3 rightward arrows pointing to a vertical bar with fracture, linked to 3 boxes labeled Matrix by rightward over leftward arrows.Change of volume fraction vs. distance displaying 2 ascending curves representing calcite (solid) and pH (dashed).

Figure 1.6 Spatial distribution of volume fraction change for precipitated minerals at 10,000 years.

From Figure 1.5 we can see that albite, calcite, K‐feldspar and chlorite dissolved significantly. A large amount of dissolved CO2 acidified the ambient water, reducing the pH. The equilibrium between calcite and formation water was broken, which resulted in the dissolution of calcite (Figure 1.5a). At the same time, plagioclase (albite) and K‐feldspar also underwent dramatic dissolutions (Figure 1.5b,c). Under acidic conditions, chlorite is vulnerable to dissolution. After injection, the chlorite began to dissolve, and the dissolution extent gradually decreased with the radial distance from the injection well (Figure 1.5d).

At the end of the simulation, maximum volume fraction change for albite was 0.07%, calcite 0.022%, K‐feldspar 0.015%, and chlorite 0.013%, respectively.

In contrast, due to release of ions from particular minerals, concentrations will increase until saturation of a specific ion‐bearing mineral is reached. As chlorite, calcite and plagioclase dissolve, sufficient Mg²+, Na+ and Ca²+ are released into the formation water to make the precipitation of smectite‐Na (Figure 1.6a), illite (Figure 1.6b), ankerite (Figure 1.6c), and smectite‐Ca (Figure 1.6e) possible

The redundant Mg²+ and concentrated HCO3− formed small amounts of magnesite (Figure 1.6 h). Additionally, sufficient Fe²+ and HCO3− facilitate the formation of siderite (Figure 1.6 g).

Due to the significant dissolution within the two‐phase mixing region, abundant ions required for precipitation of particular minerals were provided, which caused precipitation of minerals mainly within these regions. Magnesite, siderite and dawsonite deposit in the case of ample cations. This is the reason why these minerals precipitated later. By 10,000 years, the maximum precipitation of smectite‐Ca is 0.065% in volume fraction, illite 0.02%, ankerite 0.015%, quartz 0.0018%, and magnesite 0.0016%, respectively. Ankerite and dawsonite are the main minerals to solidify CO2, followed by magnesite and siderite.

The amount of CO2 trapped by minerals rose from the bottom to top of reservoir, which was caused by a large accumulation of CO2 at the top. At the end of the simulation the amount of CO2 trapped by minerals reached 8 kg/m³.

Transformation among three CO2 trapping forms over time

The injected CO2 was sequestered through three mechanisms in gaseous, liquid and solid forms, respectively. Figure 1.7 shows the amount of transformation of trapped CO2 between the three phases with time. At the beginning of the simulation, the injected CO2 was sequestered mainly by hydrodynamic trapping. With increasing simulation time, the CO2 dissolved into the formation water gradually. Due to the formation of carbonic acid, some susceptible minerals dissolved, supplying ions for the formation of carbonate minerals (dawsonite, magnesite, ankerite, siderite, etc.), resulting in mineral trapping of CO2. The amounts of CO2 trapped in different phases by the end of the simulation are summarized in Table 1.8.

Abundance of CO2 sequestration (Kg) vs. time (a) displaying 4 curves for the total (solid), gas (dashed), aqueous (dash-dotted), and mineral (short dash).

Figure 1.7 Changes in amounts of CO2 in three phases with time.

Table 1.8 Amount of CO2 trapped in the three phases at the end of the simulation.

1.3.2.4 Summary and Conclusions

When CO2 is injected into saline aquifers, it will undertake a series of physical and chemical processes for hundreds or thousands of years, or even longer. This processing mode of CO2 can affect the global carbon cycle. However, the deep saline aquifers in sedimentary basins are very complex, and processes involved in CO2 geological sequestration are complicated as well. Therefore, reactive fluid flow and geochemical transport numerical simulation is