The results classify each area into five levels: optimal, good, normal, poor and very poor. These values will help decision makers to evaluate which areas are the best considered and if it is reliable to go to the next stage.

**Figure 5.** CO2 storage site selection. A complex issue with different stakeholders

**CO2 source Maturity**

**Cultural resources (km)**

As described before, site selection is based on several criteria, values and weights. Even though the methodology proposed in this chapter allows selecting the best option in a quantitative and objective way, it is necessary to consider other points of view of the problem. Indeed, CO2 storage is a controversial way to remove CO2 from the atmosphere and even though it has been described a safe and affordable, there are many stakeholders who consider it unreliable. For those reason, and to manage the huge amount of technical information and different weight definitions of each criterion, a specific program has been developed: CO2SITE ASSESS. This software has been developed in VISUAL BASIC© (easy in terms of programming and speed in obtaining results, robust integration with Data Base, allowing operations in read/write formats). It includes the AHP algorithm (weights and values), so its interaction with the enduser is easy. Many of the technical and socioeconomic parameters can be represented in a Geographical Information System (GIS), so the CO2SITEASSESS results also generate a file

It is possible to differentiate two different databases: the first one comprises the CO2 emitters and its data (location, CO2 emission, primary energy, date of commissioning, and others), whereas the second database includes the CO2 storage location. Data to be included in this form should be the technical and socioeconomic criteria previously described, and the tool can compare the alternatives using the AHP algorithm and decision tree described in this chapter. The results classify each area into five levels: optimal, good, normal, poor and very poor. These values will help decision makers to evaluate which areas are the best considered and if it is

**Location**

>50 >20 >20 OnshoreModerate All None

**Climatology**

**Affected Infrastructure**

**New infrastructure**

**Quality of the information**

Detail data (GIS) based on deep data (wells and seismic)

**Distance (km)**

290 CO2 Sequestration and Valorization

< 25

reliable to go to the next stage.

**CO2 Quality**

Impurity content less than 1%

**4. CO2SITE ASSESS: Informatics tool**

**Population (km)**

**Table 2.** Socio-economical criteria and its best values for CO2 storage alternative.

which allows representing the results and *CO2 site storage assessment*.

**Environmental resources (km)**

**Figure 6.** CO2SITEASSES is a program to evaluate multiple criteria and to obtain georeferenced information related to potential CO2 store areas.

This tool is useful to compare many structures – even if there are many alternatives – and since the algorithm implemented in its code is based on AHP absolute mode, it is possible to compare more than seven alternatives.

## **5. Site selection: Evaluation of different options in southern Spain**

According to the Description of Work of the AVANZA CO2 project (national project supported by the Ministry of Industry and Tourism), the methodology proposed in this chapter has been applied to study specific geological structures. The company that supports this study (SACYR) considers bio-CCS technology as an option to decrease the CO2 emissions from its biomass plants located in the Guadalquivir basin (South of Spain). This study was carried out by the Technical University of Madrid, in collaboration with Gessal.

Guadalquivir basin: an area of potential interest

The Guadalquivir basin in southern Spain is an ENE–WSW elongated foreland basin devel‐ oped during the Neogene and Quaternary between the external zones of the Betic Cordillera to the south and Sierra Morena (Iberian Massif) to the north [24], which respectively forms its active and passive margins. The external zones of the Betic Foldchain are made up of Mesozoic and Cenozoic sediments that include thick calcareous and evaporitic formations, as well as siliciclastic units.

Sediments of the basin can be divided into two main stratigraphic: the lower, which includes materials deposited prior to the collision and embodies a long sedimentary process ranging from Cambrian to Permian, culminating in a strong tectonic collision known as the Hercynian orogeny; and the upper, comprising materials of the foreland basin itself, which is known as the alpine stage and begins with an intensive erosional period (Hercynian discordance), and a new subsequent sedimentary period that spreads from Upper Permian to Quaternary. The latter constitutes the proper filling of the foreland basin. It can be divided into five depositional sequences (relatively consistent set of strata, genetically related and whose roof and wall are discontinuity or continuity sequences).

It is possible to individualize two main sub-stages, which are disconnected by the alpine tectonic stage, of Burdigalian age. These sub-stages are:

In the same way as many other alpine forelands, compressive deformation seems to have been established by following a classic model of piggy back or progressive tectonic propagation, from the early active Southern System's Front to the Northern Passive Margin.

The selection of this area was made according to the following information:

**1.** The interest in the Guadalquivir Basin as a potential area for storing CO2 has been described in national and European Projects.

**2.** Technical conditions: Some structures were described in the Geocapacity Project, and others were proposed under the AVANZA CO2 Project. 1. The interest in the Guadalquivir Basin as a potential area for storing CO2 has been described in national and European Projects.

1. From Permian to Lower Miocene, sedimentation takes place over a passive or Atlantic type margin, which differentiates from North to South- platform, talus and deep water facies. These paleogeographic realms are known as Prebetic,

2. The Alpine orogeny reached its main deformation phase during the lower Miocene. We interpret this phase under a classical

In the same way as many other alpine forelands, compressive deformation seems to have been established by following a classic model

deformation model –intra continent subduction type, taking place under NW-SE and E-W compressional vectors.

Intermediate Unit, Subbetic and Flysch.

develop the bio-CCS concept.

The selection of this area was made according to the following information:

**5. Site selection: Evaluation of different options in southern Spain**

Technical University of Madrid, in collaboration with Gessal.

Guadalquivir basin: an area of potential interest

discontinuity or continuity sequences).

tectonic stage, of Burdigalian age. These sub-stages are:

described in national and European Projects.

taking place under NW-SE and E-W compressional vectors.

from the early active Southern System's Front to the Northern Passive Margin.

The selection of this area was made according to the following information:

siliciclastic units.

292 CO2 Sequestration and Valorization

and Flysch.

According to the Description of Work of the AVANZA CO2 project (national project supported by the Ministry of Industry and Tourism), the methodology proposed in this chapter has been applied to study specific geological structures. The company that supports this study (SACYR) considers bio-CCS technology as an option to decrease the CO2 emissions from its biomass plants located in the Guadalquivir basin (South of Spain). This study was carried out by the

The Guadalquivir basin in southern Spain is an ENE–WSW elongated foreland basin devel‐ oped during the Neogene and Quaternary between the external zones of the Betic Cordillera to the south and Sierra Morena (Iberian Massif) to the north [24], which respectively forms its active and passive margins. The external zones of the Betic Foldchain are made up of Mesozoic and Cenozoic sediments that include thick calcareous and evaporitic formations, as well as

Sediments of the basin can be divided into two main stratigraphic: the lower, which includes materials deposited prior to the collision and embodies a long sedimentary process ranging from Cambrian to Permian, culminating in a strong tectonic collision known as the Hercynian orogeny; and the upper, comprising materials of the foreland basin itself, which is known as the alpine stage and begins with an intensive erosional period (Hercynian discordance), and a new subsequent sedimentary period that spreads from Upper Permian to Quaternary. The latter constitutes the proper filling of the foreland basin. It can be divided into five depositional sequences (relatively consistent set of strata, genetically related and whose roof and wall are

It is possible to individualize two main sub-stages, which are disconnected by the alpine

**1.** From Permian to Lower Miocene, sedimentation takes place over a passive or Atlantic type margin, which differentiates - from North to South- platform, talus and deep water facies. These paleogeographic realms are known as Prebetic, Intermediate Unit, Subbetic

**2.** The Alpine orogeny reached its main deformation phase during the lower Miocene. We interpret this phase under a classical deformation model –intra continent subduction type,

In the same way as many other alpine forelands, compressive deformation seems to have been established by following a classic model of piggy back or progressive tectonic propagation,

**1.** The interest in the Guadalquivir Basin as a potential area for storing CO2 has been

**3.** Economic conditions: The specific area of interest is also defined due the interest of the company SACYR. This company has some power stations which use biomass as primary energy near this location. The interest of the company was evaluated in this area to develop the bio-CCS concept. 2. Technical conditions: Some structures were described in the Geocapacity Project, and others were proposed under the AVANZA CO2 Project. 3. Economic conditions: The specific area of interest is also defined due the interest of the company SACYR. This company has some power stations which use biomass as primary energy near this location. The interest of the company was evaluated in this area to

Figure 7. Detail of selected structures, and figures taken from previous studies by the IGME (Geocapacity project). **Figure 7.** Detail of selected structures, and figures taken from previous studies by the IGME (Geocapacity project).

preferred targets; both caprock and storage formations. The data include seismic reflection and refraction profiles, well logs, gravity and field observations. Finally, structural definition was done based on the interpretation of the geophysical data in each area. These interpretations allow us to define specific structures and define the CO2 capacity of each structure. Well interpretation was used to identify storage and caprock According to previous stratigraphic, petrological and petrophysical data obtained from exploration wells, it is possible to define preferred targets; both caprock and storage forma‐ tions. The data include seismic reflection and refraction profiles, well logs, gravity and field observations.

According to previous stratigraphic, petrological and petrophysical data obtained from exploration wells, it is possible to define

properties. **4.1. Sites evaluation: Application of the AHP model to this area**  Some of the structures considered were not evaluated for different conditions (shallow storage formations, lack of data or low thickness Finally, structural definition was done based on the interpretation of the geophysical data in each area. These interpretations allow us to define specific structures and define the CO2 capacity of each structure. Well interpretation was used to identify storage and caprock properties.

of the storage formations). These conditions should not eliminate these structures – indeed they have been included in the

#### CO2SITEASSESS data base – but social or economic aspects will cause them to be considered the worst areas to develop CO2 storage. **5.1. Sites evaluation: Application of the AHP model to this area**

For instance, A, B and C alternatives that have been considered. Even if the capacity calculated for each area is not enough for an industrial scale project, it could be considered for a pilot project or demonstration of the bio-CCS technology in Spain. Moreover, CO2SITEASSESS was used in another region (Duero Basin, also in Spain), where the assessment of this structures are much better than in the Guadalquivir basin. Some of the structures considered were not evaluated for different conditions (shallow storage formations, lack of data or low thickness of the storage formations). These conditions should not eliminate these structures – indeed they have been included in the CO2SITEASSESS data base – but social or economic aspects will cause them to be considered the worst areas to develop CO2 storage.

For instance, A, B and C are the alternatives that have been considered. Even if the capacity calculated for each area is not enough for an industrial scale project, it could be considered for a pilot project or demonstration of the bio-CCS technology in Spain.

Moreover, CO2SITEASSESS was used in another region (Duero Basin, also in Spain), where the assessment of this structures are much better than in the Guadalquivir basin.

**Table 3.** CO2SITEASSESS results.

All the areas defined in the present chapter have been previously defined by several hydro‐ carbon explorations: geophysical surveys and wells were considered: more than 10 wells were evaluated and hundreds of seismic studies were evaluated. Indeed, this region has active natural gas reservoirs – in two different turbidite systems: the Arenas del Guadiana Fm. in the Poseidón Gas field, in the Gulf of Cadiz, and the Arenas del Guadalquivir Fm., which produces from several small fields that are located onshore, in the Guadalquivir basin.

Moreover, outcrops were analyzed to properly evaluate mineral composition, hydrogeology and maturity conditions.

**Figure 8.** Spanish geological map and detailed area under evaluation. Source: IGME and AVANZA CO2 project.

**FUENSANTA (Alternative A):** Located in the basin's internal prebetic area (intermediate unit). As reservoir this study evaluated a carbonate rock belonging to Dogger – Lias, whereas a marls is considered as the caprock.

**ID BASIN VALUE NOTE** Duero El Gredal-Utrillas **1 POOR** Duero El Gredal-Bunt 7.26 GOOD Guadalquivir Fuensanta de Martos **1 POOR** Guadalquivir Guadalquivir H1 6.33 NORMAL Guadalquivir Nueva Carteya 6.54 NORMAL

All the areas defined in the present chapter have been previously defined by several hydro‐ carbon explorations: geophysical surveys and wells were considered: more than 10 wells were evaluated and hundreds of seismic studies were evaluated. Indeed, this region has active natural gas reservoirs – in two different turbidite systems: the Arenas del Guadiana Fm. in the Poseidón Gas field, in the Gulf of Cadiz, and the Arenas del Guadalquivir Fm., which produces

Moreover, outcrops were analyzed to properly evaluate mineral composition, hydrogeology

**Figure 8.** Spanish geological map and detailed area under evaluation. Source: IGME and AVANZA CO2 project.

from several small fields that are located onshore, in the Guadalquivir basin.

**Table 3.** CO2SITEASSESS results.

294 CO2 Sequestration and Valorization

and maturity conditions.

The anticlinal trap, from the data collected, is estimated to have a total area of 15 km2 of Structure A. The roof of the structure would be located at 1081 m.

Existence of water with low salt content has been confirmed (7 gr/l). Structure-Trap Anticline, elongated in EW direction, preserved under Subbetic materials in its western part, limited by a front thrust in its margin N, N-verging. structural closure to the "spill point" is around 400 milliseconds, which may involve around 700-800 meters.

**Figure 9.** Detailed description of the *Fuensanta structure* (called alternative A)

**GUADALQUIVIR H-1 (Alternative B):** Located within the basin internal prebetic. It is possible to define As Dogger oolitic carbonates as a storage reservoir, whereas the Malm marl may be considered as its caprock.

The trap is defined as a folding anticlinal; considering the data collected it is estimated to have a total area of 26 km2 . The roof of the store formation would be located at 1668 m. However, the target (corresponding to the Oolític Dogger *Jabalcuz* formation) has low porosity, between 2.25 and 4%.

For instance, it is possible to define a potential reservoir and a caprock.

**Figure 10.** Evaluation of alternative A, using the AHP model described above.

**Figure 11.** Detailed description of the *Guadalquivir H-1 structure* (called alternative B)

**Figure 12.** Evaluation of alternative B, using the AHP model described above.

**Figure 11.** Detailed description of the *Guadalquivir H-1 structure* (called alternative B)

**Figure 10.** Evaluation of alternative A, using the AHP model described above.

296 CO2 Sequestration and Valorization

**NUEVA CARTEYA - 1.** Within the pre-betic terminal basin, the carbonate reservoir rock (Dogger oolitic) was evaluated. Malm marls are considered to be caprock.

The trap is an anticline, with a total area of 30 km2 and estimated thickness of 160 m. The roof would be located at a depth of 1240 m.

Although a storage formation has been defined, it is necessary to study whether this structure is closed.

**Figure 13.** Detailed description of the *Guadalquivir H-1 structure* (called alternative C)

**Figure 14.** Evaluation of alternative C, using the AHP model described above.

#### **5.2. Socioeconomic values**

**Figure 13.** Detailed description of the *Guadalquivir H-1 structure* (called alternative C)

298 CO2 Sequestration and Valorization

**Figure 14.** Evaluation of alternative C, using the AHP model described above.

The area under evaluation contains different features in the same region. For instance, while the socioeconomic parameters can be considered similar for all of them, other parameters are different slightly different (i.e., distance from CO2 sources, storage area and the nearest town). It was possible to estimate the values marked in the Figure 15.

#### **Figure 15.** Evaluation of socioeconomic parameters

The area was explored during the twentieth century, so there is enough information to build a GIS and to define some of the structures (conceptual or static model). As shown in Figure 8, there are industrial CO2 sources a short distance away and the quality of the flue gases is sufficient. Some of the emitters are biomass power stations, but close to this region it might be possible to identify a larger emitter (Puente Nuevo power plant, close to Córdoba).

There are towns and cities close to each structure, but the topography can be considered favorable, and there are no environmentally protected areas close to the structures under evaluation.

**Figure 16.** Protected areas (environmental conditions: white and green areas), and towns (red dots) in the area of interest.

## **6. Conclusions**

Site selection for storing CO2 is a complex issue, especially when deep saline aquifers are under assessment. These geological structures used to be poorly characterized and the risk of unsuccessful geology exploration is high. For this reason, the Multicriteria Decision Tool can be used to evaluate related technical and socioeconomic data on different alternatives under consideration. In addition, there are different stakeholders with different points of view, so the decision maker needs to take these viewpoints into account.

AHP is the proposed multicriteria algorithm. It selects the best area in an objective way. Therefore, its use decreases the risk associated with the site selection phase, and it will easily show the strengths and weaknesses of the information or characteristics of the alternatives under study. Furthermore, it can help increase social acceptance by stakeholders.

An innovative program (CO2SITEASSESS) was been developed and validated, using some defined areas (at basin and regional scale) and structures (local scale). This software also allows obtaining georeferenced data; and the combination of both uses (georeferenced data and AHP algorithm) has never before been applied to select areas to store CO2. This combination has some advantages:

In addition, some of the criteria can change during the pre-injection phase (i.e., data available, other formations of interest, etc.) and this tool can be useful to consider at which stage each alternative is during a period of time (pre-injection: selection, characterization, static or dynamic model and engineering).

The results obtained in the High Guadalquivir basin suggest that this area is not suitable for CO2 storage on an industrial scale, but some of the structures considered in this chapter could be useful for pilot scale, especially if bio-CCS technology is applied.

Nevertheless, the CO2SITEASSES methodology has been demonstrated as robust to identify the best alternative under evaluation, and it reduces the inherent risk associated with geolog‐ ical explorations.

The AHP is applied in this study in an absolute mode, so it allows the assessment of limitless alternatives. For instance, this method and software can be useful as a standard in different regions (i.e., Spain or Europe).

Another version of the CO2SITEASSESS will be developed in the near future to relate site selection and a program characterization of each alternative. This characterization program will consider the three characterization phases: outcrop, geophysics and wells.

## **Author details**

**Figure 16.** Protected areas (environmental conditions: white and green areas), and towns (red dots) in the area of

Site selection for storing CO2 is a complex issue, especially when deep saline aquifers are under assessment. These geological structures used to be poorly characterized and the risk of unsuccessful geology exploration is high. For this reason, the Multicriteria Decision Tool can be used to evaluate related technical and socioeconomic data on different alternatives under consideration. In addition, there are different stakeholders with different points of view, so

AHP is the proposed multicriteria algorithm. It selects the best area in an objective way. Therefore, its use decreases the risk associated with the site selection phase, and it will easily show the strengths and weaknesses of the information or characteristics of the alternatives

An innovative program (CO2SITEASSESS) was been developed and validated, using some defined areas (at basin and regional scale) and structures (local scale). This software also allows obtaining georeferenced data; and the combination of both uses (georeferenced data and AHP

under study. Furthermore, it can help increase social acceptance by stakeholders.

the decision maker needs to take these viewpoints into account.

interest.

**6. Conclusions**

300 CO2 Sequestration and Valorization

B. Llamas1\*, M. Arribas1 , E. Hernandez2 and L.F. Mazadiego1

\*Address all correspondence to<EMAIL_ADDRESS>
1 Technical University of Madrid. Engineering School of Mines, Madrid, Spain

2 GESSAL, Madrid, Spain

## **References**

[13] Carlsson C., and Fullér R., Fuzzy multiple criteria decision making: Recent develop‐ ments, *Fuzzy Sets and Systems*, 1996, 78, 139-153.

**References**

302 CO2 Sequestration and Valorization

*for Global Change*, 10 (4), 693-715.

Mathiassen O.M., CO2

ards. *8* th

277-289.

United Kingdom, 2010, 27-56.

*criteria*. GEOCAPACITY Project (6th

Brown A., *D 3.1.3 DYNAMIS CO*<sup>2</sup>

Mathiassen O.M. CO2

Norway, 2006.

bridge University Press, United Kingdom. 2005.

[1] Torvanger, Asbjørn, Kristin Rypdal and Steffen Kallbekken, 2005. Geological CO2 Storage as a Climate Change Mitigation Option. *Mitigation and Adaptation Strategies*

[2] Benson, S, Cook P, Metz B. (Ed.), Davidson O. (Ed.), de Coninck H. (Ed.), Loos M. (Ed.), Meyer L. (Ed.) *IPCC Special Report: Carbon Dioxide Capture and Storage*. Cam‐

[3] Bradshaw J., Bachu S., Bonijoly D., Burruss R., Holloway S., Christensen N.P. and

[4] M. Carpenter, K. Kvien, J. Aarnes. The CO2QUALSTORE guideline for selection, characterisation and qualification of sites and projects for geological storage of CO2

[5] Bachu S., Sequestration of carbon dioxide in geological media: criteria approach for

[6] Bachu S., Screening and ranking of sedimentary basins for sequestration of CO2 in geological media in response to climate change. *Environmental Geology*, 2003, 44(3),

*Storage Technology, vol. 2*, Woodhead Energy Series 16. Woodhead Publishing Ltd,

[8] Benson S., Hepple R., Apps J., Tsang, C-F, Lippmann M., *Lessons learned from natural and industrial analogues for storage of carbon dioxide in deep geological formations*, Law‐

[9] European Parliament: *Directive 2009/31/EC on the geological storage of carbon dioxide*.

[10] Vangkilde-Pedersen, T. et al (2009): *D26 WP4 report: capacity standards and site selection*

[11] Bachu S., Bonijoly D., Bradshaw J., Burruss R., Christensen N.P., Holloway S. and

[12] De Visser E., Hendriks C., de Koeijer G., Liljemark S., Barrio M., Austegard A.,

Framework Programme). 2009

storage capacity estimation: methodology and gaps, *Interna‐*

 *quality recommendations*. DYNAMIS Project (6th

*International Journal of Greenhouse Gas Control,* 2011, 5, 942–951.

[7] Bachu S., Screening and selection criteria, and characterization for CO2

Marotovaler M (Ed.) *Developments and Innovation in Carbon Dioxide (CO*<sup>2</sup>

site selection. *Energy Conv. Manage*, 2000, 41, 953-970.

rence Berkeley National Laboratory, United States, 2002.

*tional Journal of Greenhouse Gas Control,* 2007, 1, 304-443.

Framework Programme). The Netherlands, 2007, 16-35.

Official Journal of the European Union L140/114. 2009, Brussels.

 *International Conference on Greenhouse Gas Control Technology*. Trondheim,

storage capacity estimation: Issues and development of stand‐

.

storage. In:

*) Capture and*

## **Numerical Simulation of CO2 Sequestration in Large Saline Aquifers**

Zheming Zhang and Ramesh K. Agarwal

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/57065

## **1. Introduction**

With heightened concerns over CO2 emissions from pulverized-coal power plants, there has been major emphasis in recent years on the development of safe and economical geological carbon sequestration (GCS) technology. Although it is one of the most promising technologies to address globalwarming due to anthropogenic CO2 emissions, the detailed mechanisms of GCS are not well understood. As a result, there remain many uncertainties in determining the sequestration capacity of the formation/reservoir and the safety of sequestered CO2 due to leakage. These uncertainties arise due to lack of information about the detailed interior geometry of the formation and the heterogeneity in its geological properties, such as perme‐ ability and porosity, which influence the sequestration capacity and plume migration. Furthermore, the sequestration efficiency is highly dependent on the injection strategy, which includes injection rate, injection pressure, type of injection well employed and its orientation etc. The goal of GCS is to maximize the sequestration capacity and minimize the plume migration by optimizing the GCS operation before proceeding with its large-scale deployment. In this chapter, numerical simulations of GCS are conducted using the US Department of Energy (DOE) multi-phase flow solver TOUGH2 (Transport of Unsaturated Groundwater and Heat). A multi-objective optimization code based on genetic algorithm is also developed to optimize the GCS operation for a given geological formation. It is described in Chapter titled," Optimization of CO2 Sequestration Saline Aquifers". Most of the studies are conducted for sequestration in a saline formation (aquifer). Large-scale GCS studies are conducted for the Mt.Simon, Frio and Utsira saline formations, for which some experimental data and compu‐ tations performed by other investigators are available. These simulation studies provide important insights as to the key sources of uncertainties that can influence the accuracy of simulations.

© 2014 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **2. Saline aquifer geological carbon sequestration (SAGCS)**

Studies of GCS have suggested that various geological structures can serve as potential CO2 storage sites. The major geological carbon sinks include the following structures: (1) conven‐ tional hydrocarbon reservoirs, (2) un-minable coal seams, (3) mature oil/gas reservoirs, and (4) deep saline formations. Among these candidates, this chapter focuses on carbon seques‐ tration in saline aquifers considering the following facts.

A huge geographic distribution of deep saline aquifers in North America has been identified by the DOE. The DOE-estimated storage capacity for SAGCS takes into account more than 80 percent of the overall storage capacity of all possible GCS sites. In Table 1, the low-end capacity of 2102 billion metric ton of CO2 is estimated under the condition that ineffective storage may occur due to improper and non-optimized sequestration approaches. On the other hand, the high-end capacity of 20043 billion metric ton of CO2 is estimated under the conditions that most effective and optimal storage takes place. It can be seen that the high-end estimated capacity is over nine times the low-end estimated capacity. This large difference means it is important to deploy optimized reservoir engineering techniques for effective utilization of storage potential and successful GCS practice.

## **3. The approach for numerical simulations of GCS**

The large spatial extent (of the order of kilometers) and time duration (centuries) for CO2 plume migration after injection makes the study of SAGCS very difficult at these large spatial and temporal macro-scales by using laboratory scale experiments, which can be conducted only on relatively small spatial and temporal scales varying from nanometers to a few meters and from nanoseconds to a few days/months. Conducting field tests in large-scale formations before the actual deployment takes place can be very expensive. However, numerical simu‐ lations using computation fluid dynamics (CFD) technology can be employed at industrial

**2. Saline aquifer geological carbon sequestration (SAGCS)**

tration in saline aquifers considering the following facts.

possibly suitable for GCS in the US and Canada.

of US overall GCS potential [1].

storage potential and successful GCS practice.

**3. The approach for numerical simulations of GCS**

saline aquifers.

306 CO2 Sequestration and Valorization

Studies of GCS have suggested that various geological structures can serve as potential CO2 storage sites. The major geological carbon sinks include the following structures: (1) conven‐ tional hydrocarbon reservoirs, (2) un-minable coal seams, (3) mature oil/gas reservoirs, and (4) deep saline formations. Among these candidates, this chapter focuses on carbon seques‐

**•** Concentrated locations of major sources of CO2 (such as power plants) are close to existing

**•** Geological survey has confirmed vast geological distribution of deep saline formations

**•** Preliminary estimates have suggested large storage capacity of the existing deep saline formations. The US DOE estimates an aggregate storage capacity of approximately 2102~20043 billion metric tons of CO2 for SAGCS in US, which accounts for 80~90 percent

**•** Since most of the saline formations are located deep underground, i.e., at least 800 m below

**•** Many surveys, research projects and commercial projects have already been conducted for

A huge geographic distribution of deep saline aquifers in North America has been identified by the DOE. The DOE-estimated storage capacity for SAGCS takes into account more than 80 percent of the overall storage capacity of all possible GCS sites. In Table 1, the low-end capacity of 2102 billion metric ton of CO2 is estimated under the condition that ineffective storage may occur due to improper and non-optimized sequestration approaches. On the other hand, the high-end capacity of 20043 billion metric ton of CO2 is estimated under the conditions that most effective and optimal storage takes place. It can be seen that the high-end estimated capacity is over nine times the low-end estimated capacity. This large difference means it is important to deploy optimized reservoir engineering techniques for effective utilization of

The large spatial extent (of the order of kilometers) and time duration (centuries) for CO2 plume migration after injection makes the study of SAGCS very difficult at these large spatial and temporal macro-scales by using laboratory scale experiments, which can be conducted only on relatively small spatial and temporal scales varying from nanometers to a few meters and from nanoseconds to a few days/months. Conducting field tests in large-scale formations before the actual deployment takes place can be very expensive. However, numerical simu‐ lations using computation fluid dynamics (CFD) technology can be employed at industrial

sea level, they provide great potential for secure long-term sequestration.

SAGCS, making it attractive for further research and technical contributions.

**Table 1.** Estimation of saline aquifer storage capacity of GCS for different regional carbon sequestration partnerships (RCSPs) [1]

scale to determine the fate of injected CO2 in a reservoir. With the development over the past four decades, CFD technology has now become mature and has been widely and successfully applied to various engineering problems. With the proper modeling of the storage formation and ground water transportation, CFD is capable of providing accurate enough analysis for quick estimation of reservoir performance at considerably lower cost. The governing equations of mass/energy transportation and numerical representations of the formation properties have been well explained in the TOUGH2 User's Guide [2],[3].

In a complex simulation like that of SAGCS, it is impractical to integrate all geophysical and geochemical effects into a single model while retaining acceptable computational efficiency. Therefore, careful examination of physical phenomenon of interest in SAGCS is essential to determine simplifications in modeling of features of interest [4],[5]. Another important benefit of numerical simulations is that one can investigate the effect of various injection parameters such as injection rate, injection duration, and injection well orientation and displacement on CO2 storage efficiency and plume migration in a given reservoir. The advantage of numerical simulations makes it possible to perform optimization studies of these injection parameters to achieve the highest possible storage efficiency and minimum plume migration, as described in Chapter titled," Optimization of CO2 Sequestration Saline Aquifers". Such an optimization capability can be very beneficial in successful cost effective implementation of SAGCS on industrial scale. We have employed genetic algorithm in conjunction with numerical simulator TOUGH2 for optimization of SAGCS practice. Some salient features of the genetic algorithm and its integration with TOUGH2 are described in Chapter titled," Optimization of CO2 Sequestration Saline Aquifers".

## **4. Simulation code validation using analytical and benchmark solutions**

TOUGH2 is available as source files written in Fortran77. TOUGH2 does not provide graphical user interface (GUI) of any kind. All its input files and output results are in ASCII format. TOUGH2 has very high computing efficiency when executing large-scale simulations and it is very convenient for users to make modifications to the source code if needed. However, both the problem setup and result analysis capability in TOUGH2, such as mesh generation and contour map visualization, are not as comprehensive or straightforward as available in some newer commercial multiphase flow field simulators. For any complex problem, the modeling process tends to be tedious and error-prone. To address such deficiency, a third-party GUI for TOUGH2 named PetraSim was also installed on the same machine with TOUGH2. PetraSim preserves the original TOUGH2 binary files to execute simulations, while providing a smooth interface to a user-friendly computing environment. However, PetraSim in its original form lacks the compatibility for integrating a new optimization module and some recently devel‐ oped equation-of-state modules.

Previous validation simulations on benchmark problems have shown that the simulation results obtained by TOUGH2 and PetraSim are identical [6],[7]. For our code validation purpose, we employed three widely used benchmark problems by GCS researchers world‐ wide. Simulations were conducted by both TOUGH2 and PetraSim. These three benchmark problems were first defined in the Workshop on Numerical Models for Carbon Dioxide Storage in Geological Formations at the University of Stuttgart, Germany [8],[9],[10],[11],[12]. We study the benchmark problem #1 and #3 using PetraSim, while benchmark problem #2 is simulated using the original version of TOUGH2 because of the limit on the availability of the needed equation-of-state module in PetraSim. When simulation is performed using TOUGH2, necessary post-processing programs such as Tecplot are employed for visualization and analysis.

#### **4.1. Simulation of in situ CO2 migration and comparison with analytical solution**

As a first step towards code validation, simulations for CO2 plume migration in an ideal simplified reservoir are performed. The analytical solution for this case is available [13], [14], [15] and is obtained as:

$$\frac{b\_{\rm ag}(r,t)}{B} = \frac{1}{\mu\_{\rm b} - \mu\_{\rm ag}} \left[ \sqrt{\frac{\mu\_{\rm b} \mu\_{\rm ag} V(t)}{B \pi \phi r^2}} - \mu\_{\rm ag} \right] \tag{1}$$

where *t* is the time elapsed since the inception of injection, *r* is the distance (radius) from the injection well, *bag(r,t)* is the plume thickness as a function of *r* and *t*, *B* is the total thickness of the reservoir, *φ* is the porosity of the reservoir, *µ* is the dynamic viscosity, subscripts *ag* and *b* stand for the injected CO2 and brine respectively, and *V* (*t*)=*∫Qdt* is the volume of the injected CO2 within time *t*. For a horizontal reservoir, setting *bag(r,t)* = 0 yields Eq. (2), which gives a quick evaluation of the maximum CO2 plume migration as:

$$R\_{\text{max}} = \sqrt{\frac{\mu\_b V(t)}{\mu\_{ag} B \pi \phi}}\tag{2}$$

In numerical simulations, a hypothetical deep saline reservoir of thickness of 100 m is assumed. A cylindrical computational domain is considered as shown in Figure 1. Generic hydrogeological properties are used. CO2 injection rate is set at 1 kg per year for ten years. Detailed model parameters used in the simulations are summarized in Table 1. CO2 plume migration in each year is computed by the simulation and compared with the analytical solution given by Eq. (2).

**Figure 1.** Computational domain and mesh of a generic cylindrical aquifer

is very convenient for users to make modifications to the source code if needed. However, both the problem setup and result analysis capability in TOUGH2, such as mesh generation and contour map visualization, are not as comprehensive or straightforward as available in some newer commercial multiphase flow field simulators. For any complex problem, the modeling process tends to be tedious and error-prone. To address such deficiency, a third-party GUI for TOUGH2 named PetraSim was also installed on the same machine with TOUGH2. PetraSim preserves the original TOUGH2 binary files to execute simulations, while providing a smooth interface to a user-friendly computing environment. However, PetraSim in its original form lacks the compatibility for integrating a new optimization module and some recently devel‐

Previous validation simulations on benchmark problems have shown that the simulation results obtained by TOUGH2 and PetraSim are identical [6],[7]. For our code validation purpose, we employed three widely used benchmark problems by GCS researchers world‐ wide. Simulations were conducted by both TOUGH2 and PetraSim. These three benchmark problems were first defined in the Workshop on Numerical Models for Carbon Dioxide Storage in Geological Formations at the University of Stuttgart, Germany [8],[9],[10],[11],[12]. We study the benchmark problem #1 and #3 using PetraSim, while benchmark problem #2 is simulated using the original version of TOUGH2 because of the limit on the availability of the needed equation-of-state module in PetraSim. When simulation is performed using TOUGH2, necessary post-processing programs such as Tecplot are employed for visualization and

**4.1. Simulation of in situ CO2 migration and comparison with analytical solution**

(,) 1 ( ) *ag b ag*

*b rt V t B B r*

*b ag*

max

*V t <sup>R</sup>*

m pf

m m

quick evaluation of the maximum CO2 plume migration as:

As a first step towards code validation, simulations for CO2 plume migration in an ideal simplified reservoir are performed. The analytical solution for this case is available [13], [14],

> m m

<sup>=</sup> ê ú - - ë û

where *t* is the time elapsed since the inception of injection, *r* is the distance (radius) from the injection well, *bag(r,t)* is the plume thickness as a function of *r* and *t*, *B* is the total thickness of the reservoir, *φ* is the porosity of the reservoir, *µ* is the dynamic viscosity, subscripts *ag* and *b*

stand for the injected CO2 and brine respectively, and *V* (*t*)=*∫Qdt* is the volume of the injected CO2 within time *t*. For a horizontal reservoir, setting *bag(r,t)* = 0 yields Eq. (2), which gives a

> ( ) *<sup>b</sup> ag*

*B* m

2

é ù

pf

*ag*

<sup>=</sup> (2)

(1)

m

oped equation-of-state modules.

308 CO2 Sequestration and Valorization

analysis.

[15] and is obtained as:

**Table 2.** Geometry parameters and hydro-geological properties of the generic saline aquifer

The simulation time is ten years and the CO2 migration within the aquifer is computed for each of the ten years. In Figure 2, the CO2 plume after 1, 4, 7 and 10 years since the inception of injection is shown.

As seen in Figure 2, the injected CO2 migrates upwards very rapidly and then prominently migrates underneath the caprock. A typical plume shape is already identifiable after one year

**Figure 2.** CO2 plume after 1, 4, 7 and 10 years after injection

of injection, when the farthest CO2 migration reaches about 100 m from the injection well. In the following nine years, in situ CO2 keeps migrating outwards and spreads to 300 m after the tenth year of injection. Physically, such a large radial migration of in situ CO2 is caused by gravity separation. The horizontal extent of the plume can be analytically calculated by utilizing Eq. (2). Taking the necessary values of reservoir and fluid properties from Table 2, the horizontal extent of the plume predicted by Eq. (2) for the first 10 years since injection is summarized in Table 3. The horizontal extent of the plume given by TOUGH2 simulations is also summarized in Table 3 for comparison with the analytical solution.

**Table 3.** Maximum CO2 migration underneath the caprock given by the analytical solution and TOUGH2 simulation

As shown in Table 3, TOUGH2 simulations successfully predict the maximum CO2 plume migration underneath the caprock with excellent agreement with the analytical solution given by Bachu, Nordbotten and Celia [13],[14],[15]. The insignificant difference between the numerical and analytical solutions can be explained by the fact that CO2 dissolution is accounted for in TOUGH2 while it is neglected in the derivation of the analytical solution. Since CO2 dissolution is governed by the contact area between CO2 and the ambient brine, the rate of CO2 dissolution into brine should gradually increase over time as larger contact area becomes available as the CO2 plume spreads. Nevertheless, Table 3 validates TOUGH2 as an accurate simulation tool for predicting the migration of in situ CO2.

A schematic of the CO2 plume flow is shown in Figure 3. Although based on a simplified analytical model, Figure 3 shows in situ CO2 migration due to the combined pressure-driven Darcy flow and the buoyancy-drive CO2 transport.With the injection well located on the left side of Figure 3, the CO2 plume can be identified as consisting of two distinct regions. The first region is a smaller region on the left, adjacent to the injection well, marked as region (1) in Figure 3. In this region, CO2 is distributed uniformly through the entire period of the injection interval. This implies strong hydrodynamic force caused by the pressure difference between the pressurized injection well and the unaffected aquifer. Within this region, lateral pressure gradient dominates the movement of CO2 and Darcy flow occurs, causing CO2 to migrate more radially through the aquifer. The second region is marked as region (2) in Figure 3, where the CO2 plume fully develops. In this region, buoyancy due to the density difference between CO2 and brine becomes dominant and drives the upward movement of CO2 along with lateral migration. In this region, the vertical movement of CO2 becomes dominant and results in plume flow. Since it is a phenomenon of fundamental concern in SAGCS, understanding the development of plume flow is critical for the success of SAGCS. The size of the two regions in Figure 3 can vary depending on the properties of the actual aquifer, but under most conditions, region (2) becomes dominant in size, which influences the safety and efficiency of SAGCS operations. Therefore, every effort should be made either to increase the size of region (1) or decrease the size of region (2) for successful and desirable implementation of SAGCS.

**Figure 3.** Schematic of the shape of in situ CO2 plume

of injection, when the farthest CO2 migration reaches about 100 m from the injection well. In the following nine years, in situ CO2 keeps migrating outwards and spreads to 300 m after the tenth year of injection. Physically, such a large radial migration of in situ CO2 is caused by gravity separation. The horizontal extent of the plume can be analytically calculated by utilizing Eq. (2). Taking the necessary values of reservoir and fluid properties from Table 2, the horizontal extent of the plume predicted by Eq. (2) for the first 10 years since injection is summarized in Table 3.  Maximum CO2 migration underneath the caprock given by the analytical solution and TOUGH2 simulation

**Maximum Migration based on Analytical Solution (B)**

**Deviation ((A-B)/B)**

also summarized in Table 3 for comparison with the analytical solution.

**Year Maximum Migration based on**

**Numerical Simulation (A)**

**Figure 2.** CO2 plume after 1, 4, 7 and 10 years after injection

310 CO2 Sequestration and Valorization

#### **4.2. Simulation of benchmark problem #1 – CO2 plume evolution and leakage through an abandoned well**

A three-layered formation is modeled for the first benchmark problem [9]. CO2 is injected into the deeper aquifer, shown schematically in Figure 4. It spreads in the aquifer and then rises up to a shallower aquifer upon reaching a leaky well. Quantification of the leakage rate, which depends on CO2 plume evolution and the pressure buildup in the aquifer, is the main objective of this benchmark simulation. Figure 4 shows the schematic of the problem description by providing a cross-section of the formation.

**Figure 4.** Schematic of benchmark problem #1 (cross-sectional view) [9]

The three layers in Figure 4 are identified as one aquitard layer and two (saline) aquifer layers. The lower aquifer layer is assumed to be 3000m below the ground surface. Typical saline aquifer conditions and hydrogeological properties, such as temperature, salinity, permeability, are assigned to the aquifer layers. The aquitard is assumed to be impermeable to both saline and CO2, so it is considered an ideal geological seal to flow transportation. An "abandoned well" fully penetrating the three layers is located 100 m away from the CO2 injection well. It can be either a crack in the formation or a physical abandoned well, which served as a pathway for upward CO2 migration. Supercritical CO2 is injected only into the lower aquifer through the injection well. Being less dense than brine, injected supercritical CO2 gradually migrates to the ceiling of the lower aquifer and forms a plume. The formation and migration of the plume depends on the aquifer's geometric and hydrogeological properties. Table 4 summa‐ rizes the geometric properties of the aquifer in benchmark problem #1. It should be noted that the actual geometry of the injection well and abandoned well is circular, with a radius of 0.15 m. Since the use of an unstructured grid is not supported by PetraSim, an approximation to the circular geometry is made. Maintaining an identical cross-sectional area, the original circular injection well and the leaky well are replaced by wells of square cross-section with dimension of 0.266 m × 0.266 m. Such an approximation is acceptable since the details of the flow pattern inside the wells are not critical in achieving the objective of this simulation.

**Table 4.** Geometry parameters for benchmark problem #1

**4.2. Simulation of benchmark problem #1 – CO2 plume evolution and leakage through an**

A three-layered formation is modeled for the first benchmark problem [9]. CO2 is injected into the deeper aquifer, shown schematically in Figure 4. It spreads in the aquifer and then rises up to a shallower aquifer upon reaching a leaky well. Quantification of the leakage rate, which depends on CO2 plume evolution and the pressure buildup in the aquifer, is the main objective of this benchmark simulation. Figure 4 shows the schematic of the problem description by

The three layers in Figure 4 are identified as one aquitard layer and two (saline) aquifer layers. The lower aquifer layer is assumed to be 3000m below the ground surface. Typical saline aquifer conditions and hydrogeological properties, such as temperature, salinity, permeability, are assigned to the aquifer layers. The aquitard is assumed to be impermeable to both saline and CO2, so it is considered an ideal geological seal to flow transportation. An "abandoned well" fully penetrating the three layers is located 100 m away from the CO2 injection well. It can be either a crack in the formation or a physical abandoned well, which served as a pathway for upward CO2 migration. Supercritical CO2 is injected only into the lower aquifer through the injection well. Being less dense than brine, injected supercritical CO2 gradually migrates to the ceiling of the lower aquifer and forms a plume. The formation and migration of the plume depends on the aquifer's geometric and hydrogeological properties. Table 4 summa‐ rizes the geometric properties of the aquifer in benchmark problem #1. It should be noted that the actual geometry of the injection well and abandoned well is circular, with a radius of 0.15

**abandoned well**

312 CO2 Sequestration and Valorization

providing a cross-section of the formation.

**Figure 4.** Schematic of benchmark problem #1 (cross-sectional view) [9]

#### Figure 5 shows the simulation domain and the structured mesh inside the domain.

**Figure 5.** Entire computational domain (left) and the zoomed-in-view (right) for benchmark problem #1

To accurately model the small dimensions of the wells and accurately capture the CO2 leakage rate, the mesh is highly refined in the neighborhood of the injection and leakage wells, as can be seen in the zoomed-in-view in Figure 5. Since high CO2 concentration is expected at the ceiling of the lower aquifer due to gravity segregation, the mesh in this part of the lower aquifer is also refined. The mesh in the upper aquifer is not refined since it does not affect the accuracy of simulations and results in less computational efforts. The upper aquifer is uniformly discretized with vertical discretization length of 10 m, since it is assumed that the leakage amount is small and the leakage plume shape is not of great interest. By establishing the simulation domain and the mesh in this manner, reasonably accurate results are obtained while keeping the computational effort relatively low. The hydrogeological properties of the simulation domain are summarized in Table 5.

**Table 5.** Hydrogeological parameters for benchmark problem #1

Other simulation parameters such as initial conditions and boundary conditions are summar‐ ized in Table 6.

**Table 6.** Simulation parameters for benchmark problem #1

With the properties and parameters summarized above, the numerical model of benchmark problem #1 is setup in PetraSim. A pre-injection simulation is carried out first with no injection of CO2 to achieve equilibrium condition under gravity. The equilibrium state is then imple‐ mented as an initial condition for the subsequent simulation with CO2 injection. The equili‐ brium simulation is critical to provide the simulation with CO2 injection with realistic initial conditions; this is a prerequisite procedure for all the simulations reported in this chapter. For this benchmark problem, it takes about five minutes of CPU time on the workstation for the simulation to finish. The leakage flux, pressure perturbation, and CO2 saturation distribution throughout the aquifer after 80 days of CO2 injection are examined and compared with the simulations of other investigators [12]. The leakage flux is a non-dimensional quantity defined as the ratio of CO2 leakage rate to CO2 injection rate. Detailed comparisons using various simulation codes are shown in Figure 6 and are summarized in Table 7.

**Figure 6.** CO2 leakage flux value obtained with WUSTL-TOUGH2 and other simulation codes

Aquifer Permeability 20 mDarcy Leaky Path Permeability 1 Darcy Porosity 0.15

Residual Brine Saturation 0.2 Residual CO2 Saturation 0.05 Relative Permeability Linear (*kr* = *S*) Capillary Pressure Brooks-Corey Entry Pressure 1.0 × 105 Pa

Brooks-Corey Parameter 2.0

**Thermal Condition Isothermal**

Initial CO2 Mass Fraction *XCO2* = 0 Initial Salt Mass Fraction *Xsm* = 0.20 Injection Rate 8.87 kg/s Simulation End Time 1000 days

simulation codes are shown in Figure 6 and are summarized in Table 7.

Other simulation parameters such as initial conditions and boundary conditions are summar‐

Initial Condition (Temperature) Geothermal gradient: 0.03 K/m, Initial value at 800 m: 34oC Initial Condition (Pressure) Pressure gradient: 1045 Pa/m, 30.86 MPa at 3000 m depth

With the properties and parameters summarized above, the numerical model of benchmark problem #1 is setup in PetraSim. A pre-injection simulation is carried out first with no injection of CO2 to achieve equilibrium condition under gravity. The equilibrium state is then imple‐ mented as an initial condition for the subsequent simulation with CO2 injection. The equili‐ brium simulation is critical to provide the simulation with CO2 injection with realistic initial conditions; this is a prerequisite procedure for all the simulations reported in this chapter. For this benchmark problem, it takes about five minutes of CPU time on the workstation for the simulation to finish. The leakage flux, pressure perturbation, and CO2 saturation distribution throughout the aquifer after 80 days of CO2 injection are examined and compared with the simulations of other investigators [12]. The leakage flux is a non-dimensional quantity defined as the ratio of CO2 leakage rate to CO2 injection rate. Detailed comparisons using various

Fixed-state on lateral boundaries No mass flow on top and bottom boundaries

**Table 5.** Hydrogeological parameters for benchmark problem #1

ized in Table 6.

314 CO2 Sequestration and Valorization

Boundary Conditions

**Table 6.** Simulation parameters for benchmark problem #1

In Figure 6, our result (WUSTL-TOUGH2) is shown by the large graph, while comparisons with other simulation codes are shown in the inner box.

**Table 7.** Simulation results and comparisons with other codes for benchmark problem #1

As additional comparisons, the pressure perturbation and CO2 saturation distribution after 80 days of injection are also computed and compared with those from the MUFTE numerical solver [12]. Excellent agreement is obtained, as shown in Figure 7 and Figure 8.

**Figure 7.** Pressure perturbation within the aquifer after 80 days of injection (left: WUSTL-TOUGH2; right: MUFTE)

**Figure 8.** CO2 distribution within the aquifer after 80 days of injection (left: WUSTL-TOUGH2; right: MUFTE)

As seen in Figure 8, the CO2 plume is nicely captured in the simulation.

The simulation of benchmark problem #1 is very instructive. Three conclusions can be made: 1) Small variations among the results from different numerical simulators with different users are un-avoidable. Such variations are expected because some parameters are intentionally left unspecified. 2) Our results are in satisfactory agreement with the results of other investigators. 3) The most important CO2 behavior under reservoir condition, i.e., the plume flow, is well captured and understood by the simulations. This simulation and others for benchmark problems #2 and #3 not only validate our numerical solver but also provide insights needed for the optimization studies reported in Chapter titled," Optimization of CO2 Sequestration Saline Aquifers".

#### **4.3. Simulation of benchmark problem #2 – enhanced CH4 recovery in combination with CO2 sequestration in depleted gas reservoirs**

For decades, the oil and gas industry has been using a reservoir engineering technique to increase the oil/gas production from mature reservoirs, known as enhanced oil/gas recovery (EOR/EGR). As the original formation fluid (oil or natural gas) gets extracted, pressure in the reservoir gradually decreases. Such de-pressurization process makes it increasingly difficult to maintain the desired production rate. The reservoir needs to be re-pressurized to mitigate the drop of oil/gas production. One of the means to do this is to inject CO2 into the mature reservoir. With void space being occupied by the injected CO2, remaining oil/gas is pushed out of the reservoir. Meanwhile, the depleted reservoir becomes an ideal carbon sink for longterm storage. EOR/EGR with CO2 sequestration, also known as CSEOR/CSEGR, has been frequently used by the industry due to its strong economic merits.

In benchmark problem #2, a five-spot pattern domainis considered for modeling. The five-spot pattern is a common configuration for oil/gas production. A schematic of the reservoir is shown in Figure 3.

Natural gas is produced at the four upper corners of the reservoir, while CO2 is injected in the middle of the domain at the bottom-most part. This is a direct result of CO2 being heavier than CH4 under the reservoir conditions. Injection of CO2 at the bottom avoids gas mixing and creates better sweep efficiency. The main goal of this benchmark simulation is to identify the gas recovery factor, defined as the ratio of enhanced CH4 production to the original remaining CH4 amount until the shutdown of the production well. Additionally, the time until production Numerical Simulation of CO2 Sequestration in Large Saline Aquifers http://dx.doi.org/10.5772/57065 317

**Figure 9.** Schematic of the 3D five-spot pattern for benchmark problem #2

**Figure 8.** CO2 distribution within the aquifer after 80 days of injection (left: WUSTL-TOUGH2; right: MUFTE)

The simulation of benchmark problem #1 is very instructive. Three conclusions can be made: 1) Small variations among the results from different numerical simulators with different users are un-avoidable. Such variations are expected because some parameters are intentionally left unspecified. 2) Our results are in satisfactory agreement with the results of other investigators. 3) The most important CO2 behavior under reservoir condition, i.e., the plume flow, is well captured and understood by the simulations. This simulation and others for benchmark problems #2 and #3 not only validate our numerical solver but also provide insights needed for the optimization studies reported in Chapter titled," Optimization of CO2 Sequestration

**4.3. Simulation of benchmark problem #2 – enhanced CH4 recovery in combination with**

For decades, the oil and gas industry has been using a reservoir engineering technique to increase the oil/gas production from mature reservoirs, known as enhanced oil/gas recovery (EOR/EGR). As the original formation fluid (oil or natural gas) gets extracted, pressure in the reservoir gradually decreases. Such de-pressurization process makes it increasingly difficult to maintain the desired production rate. The reservoir needs to be re-pressurized to mitigate the drop of oil/gas production. One of the means to do this is to inject CO2 into the mature reservoir. With void space being occupied by the injected CO2, remaining oil/gas is pushed out of the reservoir. Meanwhile, the depleted reservoir becomes an ideal carbon sink for longterm storage. EOR/EGR with CO2 sequestration, also known as CSEOR/CSEGR, has been

In benchmark problem #2, a five-spot pattern domainis considered for modeling. The five-spot pattern is a common configuration for oil/gas production. A schematic of the reservoir is shown

Natural gas is produced at the four upper corners of the reservoir, while CO2 is injected in the middle of the domain at the bottom-most part. This is a direct result of CO2 being heavier than CH4 under the reservoir conditions. Injection of CO2 at the bottom avoids gas mixing and creates better sweep efficiency. The main goal of this benchmark simulation is to identify the gas recovery factor, defined as the ratio of enhanced CH4 production to the original remaining CH4 amount until the shutdown of the production well. Additionally, the time until production

As seen in Figure 8, the CO2 plume is nicely captured in the simulation.

Saline Aquifers".

316 CO2 Sequestration and Valorization

in Figure 3.

**CO2 sequestration in depleted gas reservoirs**

frequently used by the industry due to its strong economic merits.

shut-down, which is defined as the moment when the production contains up to 20% of CO2 by mass, needs to be determined.

Due to symmetry, only a quarter of the domain is modeled, as shown in Figure 10 as the volume bounded by the solid lines. The dimension of the modeled domain is 201.19 m ×201.19 m with thickness of 45.72 m. Due to the relatively strong diffusion, the discretization length has a strong influence on the gas mixing [10].It is therefore strictly specified as 4.572 m for both vertical and horizontal direction.Figure 10 shows the CFD model and its mesh in the quarter computational domain of the five-spot reservoir.

The hydrogeological properties assigned to the model are summarized in Table 8.

**Table 8.** Hydrogeological properties of the domain for benchmark problem #2

Initial conditions and boundary conditions and some other parameters of the domain are given in Table 9.

**Table 9.** Simulation parameters for benchmark problem #2

The termination of the simulation depends solely on the mass fraction of CO2 in the reservoir. It takes about 30 minutes of CPU time to run 2,000 days of simulation before major CO2 contamination occurs. The recovery factor, production shut-down time, pressure and CO2 saturation distribution in the reservoir are investigated and compared with the results of other investigators.

In Figure 11, our results using TOUGH2 are shown as the large graph, while results of simulations from other investigators are shown in the inner box. Table 10 provides compari‐ sons for recovery factor and production well shut-down time with other investigators' simulations [12].

**Table 10.** Comparisons of recovery factor and production shut-down time

**Figure 11.** History of enhanced CH4 recovery for benchmark problem #2

Permeability Horizontal: 50 mDarcy, Vertical: 5 mDarcy

Initial conditions and boundary conditions and some other parameters of the domain are given

Boundary Conditions No mass flow at all boundaries; Constant pressure at CH4 production

well

**Recovery Factor Production Well Shut-down Time**

Porosity 0.23

Relative Permeability Liquid: immobile, Gas: linear (*kr* = *S*)

Residual Brine Saturation 0

**Table 8.** Hydrogeological properties of the domain for benchmark problem #2

**Table 9.** Simulation parameters for benchmark problem #2

in Table 9.

318 CO2 Sequestration and Valorization

investigators.

simulations [12].

Capillary Pressure None

**Thermal Condition Isothermal** Initial Condition (Temperature) 66.7 oC Initial Condition (Pressure) 3.55 MPa

Initial CO2 Mass Fraction *XCO2* = 0 Initial CH4 Mass Fraction *Xsm* = 1

Injection Rate 0.1 kg/s until shut-down

The termination of the simulation depends solely on the mass fraction of CO2 in the reservoir. It takes about 30 minutes of CPU time to run 2,000 days of simulation before major CO2 contamination occurs. The recovery factor, production shut-down time, pressure and CO2 saturation distribution in the reservoir are investigated and compared with the results of other

In Figure 11, our results using TOUGH2 are shown as the large graph, while results of simulations from other investigators are shown in the inner box. Table 10 provides compari‐ sons for recovery factor and production well shut-down time with other investigators'

TOUGH2 (WUSTL) 61.4% 2063 days TOUGH2 (CO2/CRC) 58% 1987 days MUFTE (U. Stuttgart) 53% 1894 days IPARS (U. Texas) 55% 1891 days

**Table 10.** Comparisons of recovery factor and production shut-down time

To visualize how the displacement process of CH4 by CO2 works, the density and CO2 mass fraction profiles at production shut-down are examined in Figure 12.

**Figure 12.** (a) Density profile (b) CO2 mass fraction profile at production shut-down for benchmark problem #2

Figure 11, Table 10 and Figure 12 lead to the following four conclusions: 1) Small variations in the results among different simulations with different users are unavoidable. Such variations are expected because some parameters are intentionally left unspecified. 2) Our results are in satisfactory agreement with those of other investigators. 3) It can be seen that the injected CO2 migrates from the near lower corner to the far upper corner in a semi-spherical fashion. Unlike SAGCS, in situ CO2 tends to sink to the bottom of the reservoir. This indicates strong gravity segregation caused by the density difference. 4) Production gas contamination caused by upward movement of CO2 occurs at the production well despite the gravity segregation. This is due to the strong convective flow near the production well and mass diffusion.

#### **4.4. Simulation of benchmark problem #3 – CO2 injection in a heterogeneous geological formation**

Accurate estimation of in situ CO2 dissolution into the ambient brine is another important aspect of SAGCS simulation, since CO2 becomes securely sequestrated once dissolved. Overestimation of CO2 dissolution could lead to underestimating the potential leakage; on the other hand, underestimation of CO2 dissolution would result in inefficient utilization of the formation's storage potential. In the meantime, it is instructive to model a realistic geological reservoir with heterogeneous hydrogeological properties for more realistic estimation of CO2 storage capacity. In benchmark problem #3, part of the Johanson formation off the Norwegian coast is modeled for SAGCS [11]. The Johanson formation is a highly heterogene‐ ous formation, especially in its porosity and permeability, as shown in Figure 13. CO2 is injected in the middle of the modeled formation at 50 m from the bottom.  The goal of this benchmark study is to identify the amount of dissolved CO2, the amount of CO2 still in gaseous phase, and how these amounts evolve with respect to time. This study is very instructive to understand the dissolution process of injected CO2 under reservoir conditions.

**Figure 13.** Johanson formation's porosity heterogeneity for benchmark problem #3[11]

The dimension of the modeled portion of the Johanson formation is 9600 m × 8900 m with varying thickness between 50 m to 150 m. The injection well is located at (*x* = 5440 m, *y* = 3300 m) over the bottom 50 m of the formation. The coordinates of vertices of 54756 hexahedral cells in Cartesian system have been provided for geometry construction [11]. Figure 14 shows the final CFD model of the Johanson formation.

**Figure 14.** (a) Front view and (b) Rear view of the modeled Johanson formation

by upward movement of CO2 occurs at the production well despite the gravity segregation. This is due to the strong convective flow near the production well and mass diffusion.

**4.4. Simulation of benchmark problem #3 – CO2 injection in a heterogeneous geological**

Accurate estimation of in situ CO2 dissolution into the ambient brine is another important aspect of SAGCS simulation, since CO2 becomes securely sequestrated once dissolved. Overestimation of CO2 dissolution could lead to underestimating the potential leakage; on the other hand, underestimation of CO2 dissolution would result in inefficient utilization of the formation's storage potential. In the meantime, it is instructive to model a realistic geological reservoir with heterogeneous hydrogeological properties for more realistic estimation of CO2 storage capacity. In benchmark problem #3, part of the Johanson formation off the Norwegian coast is modeled for SAGCS [11]. The Johanson formation is a highly heterogene‐ ous formation, especially in its porosity and permeability, as shown in Figure 13. CO2 is injected in the middle of the modeled formation at 50 m from the bottom.  The goal of this benchmark study is to identify the amount of dissolved CO2, the amount of CO2 still in gaseous phase, and how these amounts evolve with respect to time. This study is very instructive to understand

the dissolution process of injected CO2 under reservoir conditions.

**Figure 13.** Johanson formation's porosity heterogeneity for benchmark problem #3[11]

final CFD model of the Johanson formation.

The dimension of the modeled portion of the Johanson formation is 9600 m × 8900 m with varying thickness between 50 m to 150 m. The injection well is located at (*x* = 5440 m, *y* = 3300 m) over the bottom 50 m of the formation. The coordinates of vertices of 54756 hexahedral cells in Cartesian system have been provided for geometry construction [11]. Figure 14 shows the

**formation**

320 CO2 Sequestration and Valorization

The porosity and permeability of the modeled Johanson formation areshown in Figure 15.

**Figure 15.** (a) Porosity and (b) Permeability of the modeled Johanson formation

Hydrogeological properties of modeled Johanson formation are summarized in Table 11.

**Table 11.** Hydrogeological properties of the modeled Johanson formation

Initial conditions, boundary conditions and other parameters of the modeled Johanson formation are summarized in Table 12.

**Table 12.** Simulation parameters for the modeled Johanson formation

Both gaseous and aqueous CO2 accumulations after 50 years are considered as benchmark criteria for making comparisons with the simulations of other investigators. In Figure 16 our results using TOUGH2 are shown as the large graph, while results from other simulations are shown in the inner box.

**Figure 16.** Gaseous and aqueous CO2 accumulations for 50 years

Table 13 provides additional quantitative comparisons [12].

Initial conditions, boundary conditions and other parameters of the modeled Johanson

**Thermal Condition Isothermal**

Initial CO2 Mass Fraction *XCO2* = 0 Initial Salt Mass Fraction *Xsm* = 0.1

**Table 12.** Simulation parameters for the modeled Johanson formation

**Figure 16.** Gaseous and aqueous CO2 accumulations for 50 years

Initial Temperature 0.03 oC/m; 100 oC at 3000 m depth Initial Pressure 1075 Pa/m, 30.86 MPa at 3000 m depth

Injection Rate 15 kg/s (for 25 years), 0 kg/s thereafter Discretization Number of computational grids: 18804

Both gaseous and aqueous CO2 accumulations after 50 years are considered as benchmark criteria for making comparisons with the simulations of other investigators. In Figure 16 our results using TOUGH2 are shown as the large graph, while results from other simulations are

Fixed-state on lateral boundaries No mass flow on fault, top and bottom boundaries

Non-uniform for *x*-, *y*-, *z*-directions

formation are summarized in Table 12.

322 CO2 Sequestration and Valorization

Boundary Conditions

shown in the inner box.

**Table 13.** Comparisons of gaseous and aqueous CO2 accumulations after 50 years

A comparison of the CO2 migration after 50 years is given in Figure 17.

**Figure 17.** CO2 saturation in the formation after 50 years, plan view

Figure 16, Table 13, and Figure 17 lead to four conclusions: 1) Small variations in the results among different simulations with different users are unavoidable. Such variations are expected because some parameters are intentionally left unspecified. 2) Our results are in satisfactory agreement with the results of other investigators. 3) CO2 dissolution into the ambient pore‐ water is a process that takes place very slowly. 4) The greater slope of aqueous CO2 during the first 25 years (when injection continues) implies enhanced carbon dissolution due to convec‐ tion during the first 25 years.

#### **4.5. Conclusions of benchmark problems**

The benchmark simulations presented in this chapter demonstrate that the TOUGH2 numer‐ ical simulator is capable of producing accurate and consistent results for various types of problems related to GCS. These simulations allow us to conduct simulations of large-scale SAGCS in identified saline formations with confidence, and proceed towards the development of a numerical optimization module for TOUGH2 and perform optimization designs of innovative reservoir engineering techniques for enhanced SAGCS safety and storage efficien‐ cy, as reported in Chapter titled," Optimization of CO2 Sequestration Saline Aquifers".

## **5. Simulations of GCS in identified large scale saline aquifers**

Accurate large-scale simulations of existing (completed/continuing) SAGCS projects for identified known aquifers are crucial to create confidence in the future deployment of SAGCS projects. Although detailed history-matching simulations of existing SAGCS projects are challenging due to various uncertainties, e.g., in the reservoir topography and hydrogeology, the simulations can still provide informative insights into several aspects of SAGCS, such as the variance in multiphase flow properties, integrity of the geological seals, and the mechanism of CO2 trapping. Such insights are essential for better understanding of the nature of SAGCS and its best practices for deployment. Detailed history-matching simulations have always been an important part in the SAGCS research activity. In the following sections, the results of SAGCS simulations for three large-scale identified aquifers are described.

#### **5.1. SAGCS simulation for Mt.Simon formation**

Located in the Illinois basin, the Mt.Simon sandstone formationis a huge saline aquifer that covers most of Illinois, southwestern Indiana, southern Ohio and western Kentucky. The estimated storage capacity of the Mt.Simon formation ranges from 27000 to 109000 million tons of CO2[16],[17]. The Midwest Geological Sequestration Consortium (MGSC) is the regional consortium conducting studies of the possibility of large-scale GCS throughout the Illinois basin. ADM-Decatur GCS Project and FutureGen 2.0 Project are the two best-known SAGCS projects being currently carried out in theMt.Simon formation.

The depth of the Mt.Simon formation varies significantly throughout its coverage [18],[19]. In the southern part, it reaches as deep as 4300 m below mean sea level (MSL) while it increases to 80 m below MSL in the north. Consequently, a south-north geological slope of approximately 8 m/km has been estimated. The thickness of the Mt.Simon formation also changes signifi‐ cantly.A maximum thickness of 800 m in the north has been measured while it diminishes to zero further to the south. Other than the variance in topography, analysis of rock samples has suggested strong anisotropy in the formation's hydrogeological properties, with porosity ranging from 0.062 to 0.2 and permeability ranging from 5 mDarcyto 1000 mDarcy. Lowpermeable Eau Claire shale, which sits above the Mt.Simon formation, serves as the caprock. Except for some small regions near the Mississippi River, Eau Claire shale is very thick (more than 90 m) throughout most of the Illinois basin. The security of SAGCS over Mt.Simon formation is therefore assured by the continuous coverage of Eau Claire shale. A Precambrian granite formation stretches beneath the Mt.Simon saline aquifer.

A recent geological survey has suggested an area in the center of the Mt.Simon formation to be the core injection area – an area in which future storage sites are likely to be located. This core injection area is indicated as the area compressed by the white boundary in Figure 12, along with the elevation information of the Mt.Simon formation. As can be seen from Figure 12, both the ADM and FutureGen 2.0 projects are located in the core injection area.

Mt. Simon sandstone is a typical stratified saline formation. According to the geological survey, strong anisotropy in porosity, permeability and capillary pressure exists through the entire

**Figure 18.** Core injection area and elevation of Mt.Simon sandstone

**5. Simulations of GCS in identified large scale saline aquifers**

SAGCS simulations for three large-scale identified aquifers are described.

**5.1. SAGCS simulation for Mt.Simon formation**

324 CO2 Sequestration and Valorization

projects being currently carried out in theMt.Simon formation.

granite formation stretches beneath the Mt.Simon saline aquifer.

Accurate large-scale simulations of existing (completed/continuing) SAGCS projects for identified known aquifers are crucial to create confidence in the future deployment of SAGCS projects. Although detailed history-matching simulations of existing SAGCS projects are challenging due to various uncertainties, e.g., in the reservoir topography and hydrogeology, the simulations can still provide informative insights into several aspects of SAGCS, such as the variance in multiphase flow properties, integrity of the geological seals, and the mechanism of CO2 trapping. Such insights are essential for better understanding of the nature of SAGCS and its best practices for deployment. Detailed history-matching simulations have always been an important part in the SAGCS research activity. In the following sections, the results of

Located in the Illinois basin, the Mt.Simon sandstone formationis a huge saline aquifer that covers most of Illinois, southwestern Indiana, southern Ohio and western Kentucky. The estimated storage capacity of the Mt.Simon formation ranges from 27000 to 109000 million tons of CO2[16],[17]. The Midwest Geological Sequestration Consortium (MGSC) is the regional consortium conducting studies of the possibility of large-scale GCS throughout the Illinois basin. ADM-Decatur GCS Project and FutureGen 2.0 Project are the two best-known SAGCS

The depth of the Mt.Simon formation varies significantly throughout its coverage [18],[19]. In the southern part, it reaches as deep as 4300 m below mean sea level (MSL) while it increases to 80 m below MSL in the north. Consequently, a south-north geological slope of approximately 8 m/km has been estimated. The thickness of the Mt.Simon formation also changes signifi‐ cantly.A maximum thickness of 800 m in the north has been measured while it diminishes to zero further to the south. Other than the variance in topography, analysis of rock samples has suggested strong anisotropy in the formation's hydrogeological properties, with porosity ranging from 0.062 to 0.2 and permeability ranging from 5 mDarcyto 1000 mDarcy. Lowpermeable Eau Claire shale, which sits above the Mt.Simon formation, serves as the caprock. Except for some small regions near the Mississippi River, Eau Claire shale is very thick (more than 90 m) throughout most of the Illinois basin. The security of SAGCS over Mt.Simon formation is therefore assured by the continuous coverage of Eau Claire shale. A Precambrian

A recent geological survey has suggested an area in the center of the Mt.Simon formation to be the core injection area – an area in which future storage sites are likely to be located. This core injection area is indicated as the area compressed by the white boundary in Figure 12, along with the elevation information of the Mt.Simon formation. As can be seen from Figure

Mt. Simon sandstone is a typical stratified saline formation. According to the geological survey, strong anisotropy in porosity, permeability and capillary pressure exists through the entire

12, both the ADM and FutureGen 2.0 projects are located in the core injection area.

depth of the formation. Based on variance of porosity, the Mt. Simon formation can be distinguished as four subunits, namely an upper unit with sandstone and shale tidally influenced and deposited, a middle unit with relatively clean sandstone, an Arkosic unit with highly porous and permeable sandstone, and a lower unit with decreased porosity and permeability. The high porosity and permeability of the Arkosic unit makes it an ideal candidate for the injection site. When modeling, these four subunits of Mt.Simon are further divided into 24 layers, each of which has a layer-averaged porosity and permeability value [16],[19]. With the consideration of data availability, a candidate site for future sequestration project, the Weaver-Horn #1 well (WH #1 well is shown as the red dot in Figure 18), has been chosen for our simulation study. The detailed well log of WH #1 is summarized in Table 14[20]. It is desired to model the anisotropy of hydrogeological properties as accurately as possible to capture its effect on in situ CO2 transport. It should be noted that the lower unit of the Mt.Simon formation is not considered in the model due to its absence near WH #1 well. Both Eau Claire shale and Precambrian granite are modeled as impermeable formations.

A cylindrical model of the Mt.Simon formation is constructed. For thermal condition,the model uses calculated values with a thermal gradient of 9.2°C/km. The reservoir pressure is assumed to be hydrostatic pressure with a gradient of about 10.8 MPa/km from the ground surface. Salinity is assumed to increase with the depth, starting from 235 mg/L at 450 m below ground surface with a gradient of 12.8 mg/L per meter in depth. A north-south geological slope of 8 m/km is also considered in the model. A "no-flux" boundary condition is applied at the top and bottom of the model, representing the impermeable upper and lower bounding forma‐ tions. A "fixed-state" boundary condition is imposed at the lateral boundary to represent an essentially "open" reservoir. The permeability and porosity of the 24 sublayers can be seen in Figure 19.

**Table 14.** Porosity, permeability and characteristic capillary pressure of the 24 layers of Mt.Simon at injection site WH #1 [20]

**Figure 19.** (a) Permeability, (b) porosity, and computational mesh of the 24 sublayers of the Mt.Simon formation mod‐ el at WH #1 well

Due to the relatively high porosity and permeability, CO2 injection is assigned at the bottom Arkosic unit (bottom three sublayers). The injection rate is assigned to be 5 million tons per year and injection lasts for 50 years. CO2 footprint at 5th year, 25th year and 50th year since the beginning of injection is examined.

**Figure 20.** Saturation of gaseous CO2 at (a) 5th (b) 25th and (c) 50th year of injection

As seen in Figure 20, the CO2 plume evolves with a complex spatial pattern during the 50 years of injection. Within the Arkosic unit where the injector is located, extensive lateral migration with relatively higher concentration of gaseous CO2 is observed. In the overlying sublayers, however, a strong secondary sealing effect that retards the vertical migration of gaseous CO2 is observed, as can be seen from the pyramid-shaped sub-plume. Detailed analysis of secon‐ dary sealing effect is done as follows. As shown in Figure 20, the injected CO2 migrates laterally away from the injector within the highly permeable Arkosic unit in the first 5 years after the start of injection. Simultaneously, buoyancy also leads to upward movement of CO2 until it encounters the immediately overlying low-permeability sublayer (sublayer #20). The low permeability of sublayer #20 directly results in higher capillary pressure experienced by mobile CO2, and thus a stronger vertical pressure gradient is required for mobile CO2 to penetrate sublayer #20. When the capillary pressure is greater than the phase pressure of CO2, sublayer #20 appears to be "impermeable" to the underlying CO2 plume. Consequently, gaseous CO2 accumulates under this layer and continues spreading out laterally, finally reaching a maxi‐ mum extent of approximately 3000 m. Meanwhile, the increased CO2 column under sublayer #20 brings up its phase pressure. Once the phase pressure of CO2 exceeds the entry pressure of sublayer #20, mobile CO2 breaks the capillary barrier of its overlying layer and starts to penetrate it. Such accumulation-penetration-breakthrough behavior of gaseous CO2 occurs each time the upward migrating CO2 encounters an overlying sublayer with lower permea‐ bility. Because the high capillary entry pressure of the overlying layer temporarily prevents CO2 from migrating upwards,this phenomenon is identified as the secondary sealing effect. As can be seen from Figure 20, this effect is a very effective means to retard the upward migration of in situ CO2. Its contribution makes gaseous CO2 barely reach the Eau Clair shale even after 50 years of injection.

#### **5.2. SAGCS simulation of Frio formation**

**Figure 19.** (a) Permeability, (b) porosity, and computational mesh of the 24 sublayers of the Mt.Simon formation mod‐

Due to the relatively high porosity and permeability, CO2 injection is assigned at the bottom Arkosic unit (bottom three sublayers). The injection rate is assigned to be 5 million tons per

el at WH #1 well

**Sub-Unit Sublayer Layer Depth (m)**

326 CO2 Sequestration and Valorization

Upper Unit

Middle Unit

Arkosic Unit

#1 [20]

**Mean Porosity**

 2140 – 2150 0.061 5 0.37 2150 – 2182 0.109 300 0.06 2182 – 2197 0.074 10 0.28 2197 – 2203 0.083 3.6 0.4875 2203 – 2230 0.195 110 0.1 2230 – 2232 0.071 1.1 0.8 2232 – 2280 0.13 210 0.083

 2280 – 2322 0.083 5.4 0.4125 2322 – 2331 0.24 150 0.0875 2331 – 2340 0.088 8 0.35 2340 – 2350 0.156 800 0.095 2350 – 2370 0.25 80 0.125 2370 – 2378 0.163 900 0.095 2378 – 2385 0.195 105 0.1007 2385 – 2399 0.163 800 0.05 2399 – 2406 0.136 72 0.1167 2406 – 2412 0.156 700 0.05 2412 – 2424 0.129 160 0.09 2424 – 2430 0.161 850 0.05 2430 – 2462 0.128 60 0.15

21 2462 – 2500 0.202 1000 0.05 22 2500 – 2502 0.14 190 0.09 23 2502 – 2537 0.151 1000 0.04

**Table 14.** Porosity, permeability and characteristic capillary pressure of the 24 layers of Mt.Simon at injection site WH

**Mean Permeability (mDarcy)**

**Characteristic Capillary Pressure (bar)**

> The SAGCS pilot project for the Frio deep saline formation near the Mexican Gulf Coast is the subject of study in this simulation. The Frio project has two characteristics that make it

attractive for numerical study. First, it is a completed pilot project with detailed field data available; secondly, hysteresis information of relative permeability and capillary pressure has been obtained from the core sample of the Frio saline formation. The hysteresis effect is an important factor in obtaining accurate estimation of CO2 migration and dissolution for fullterm SAGCS simulation. A full-term simulation refers to a simulation that investigates the fate of in situ CO2 through the entire life cycle of a SAGCS project, which consists of both injection and post-injection periods.

The Frio SAGCS pilot project was conducted at the South Liberty oil field operated by Texas American Resources in Dayton, Texas. Starting from October 4, 2004, 1600 tons of CO2 was injected into the Frio formation about 1500 m below the ground surface within 10 days. The Frio formation consists of brine-bearing sandstone with high permeability beneath the Gulf Coast. It is a relatively thin sandstone layer only 23 m thick. A steep geological slope of 16° from south to north has been identified for Frio formation [21]. The Frio pilot project employed one injection well and one observation well about 33 m to its north. Other than the conventional pre-injection geological surveys, laboratory analysis of core samples suggested the hysteresis behavior of relative permeability and capillary pressure in the Frio formation. The hysteresis has been considered in our simulations.

The hydrogeological parameters of the modeled Frio formation are summarized in Table 15.

**Table 15.** Geometry and hydrogeological parameters for Frio formation

Characteristics of capillary pressure and relative permeability have been obtained from mercury-injection laboratory experiments on core samples from Frio formation, given in Figure 21. Hysteresis in both capillary pressure and relative permeability can be clearly observed. Drainage (of pore-water) curves are marked red and imbibition (of pore-water) curves are marked blue. When multiple drainage-imbibition cycles occur, different imbibition curves represent different orders of drainage-imbibition cycles. The primary imbibition curve, i.e., when brine imbibition occurs for the first time, is depicted as a bold solid curve. Since only one drainage-imbibition cycle takes place when continuous CO2 injection is imposed before it is permanently shut down, only the primary imbibition curve needs to be considered in the model.

attractive for numerical study. First, it is a completed pilot project with detailed field data available; secondly, hysteresis information of relative permeability and capillary pressure has been obtained from the core sample of the Frio saline formation. The hysteresis effect is an important factor in obtaining accurate estimation of CO2 migration and dissolution for fullterm SAGCS simulation. A full-term simulation refers to a simulation that investigates the fate of in situ CO2 through the entire life cycle of a SAGCS project, which consists of both injection

The Frio SAGCS pilot project was conducted at the South Liberty oil field operated by Texas American Resources in Dayton, Texas. Starting from October 4, 2004, 1600 tons of CO2 was injected into the Frio formation about 1500 m below the ground surface within 10 days. The Frio formation consists of brine-bearing sandstone with high permeability beneath the Gulf Coast. It is a relatively thin sandstone layer only 23 m thick. A steep geological slope of 16° from south to north has been identified for Frio formation [21]. The Frio pilot project employed one injection well and one observation well about 33 m to its north. Other than the conventional pre-injection geological surveys, laboratory analysis of core samples suggested the hysteresis behavior of relative permeability and capillary pressure in the Frio formation. The hysteresis

The hydrogeological parameters of the modeled Frio formation are summarized in Table 15.

Permeability Isotropic, 1 mDarcy

Relative Permeability van Genuchten-Mualem with hysteresis Capillary Pressure van Genuchten-Mualem with hysteresis

Initial Conditions *P* = 15.2 MPa, *T* = 59oC

Characteristics of capillary pressure and relative permeability have been obtained from mercury-injection laboratory experiments on core samples from Frio formation, given in Figure 21. Hysteresis in both capillary pressure and relative permeability can be clearly observed. Drainage (of pore-water) curves are marked red and imbibition (of pore-water) curves are marked blue. When multiple drainage-imbibition cycles occur, different imbibition curves represent different orders of drainage-imbibition cycles. The primary imbibition curve, i.e., when brine imbibition occurs for the first time, is depicted as a bold solid curve. Since only

Initial CO2 Mass Fraction *XCO2* = 0 Initial Salt Mass Fraction *Xsm* = 0.093

**Table 15.** Geometry and hydrogeological parameters for Frio formation

Porosity 0.28

Residual Brine Saturation 0.15 Residual CO2 Saturation 0.2

and post-injection periods.

328 CO2 Sequestration and Valorization

has been considered in our simulations.

**Figure 21.** Capillary pressure and relative permeability characteristics of Frio formation[22],[23]

Following Doughty et al.'s work [22], a rectangular portion of Frio formation with dimension of 2500 m northwest-southwest, 800 m northeast-southeast, and 23 m thickness is modeled, as shown in Figure 22. The injection well is located at a point with coordinate (*x*=560 m, *y*=800 m) from the lower left corner of the computational domain. Although the formation is 23 m thick, injection takes place only over the first 8 m from the caprock. An observation well is located 33 m to the north. Because flow transport is most intense near the injection and observation wells, it is evolved in a computational domain of dimension 30 m × 30 m whose mesh is refined for accurate capture of the flowpattern. The injection and observation well locations, well depth, computational mesh, and north-south slope of the numerical model are all shown in Figure 22.

The simulation time is set at 10 days to match the actual duration of injection. It takes approx‐ imately 12 hours of CPU time for the simulation to finish. The profiles of gaseous phase CO2 at the end of injection in the vertical cross-section containing both injection and observation wells are shown in Figure 23. Doughty et al.'s result [22] is also shown in Figure 23 for comparison.

Additionally, Figure 24 shows the CO2 saturation profiles at the injection and observation wells obtained by our simulation. They are compared with those given by Doughty et al. and the reservoir saturation tool (RST) logs [22]. The RST well logs are actual measurements in the field during the pilot project.

**Figure 22.** Model geometry and mesh in a portion of Frio formation and zoomed-in side view of the injection and observation wells

**Figure 23.** CO2 footprint on 10th day when injection stops (comparison with Doughty et al.'s work[22])

As seen from Figure 23, the highly asymmetric CO2 plume suggests a strong tendency of movement upwards in the direction of the geological slope. Unlike the case of Mt. Simon SAGCS, the CO2 plume in the Frio project is shaped like an inverted pyramid, which implies the lack of secondary sealing effect. Both the asymmetric migration and inverted pyramidshaped plume indicate strong evidence of the dominant role of gravity segregation in deter‐ mining the in situ CO2 migration. Considering the relatively short-term injection (10 days) for the Frio SAGCS project, this implies that in situ CO2 migrates mostly convectively. Further‐ more, it demonstrates that the poor permeability caprock layer above the injection serves quite well as the CO2 barrier.

Numerical Simulation of CO2 Sequestration in Large Saline Aquifers http://dx.doi.org/10.5772/57065 331

**Figure 24.** CO2 saturation profiles given by the numerical simulations and the RST logs

Comparing our results with those from Doughty et al. [22], the following conclusions can be drawn: (1) Overall our results are in good agreement with those of Doughty et al. for the plume shape, tendency of plume migration induced by the slope, distance of plume migration, and gaseous CO2 saturation, etc. (2) Discrepancy still exists at the detailed simulation level. The results show that in our simulations, CO2 saturation at the injection well reaches a maximum of 0.8 by the 10th day of injection. Although being consistent with Doughty et al.'s work, it differs from the field data. Results from the RST measurement suggest a CO2 saturation value of 1.0, i.e., dry-out of brine, occurs adjacent to the injection well. The occurrence of brine dryout is fairly common near the injection well due to the strong pressure gradient. However, the absence of brine dry-out in both our and Doughty et al.'s simulations can be explained by the designated brine residual saturation value. In our TOUGH2 simulations, a brine residual saturation value of 0.15 is pre-assigned to the entire computational grid including the grid near the injection wells. Since residual saturation describes the minimum saturation value of a certain phase being displaced, it means that at least 15% of the pore space will remain occupied by brine regardless of the pressure gradient. A direct result is the capped CO2 saturation value at 0.85 and the absence of brine dry-out. (3) Our simulation shows quicker decrease in gas saturation during the injection interval. In Doughty et al.'s work, the gas saturation only drops from 0.8 to 0.65 for the upper 5 m of injection depth. In contrast, it drops from 0.8 to 0.4 in our simulation. This implies stronger buoyancy in our simulation, and thus results in a steeper inclined CO2-brine interface. This also explains the slight overshoot in the plume migration to the north in our simulation.

#### **5.3. SAGCS simulation for Utsira formation**

As seen from Figure 23, the highly asymmetric CO2 plume suggests a strong tendency of movement upwards in the direction of the geological slope. Unlike the case of Mt. Simon SAGCS, the CO2 plume in the Frio project is shaped like an inverted pyramid, which implies the lack of secondary sealing effect. Both the asymmetric migration and inverted pyramidshaped plume indicate strong evidence of the dominant role of gravity segregation in deter‐ mining the in situ CO2 migration. Considering the relatively short-term injection (10 days) for the Frio SAGCS project, this implies that in situ CO2 migrates mostly convectively. Further‐ more, it demonstrates that the poor permeability caprock layer above the injection serves quite

**Figure 23.** CO2 footprint on 10th day when injection stops (comparison with Doughty et al.'s work[22])

**Figure 22.** Model geometry and mesh in a portion of Frio formation and zoomed-in side view of the injection and

well as the CO2 barrier.

observation wells

330 CO2 Sequestration and Valorization

The Sleipner project near the Norwegian coast on the North Sea is probably the most presti‐ gious, important and successful SAGCS demonstration so far. It has the most complete topographic description, industrial-scale injection amount, and long-term monitoring data. Nevertheless, great uncertainties still exist for accurate reservoir-scale simulation of the Sleipner SAGCS project. Simulation studies of this project can provide helpful insights in understanding the transport behavior of in situ CO2 and the reservoir performance.

Starting from 1996, the Sleipner field in the North Sea has been the host of the world's first commercial SAGCS project. CO2 is captured from the gas mixture produced from a nearby deeper natural gas reservoir. To date, approximately 1 million tons of supercritical CO2 has been sequestered annually. Utsira saline aquifer is the target formation for permanent carbon sequestration for the Sleipner SAGCS project. The Utsira formation is located at a depth of 800 m – 1100 m from the seabed with thickness of about 200 m – 250 m. The injection site is located at the southern portion of the Utsira formation. A 250 m – 330 m thick shale layer known as the Nordland Formation serves as the caprock, and core testing has suggested its potential for bearing a CO2 column of at least 100 m but perhaps up to 400 m (depending on the in situ conditions). It is estimated that the Utsira formation has permeability of about 1– 8 mDarcy, porosity of about 0.35 – 0.4, and temperature of about 34°C – 37°C. It is also estimated that the reservoir bears hydrostatic pressure from its overburden formations. Similar to the Mt.Simon formation, the Utsira formation is also highly stratified, consisting of sublayers with highpermeability sandstone and low-permeability shale. Therefore, it is expected that the secon‐ dary sealing effect will occur. Figure 25 shows a 2-D seismic image taken in 2008 revealing the CO2 plume in the Utsira formation. Multiple layers can be distinctly identified from the seismic image.

Two numerical models have been constructed to study the Sleipner SAGCS project. The first model is a generalized axisymmetric layered model for estimating the ballpark migration of in situ CO2. The purpose of this simulation is to determine the secondary sealing effect and gain an overview of the plume migration within the Utsira formation. The second model describers a 48 km2 area of detailed topmost sandstone layer (marked as Layer #9 in Figure 26). Layer #9 is of particular interest regarding the safety of the sequestration project, as it is the layer within which the highest concentration of gaseous CO2 exists and the most significant plume migration occurs. Detailed topography of Layer #9 is shown in model #2, making it a complicated 3D model. The 3D Layer #9 model is introduced to investigate the effect of actual topography on in situ CO2 migration, while avoiding intensive computational effort associated with full 3D modeling and simulation of the entire Utsira formation.

#### *5.3.1. Model #1 – Generalized stratified model of the Utsira formation*

A pre-injection geological survey unveiled the layered structure of the Utsira formation. The majority of the formation can be identified as an eight-layered structure, but one extra layer needs to be added to the structure near the injection site due to the existence of a sand wedge. Therefore, a cylindrical domain with nine alternating sandstone and shale layers is constructed as shown in Figure 26. According to the seismic survey, it is assumed that all four shale layers have identical thickness of 5 m, four shallower sandstone layers have identical thickness of 25 m, and the bottom sandstone layer has a thickness of 60 m. This adds up to a total 180 m thickness for the modeled Utsira formation. The lateral radius of the generalized cylindrical model reaches 100 km, which is about the same as the actual extent of the southern part of the Utsira formation. According to Audigane et al. [25], all sandstone layers have identical and

**Figure 25.** eismic image of Utsira formation after 9 years of injection, S-N cross-section [24]

Sleipner SAGCS project. Simulation studies of this project can provide helpful insights in

Starting from 1996, the Sleipner field in the North Sea has been the host of the world's first commercial SAGCS project. CO2 is captured from the gas mixture produced from a nearby deeper natural gas reservoir. To date, approximately 1 million tons of supercritical CO2 has been sequestered annually. Utsira saline aquifer is the target formation for permanent carbon sequestration for the Sleipner SAGCS project. The Utsira formation is located at a depth of 800 m – 1100 m from the seabed with thickness of about 200 m – 250 m. The injection site is located at the southern portion of the Utsira formation. A 250 m – 330 m thick shale layer known as the Nordland Formation serves as the caprock, and core testing has suggested its potential for bearing a CO2 column of at least 100 m but perhaps up to 400 m (depending on the in situ conditions). It is estimated that the Utsira formation has permeability of about 1– 8 mDarcy, porosity of about 0.35 – 0.4, and temperature of about 34°C – 37°C. It is also estimated that the reservoir bears hydrostatic pressure from its overburden formations. Similar to the Mt.Simon formation, the Utsira formation is also highly stratified, consisting of sublayers with highpermeability sandstone and low-permeability shale. Therefore, it is expected that the secon‐ dary sealing effect will occur. Figure 25 shows a 2-D seismic image taken in 2008 revealing the CO2 plume in the Utsira formation. Multiple layers can be distinctly identified from the seismic

Two numerical models have been constructed to study the Sleipner SAGCS project. The first model is a generalized axisymmetric layered model for estimating the ballpark migration of in situ CO2. The purpose of this simulation is to determine the secondary sealing effect and gain an overview of the plume migration within the Utsira formation. The second model

26). Layer #9 is of particular interest regarding the safety of the sequestration project, as it is the layer within which the highest concentration of gaseous CO2 exists and the most significant plume migration occurs. Detailed topography of Layer #9 is shown in model #2, making it a complicated 3D model. The 3D Layer #9 model is introduced to investigate the effect of actual topography on in situ CO2 migration, while avoiding intensive computational effort associated

A pre-injection geological survey unveiled the layered structure of the Utsira formation. The majority of the formation can be identified as an eight-layered structure, but one extra layer needs to be added to the structure near the injection site due to the existence of a sand wedge. Therefore, a cylindrical domain with nine alternating sandstone and shale layers is constructed as shown in Figure 26. According to the seismic survey, it is assumed that all four shale layers have identical thickness of 5 m, four shallower sandstone layers have identical thickness of 25 m, and the bottom sandstone layer has a thickness of 60 m. This adds up to a total 180 m thickness for the modeled Utsira formation. The lateral radius of the generalized cylindrical model reaches 100 km, which is about the same as the actual extent of the southern part of the Utsira formation. According to Audigane et al. [25], all sandstone layers have identical and

with full 3D modeling and simulation of the entire Utsira formation.

*5.3.1. Model #1 – Generalized stratified model of the Utsira formation*

area of detailed topmost sandstone layer (marked as Layer #9 in Figure

understanding the transport behavior of in situ CO2 and the reservoir performance.

image.

describers a 48 km2

332 CO2 Sequestration and Valorization

isotropic hydrogeological properties, as do the shale layers. Figure 26 shows the layered structure and computational mesh of the modeled Utsira formation as well as the location of CO2 injection. The permeability and porosity are 3 mDarcy and 0.42 for the sand layer and 10 mDarcy and 0.1025 for the shale layer. A temperature of 37o C and pressure of 11 MPa are applied as the pre-injection conditions. van Genuchten-Mualem functions are used to describe both relative permeability and capillary pressure. CO2 injection of 30 kg/s is assigned as a point source at the middle of the bottom-most sand layer.

**Figure 26.** Computational mesh and layered structure of the generalized 9-layered model of the Utsira formation

The simulation time is set at 15 years and the CO2 plume profile is examined for each year. Figure 27 shows the cross-sectional view of gaseous CO2 in the reservoir for ten consecutive years since the start of injection.

The results shown in Figure 27 provide evidence of strong secondary sealing effect for migration of in situ CO2. Similar to the case of Mt. Simon SAGCS, the injected CO2 first migrates upwards driven by buoyancy until it reaches the first shale layer. Due to the low permeability and high capillary entry pressure, CO2 is confined by this shale layer and is forced to migrate radially.  The accumulation-penetration-breakthrough takes place each time the CO2 plume encounters a new shale layer and forms an inverted pyramid shaped sub-plume, as documented clearly by the first and second year plume shapes in Figure 27. Due to the secondary sealing effect, in situ CO2 has very limited contact with the caprock of the Utsira formation by the third year of injection.

**Figure 27.** In situ CO2 distribution for 15 years of injection

Additionally, ten-year CO2 flux analysis has been made for the topmost sandstone layer (Layer #9) since it is critical to identify the accumulation of CO2 underneath the caprock. As shown in Figure 30, excellent agreement between our simulation and the seismic amplitudes analysis [26] is observed, suggesting the overall accuracy of our model despite some discrepancy at the detailed level. The flux analysis shown in Figure 28 also implies that the accumulation rate of CO2 in the topmost sandstone layer tends to increase until it becomes stabilized. This fact can be explained by mechanism of secondary sealing effect.

#### *5.3.2. Model #2 – Detailed 3D model of the Utsira Layer #9 formation*

In situ CO2 possesses strong potential to migrate upward due to buoyancy, and thus accumulates under the caprock unless capillary barrier is compromised. Previous experi‐ ence has demonstrated that the accumulation of CO2 under the caprock occurs in a relatively short period compared to the entire time span of SAGCS projects, and it is a major concern for storage security. Therefore, it is critical for a SAGCS project to identify the accumula‐ tion of CO2 and its tendency for migration underneath the caprock. With such informa‐ tion available, precautionary treatments can be deployed to avoid potential leakage. The Utsira formation near the injection site has been identified as a nine-layer structure as shown in Figure 29. The topmost sandstone layer, Layer #9, is of most interest since it has the

**Figure 28.** Gaseous CO2 accumulation in the topmost sandstone layer

and high capillary entry pressure, CO2 is confined by this shale layer and is forced to migrate radially.  The accumulation-penetration-breakthrough takes place each time the CO2 plume encounters a new shale layer and forms an inverted pyramid shaped sub-plume, as documented clearly by the first and second year plume shapes in Figure 27. Due to the secondary sealing effect, in situ CO2 has very limited contact with the caprock of the Utsira formation by the third year of

Additionally, ten-year CO2 flux analysis has been made for the topmost sandstone layer (Layer #9) since it is critical to identify the accumulation of CO2 underneath the caprock. As shown in Figure 30, excellent agreement between our simulation and the seismic amplitudes analysis [26] is observed, suggesting the overall accuracy of our model despite some discrepancy at the detailed level. The flux analysis shown in Figure 28 also implies that the accumulation rate of CO2 in the topmost sandstone layer tends to increase until it becomes stabilized. This fact can

In situ CO2 possesses strong potential to migrate upward due to buoyancy, and thus accumulates under the caprock unless capillary barrier is compromised. Previous experi‐ ence has demonstrated that the accumulation of CO2 under the caprock occurs in a relatively short period compared to the entire time span of SAGCS projects, and it is a major concern for storage security. Therefore, it is critical for a SAGCS project to identify the accumula‐ tion of CO2 and its tendency for migration underneath the caprock. With such informa‐ tion available, precautionary treatments can be deployed to avoid potential leakage. The Utsira formation near the injection site has been identified as a nine-layer structure as shown in Figure 29. The topmost sandstone layer, Layer #9, is of most interest since it has the

injection.

334 CO2 Sequestration and Valorization

**Figure 27.** In situ CO2 distribution for 15 years of injection

be explained by mechanism of secondary sealing effect.

*5.3.2. Model #2 – Detailed 3D model of the Utsira Layer #9 formation*

highest concentration of gaseous CO2 and has direct contact with the overlying caprock formation. Seismic surveys have shown striking growth of CO2 accumulation in Layer #9 between 1999 and 2006 as shown in Figure 29 [26].

The black dot in Figure 29 marks the location of the injection well, which is roughly 200 m under Layer #9. Two distinct local CO2 accumulations appeared after about three years of injection (recall that injection began in 1996), indicating CO2 began to accumulate under the caprock. However, CO2 migration in Layer #9 was not symmetric due to the topography of the caprock. The northward migration of initially impacted CO2, seen as the "body" of the plume in Figure 29, implies a local topographic dome; a prominent north-tending migration, seen as the "finger" of the plume in Figure 30. It implies the spill of local structurally trapped CO2 along a north-tending topographic ridge. CO2 migration along the north-tending ridge was rather fast at about 1 m/day between 2001 and 2004 [24].

In order to examine the plume's evolution within the topmost layer more closely, a 3D model of Utsira Layer #9 is created with detailed topography. It should be noted that only Layer #9, not the entire depth, is modeled because of the following considerations. To ensure the accurate capture of topographic effect on plume shaping, a computational domain with considerable fine mesh resolution has to be modeled based on geological survey data. The computational effort and thus the feasibility of highly detailed model of the entire Utsira formation is very intensive and time consuming. Secondly, CO2 has to breakthrough several layers of relatively low-permeability shale prior to reaching the topmost layer. While it is difficult to quantify the breakthrough of gaseous CO2, the quantification of CO2 feeding into the topmost layer (Layer

**Figure 29.** Amplitude maps of Layer #9 from 1999 to 2008, from Singh et al. [26]

#9) is rather reliable. Therefore, a model of only the topmost layer (Layer #9) could provide an ideal platform to investigate the effect of various parameters such as topography on the shaping of the CO2 plume, and also could be used for optimization purpose to achieve high efficiency sequestration while maintaining an affordable computational effort and cost.

A reservoir model with dimension of 1600 m × 4900 m with varying thickness was constructed. It covers the portion of the Utsira formation where the plume, shown in Figure 31, is located. As mentioned earlier, the topography of this portion of the Utsira formation is accurately modeled based on seismic geological survey data (provided by Zhu and Lu of the University of Indiana [27]) with 50 m × 50 m mesh resolution. Because only Layer #9 is modeled, the thickness of the computational domain varies from 3.5 m to 26.3 m with an average thickness of 11.3 m. However, to accurately capture the accumulation and upward and lateral movement of CO2, 37 layers are used along the thickness. The topmost layer and the bottom two layers are designated to represent the low permeability shale, while the 34 layers in the middle are assigned the properties of mudstone. In the 3D Layer #9 model, permeability anisotropy is considered with west-east permeability of 2 mDarcy, north-south permeability of 10 mDarcy, and vertical permeability of 200 mDarcy. A 3D overview of the Layer #9 model is shown in Figure 24.

**Figure 30.** 3D overview and plan view of the 3D Layer #9 model of Utsira indicating feeder locations (black dot: main feeder; cyan square: secondary feeder)

Table 16 summarizes the hydrogeological properties of the Layer #9 model.

**Table 16.** Hydrogeological properties of the Utsira Layer #9 model

#9) is rather reliable. Therefore, a model of only the topmost layer (Layer #9) could provide an ideal platform to investigate the effect of various parameters such as topography on the shaping of the CO2 plume, and also could be used for optimization purpose to achieve high efficiency sequestration while maintaining an affordable computational effort and cost.

**Figure 29.** Amplitude maps of Layer #9 from 1999 to 2008, from Singh et al. [26]

A reservoir model with dimension of 1600 m × 4900 m with varying thickness was constructed. It covers the portion of the Utsira formation where the plume, shown in Figure 31, is located. As mentioned earlier, the topography of this portion of the Utsira formation is accurately modeled based on seismic geological survey data (provided by Zhu and Lu of the University of Indiana [27]) with 50 m × 50 m mesh resolution. Because only Layer #9 is modeled, the thickness of the computational domain varies from 3.5 m to 26.3 m with an average thickness of 11.3 m. However, to accurately capture the accumulation and upward and lateral movement of CO2, 37 layers are used along the thickness. The topmost layer and the bottom two layers are designated to represent the low permeability shale, while the 34 layers in the middle are assigned the properties of mudstone. In the 3D Layer #9 model, permeability anisotropy is considered with west-east permeability of 2 mDarcy, north-south permeability of 10 mDarcy, and vertical permeability of 200 mDarcy. A 3D overview of the Layer #9 model is shown in

**Figure 30.** 3D overview and plan view of the 3D Layer #9 model of Utsira indicating feeder locations (black dot: main

Table 16 summarizes the hydrogeological properties of the Layer #9 model.

Figure 24.

336 CO2 Sequestration and Valorization

feeder; cyan square: secondary feeder)

It should be noted that in the 3D Layer #9 model, the source of CO2 is identified as "feeder" but not "injector" to emphasize that CO2 is supplied from the lower aquifer through leakage pathways rather than by direct injection. Since the actual CO2 injector is located at about 200 m under Layer #9, information on the injection rate recorded at the injector is not applicable for the CO2 feeders in the Layer #9 model. To determine the CO2 feeding rate to Layer #9, seismic surveys of CO2 distribution are used to obtain its volume under in situ conditions, and then converted to mass flow rate. Information of CO2 accumulative mass provided by Zhu and Lu [28] is summarized inTable 17.

**Table 17.** Accumulated CO2 mass in Layer #9, 1999-2008

Table 17 gives the CO2 accumulation in Layer #9 from 1999 to 2008. It can be seen that CO2 feeding rate to Layer #9 keeps on increasing for the recorded nine years. Recalling the analysis of secondary sealing effect given for the previous cases, it is the pressure gradient between the gaseous CO2 phase pressure at lower aquifer and the capillary pressure of the overlying shale layer that determines the breakthrough of CO2 and its flow rate. When breakthrough first occurs, the pressure gradient just breaks the equilibrium state, resulting in relatively low breakthrough mass flux to Layer #9. However, as more CO2 accumulates, the pressure gradient gradually increases and leads to increasing breakthrough mass flux, as depicted in Table 17. A nine-year average feeding rate of about 2.89 kg/s can be obtained from Table 17. Both the nine-year average value and the values in Table 17 have been considered in the simulations.

The significant north-tending plume finger is rather perplexing for regular pressure-gradient driven Darcy flow. Analysis suggests three possible explanations for the cause of the promi‐ nent north-tending CO2 fingering along the ridge: (1) significantly higher permeability applied to the ridge; (2) existence of north-south geological slope which enhances the buoyancy-drive migration along the ridge; and (3) existence of a secondary (or even multiple) CO2 pathways under the ridge. The hypothesis of significantly higher permeability at the ridge can be easily ruled out since no such evidence is obtained from the geological survey. It is still under debate whether geological slope should be considered when analyzing the CO2 migration in the Utsira formation. Chadwick and Noy's work [24] suggested two possible values of geological slope based on the seismic images of the cross-section of the Utsira formation, which are 8.2 m/km and 5.8 m/km.

The simulation time is set at nine years, which corresponds to the injection period of 1999~2008. CO2 plume migration at the topmost layer is examined for each year. Considering all the uncertainties mentioned above, a total of nine simulations are performed until good historymatching is obtained, as shown in Figure 31.

**Figure 31.** Simulated CO2 migration in Layer #9, 2000 ~ 2008

Both the 2D generalized Utsira formation model and the 3D detailed Utsira Layer #9 model have generated satisfactory simulation results, as seen by history-matching. In summary, five major conclusions can be reached, as follows. First, the simulations show that the permeability anisotropy should be accurately modeled. Vertical-to-horizontal anisotropy close to 10:1 is needed to accurately capture the upward migration of CO2. Horizontal anisotropy of 2:10 is needed to capture the northern spill of CO2 into the north-tending ridge. Second, a secondary feeder is likely to exist directly under the north-tending ridge to generate sufficient plume migration along the ridge. This suggests multiple pathways for CO2 breakthrough from the lower aquifer structure. Third, the fact that the injection gas is a CO2-methane mixture is very important in modeling since the presence of 2% methane enhances the buoyancy. Fourth, it is critical that time-dependent feeding of CO2 be modeled. This is consistent with the behavior of CO2 path flow breaking the capillary pressure barrier, as noted before for the secondarysealing effect in the case of the Mt.Simon formation. And finally, simulation results suggest strong mobility of gaseous CO2 under the caprock (shale) without major leakage, implying that the caprock serves quite well as the non-permeable CO2 barrier while exerting little resistance for the lateral flow of CO2 underneath.

The simulation studies of the three actual large-scale identified deep saline aquifers conclude this chapter. These studies have provided important insights and best practices for obtaining accurate simulations of GCS using TOUGH2. In Chapter titled," Optimization of CO2 Sequestration Saline Aquifers", innovative reservoir techniques and their optimization, for more efficient and secured SAGCS operations, are discussed.

## **6. Concluding remarks**

Table 17 gives the CO2 accumulation in Layer #9 from 1999 to 2008. It can be seen that CO2 feeding rate to Layer #9 keeps on increasing for the recorded nine years. Recalling the analysis of secondary sealing effect given for the previous cases, it is the pressure gradient between the gaseous CO2 phase pressure at lower aquifer and the capillary pressure of the overlying shale layer that determines the breakthrough of CO2 and its flow rate. When breakthrough first occurs, the pressure gradient just breaks the equilibrium state, resulting in relatively low breakthrough mass flux to Layer #9. However, as more CO2 accumulates, the pressure gradient gradually increases and leads to increasing breakthrough mass flux, as depicted in Table 17. A nine-year average feeding rate of about 2.89 kg/s can be obtained from Table 17. Both the nine-year average value and the values in Table 17 have been considered in the simulations. The significant north-tending plume finger is rather perplexing for regular pressure-gradient driven Darcy flow. Analysis suggests three possible explanations for the cause of the promi‐ nent north-tending CO2 fingering along the ridge: (1) significantly higher permeability applied to the ridge; (2) existence of north-south geological slope which enhances the buoyancy-drive migration along the ridge; and (3) existence of a secondary (or even multiple) CO2 pathways under the ridge. The hypothesis of significantly higher permeability at the ridge can be easily ruled out since no such evidence is obtained from the geological survey. It is still under debate whether geological slope should be considered when analyzing the CO2 migration in the Utsira formation. Chadwick and Noy's work [24] suggested two possible values of geological slope based on the seismic images of the cross-section of the Utsira formation, which are 8.2 m/km

The simulation time is set at nine years, which corresponds to the injection period of 1999~2008. CO2 plume migration at the topmost layer is examined for each year. Considering all the uncertainties mentioned above, a total of nine simulations are performed until good history-

Both the 2D generalized Utsira formation model and the 3D detailed Utsira Layer #9 model have generated satisfactory simulation results, as seen by history-matching. In summary, five major conclusions can be reached, as follows. First, the simulations show that the permeability

and 5.8 m/km.

338 CO2 Sequestration and Valorization

matching is obtained, as shown in Figure 31.

**Figure 31.** Simulated CO2 migration in Layer #9, 2000 ~ 2008

In this chapter, some key factors relevant to saline aquifer geological carbon sequestration (SAGCS) have been investigated. First, numerical simulations have been performed for completed/ongoing SAGCS projects in three large scale identified saline formations using the DOE numerical simulator TOUGH2. Before performing these studies, TOUGH2 was validated against the available analytical solutions and the benchmark numerical test cases. The three case studies of SAGCS in large-scale saline formations have provided important insights into the reservoir performance and sequestration uncertainties.

## **Acknowledgements**

Financial support for this work was provided by the Consortium for Clean Coal Utilization (CCCU) at Washington University in St. Louis, MO, USA.

## **Author details**

Zheming Zhang and Ramesh K. Agarwal

Department of Mechanical Engineering and Materials Science, Washington University in St. Louis, St. Louis, USA

## **References**

[12] Class H., Ebigbo A., Rainer H., Dahle H.K., Nordbotten J.M., Celia M.A., Audigane P., Darcis M., Ennis-King J., and Fan Y. A Benchmark Study on Problems Related to CO2 Storage in Geologic Formations. Computational Geosciences 2009; 13409-434.

**References**

340 CO2 Sequestration and Valorization

May 2013).

shop/(accessed 24 May 2013).

lem1.pdf(accessed 24 May 2013).

co2-workshop/Problem2.pdf(accessed 24 May 2013).

workshop/Problem3.pdf (accessed 24 May 2013).

atlasIV/ (accessed 24 May 2013).

[1] US Department of Energy. The 2012 United States Carbon Utilization and Storage At‐ las – Fourth Edition; 2008.http://www.netl.doe.gov/technologies/carbon\_seq/refshelf/