**The Role of Simulation Models in Monitoring Soil Organic Carbon Storage and Greenhouse Gas Mitigation Potential in Bioenergy Cropping Systems**

Manyowa N. Meki, James R. Kiniry, Kathrine D. Behrman, Meghan N. Pawlowski and Susan E. Crow

Additional information is available at the end of the chapter

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

**1. Introduction**

There is an increased demand on agricultural systems worldwide to provide food, fiber, and feedstock for the emerging bioenergy industry, raising legitimate concerns about the associ‐ ated impacts of such intensification on the environment [1, 2]. In the US, the revised National Renewable Fuel Standard program has mandated 30 billion gallons of renewable fuel by 2020 [3]. This aggressive promotion of the bioenergy industry in the US, and elsewhere in the world, is being fueled by several factors: the opportunity to reduce dependence on fossil fuels through renewable energy; the search for energy independence or security; its potential to reduce net greenhouse gas (GHG) emissions and hence provide mitigation options to combat climate change; and the possibility to improve farmers' incomes and rural economies, reduce national budget deficits and trade imbalances [4]. There is thus immense pressure on agricultural systems to provide the much needed feedstocks for this fast emerging industry. Of critical concern, however, is the fact that bioenergy feedstock production depends on finite resources such as land and water, and if not sustainably managed could have huge negative impacts on ecosystem services.

Of the many ecosystem services that could be impacted by the large-scale production of bioenergy feedstocks, soil organic carbon (SOC) plays a pivotal role in maintaining and enhancing the natural soil resource base and the GHG emission balance. According to the [5], periodic checks of SOC should be conducted to ensure that soil quality is not sacrificed for renewable bioenergy. Soil degradation causes carbon (C) loss and release of GHGs as a result

© 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.

of accelerated decomposition from land use change or unsustainable land management practices [6]. To date, there is little or no information that can be used to assess the long-term impacts of bioenergy production on SOC sequestration and storage.

In this chapter we discuss and examine the potential application of process-based simula‐ tion models in monitoring SOC and its GHG mitigation potential. The first section highlights the critical functions and benefits of SOC in maintaining soil quality, and hence productiv‐ ity and sustainability of cropping systems. In addition, we examine the ability of bioener‐ gy cropping systems to sequester SOC and hence mitigate climate change. In the second part, we focus on the need for periodic checks and long-term monitoring of SOC. We point out the challenges of measuring and quantifying SOC at the field level and review selected modeling studies that illustrate how process-based simulation models can be applied to monitor SOC and GHG emissions. A discussion of the importance of life cycle assess‐ ment (LCA) in evaluating a bioenergy cropping system's effect on the GHG emission balance or global warming potential (GWP) is presented. Finally, we offer our perspec‐ tives on how simulation models could provide the tools to track SOC for the Carbon Credit Markets (CCMs), and briefly examine government policies that are necessary to ensure success.

### **2. Soil organic matter and soil organic carbon – Functions and benefits**

Soils contain C in both organic and inorganic forms. In most soils the majority of C is held as SOC. The term soil organic matter (SOM) is used to describe the organic constituents in the soil, while SOC refers to the C occurring in the SOM [7]. The terms SOM and SOC are often used interchangeably because SOC is the main constituent of SOM, which is typically around 57% C. Besides carbon, SOM also contains other important elements such as oxygen (O, ~40%), nitrogen (N, ~3%), as well as smaller amounts of phosphorus (P), potassium (K), calcium (Ca), magnesium (Mg), sulfur (S) and micronutrients [8].

Soil organic matter is composed of dead plant and animal material at various stages of decomposition. These materials are mainly plant residues and leaf litter, animal wastes and remains, plant roots and exudates, and soil biota, e.g., microbes and earthworms. Soil biota are continually breaking down the SOM through physical and biochemical reactions, releasing C and nutrients back to the soil, and GHGs such as carbon dioxide (CO2), nitrous oxide (N2O) and methane (CH4) to the atmosphere. Soils also contain inorganic C in the form of biochar and carbonate minerals, such as calcite, dolomite and gypsum.

Soil organic carbon plays a very important function in the growth of plants through mainte‐ nance as well as improvement of many soil properties. These functions of SOC can be classified into three broad categories: biological, physical, and chemical. According to [9], these catego‐ ries are not static entities; instead, a change in one property will likely affect the other soil properties as well. An overview of the principal functions and benefits of SOM in soils is given in Figure 1.

The Role of Simulation Models in Monitoring Soil Organic Carbon Storage and Greenhouse Gas... http://dx.doi.org/10.5772/57177 253

increasingly highlighted.

Whether or not countries implement policies that stimulate increased C sequestration through bioenergy cropping, the importance of enhancing the soil resource base should be

**10** 39 [84], [84]

**13** 23 monitoring simulation

**16** 40 Foremost, Some suggest that

**16** 36 recently, recent

**17** should could

**17-18** 38-2 Finally, there is need for governments to put in place policies that stimulate increased C sequestration through bioenergy cropping, while ensuring that the soil resource base is

enhanced rather than degraded.

of accelerated decomposition from land use change or unsustainable land management practices [6]. To date, there is little or no information that can be used to assess the long-term

In this chapter we discuss and examine the potential application of process-based simula‐ tion models in monitoring SOC and its GHG mitigation potential. The first section highlights the critical functions and benefits of SOC in maintaining soil quality, and hence productiv‐ ity and sustainability of cropping systems. In addition, we examine the ability of bioener‐ gy cropping systems to sequester SOC and hence mitigate climate change. In the second part, we focus on the need for periodic checks and long-term monitoring of SOC. We point out the challenges of measuring and quantifying SOC at the field level and review selected modeling studies that illustrate how process-based simulation models can be applied to monitor SOC and GHG emissions. A discussion of the importance of life cycle assess‐ ment (LCA) in evaluating a bioenergy cropping system's effect on the GHG emission balance or global warming potential (GWP) is presented. Finally, we offer our perspec‐ tives on how simulation models could provide the tools to track SOC for the Carbon Credit Markets (CCMs), and briefly examine government policies that are necessary to ensure

**2. Soil organic matter and soil organic carbon – Functions and benefits**

magnesium (Mg), sulfur (S) and micronutrients [8].

and carbonate minerals, such as calcite, dolomite and gypsum.

Soils contain C in both organic and inorganic forms. In most soils the majority of C is held as SOC. The term soil organic matter (SOM) is used to describe the organic constituents in the soil, while SOC refers to the C occurring in the SOM [7]. The terms SOM and SOC are often used interchangeably because SOC is the main constituent of SOM, which is typically around 57% C. Besides carbon, SOM also contains other important elements such as oxygen (O, ~40%), nitrogen (N, ~3%), as well as smaller amounts of phosphorus (P), potassium (K), calcium (Ca),

Soil organic matter is composed of dead plant and animal material at various stages of decomposition. These materials are mainly plant residues and leaf litter, animal wastes and remains, plant roots and exudates, and soil biota, e.g., microbes and earthworms. Soil biota are continually breaking down the SOM through physical and biochemical reactions, releasing C and nutrients back to the soil, and GHGs such as carbon dioxide (CO2), nitrous oxide (N2O) and methane (CH4) to the atmosphere. Soils also contain inorganic C in the form of biochar

Soil organic carbon plays a very important function in the growth of plants through mainte‐ nance as well as improvement of many soil properties. These functions of SOC can be classified into three broad categories: biological, physical, and chemical. According to [9], these catego‐ ries are not static entities; instead, a change in one property will likely affect the other soil properties as well. An overview of the principal functions and benefits of SOM in soils is given

impacts of bioenergy production on SOC sequestration and storage.

success.

252 CO2 Sequestration and Valorization

in Figure 1.

**Figure 1.** The biological, physical, and chemical functions of soil organic matter (SOM). Source: adapted from [10].

*Biological functions –* Organic matter is a source of energy for many soil biota. Soil fauna play an important role in the initial breakdown of complex and large pieces of organic matter, making it easier for soil microorganisms to release C and plant nutrients during the process of decomposition. While microorganisms only make up less than 5% of SOM [11], they are the key driver of important processes such as decomposition, nutrient cycling, N fixation, and symbiotic plant relationships. These processes have direct and indirect effects on plant nutrient availability. It is estimated that 90-95% of the soil N, and about 2% of available P and S in natural ecosystems, come from organic matter.

2 *Chemical functions* – Humus is finely divided stable organic matter. Due to the large exposed surface area, humus contributes 30-70% of the cation exchange sites that adsorb plant nutrients that are eventually taken up by plants. Humus chelates micronutrients (usually iron, zinc, copper, or manganese) into soluble compounds, hence increasing their general mobility in soils and availability to plants. Humus also plays an important role in buffering the soil against any changes in acidity, salinity, and damage by pesticides. The exchange sites are also believed to help in cleaning up contaminated water by temporarily adsorbing heavy metal pollutants such as lead, cadmium, copper and zinc.

*Physical functions* – Plant roots, microbially-derived polysaccharides, and fungal hyphae cement soil particles into aggregates, thereby improving soil aeration, rate of water infiltration, and hence water holding capacity [12]. Build-up of organic matter lowers the soil bulk density, which allows the soil to be tilled with less horsepower, and improves plant rooting. Organic matter left on the soil surface, such as plant residues, also enhances rainfall capture, which reduces sediment and nutrient losses by runoff. Furthermore, organic matter shades the soil, which reduces excessive evapotranspiration and soil temperature in the summer, while keeping the soil warmer in winter.

A review of the literature shows that the productivity and sustainability of cropping systems are positively related to SOC content [13]. The ability to increase SOC in agricultural ecosys‐ tems is of interest both for sequestering atmospheric CO2 and for restoring organic-matter pools important to soil health [14]. The challenge is to identify soil management practices that promote SOC sequestration while ensuring productivity and profitability for farmers.

#### **2.1. Bioenergy crop yields, impacts on soil organic carbon storage and greenhouse gas mitigation potential**

Although most current bioenergy crop production is from annual row crops such as corn (*Zea mays* L.) and soybean [*Glycine max* (L.) Merr.], we do not discuss these here because they offer few or no environmental benefits over current agricultural practises [16]. Instead, we focus on second-generation bioenergy crops (Table 1) as identified by the US Department of Energy (DOE) [15, 16]: (1) the highly productive perennial warm-season grasses (WSGs) as emerging biofuel feedstocks for cellulosic ethanol production [e.g. switchgrass (*Panicum virgatum L.*), miscanthus (*Miscanthus giganteus*), Napiergrass (elephantgrass) (*Pennisetum purpureum* Schum.), and energy cane (*Saccharum officinarum* L.)]; and (2) the short-rotation woody crops (SRWCs) [e.g. hybrid poplar (*Populus* spp.) and willow (*Salix* spp.)].

In the US, average yields of selected WSGs ranged from 4.0 to 14.0 Mg ha-1 yr-l [15]. More recent results indicate the potential to reach 20.0 Mg ha-1 yr-1 in several locations. Energy cane and Napiergrass can produce yields from 12.2 to 32.4 Mg ha-1 yr-1. However, these species are very sensitive to frost and are limited to the southern regions of the US. Short rotation woody crop yields range from 10.0 to 17.0 Mg ha-1 yr-1. This is 2-3 times the yields normally achieved by traditional forest management. Observed positive attributes of these bioenergy crops include: being perennials, they require less tillage than conventional row crops (tillage is costly and is detrimental to maintenance of SOC stocks); their deep and extensive rooting patterns effec‐ tively exploit and extract water and soil nutrients deep-down the soil profile; they require little or no external nutrient inputs due to nutrient retention and efficient nutrient cycling between growing years or harvests; most are adapted to low-quality land (nutrient-depleted, compact‐ ed, poorly drained, and eroded); and they are tolerant to a host of abiotic stresses, such as drought, heat, cold, low pH, metal toxicity and salinity. The perennial grasses switchgrass and miscanthus have the additional advantage of the C4 pathway of photosynthesis, which is associated with high water use efficiency (WUE).

To date, there are a few empirical studies of SOC sequestration under both WSGs and SRWCs cropping systems. Therefore, the following discussion is based on data from studies carried out mainly in the US and Europe. More data is needed to provide a comprehensive under‐ standing for these cropping systems.

#### *2.1.1. Impacts of Warm Season Grasses*

Growing WSGs offers promise to displace fossil fuels and reduce net CO2 emissions through SOC sequestration [17]. When compared to row crops, the greater SOC storage under WSGs is attributed to greater aboveground biomass and root biomass concentration [18, 19]. Switch‐ grass produces about 6.7 Mg ha−1 yr−1 of total root biomass and half of this is found below the 30-cm soil depth [19], while about 10% is found below the 60-cm soil depth. Warm-season grass root systems, which often extend to nearly 3 m in the soil profile, enhance translocation


**Table 1.** Model bioenergy crops identified by the US Department of Energy (DOE).

pools important to soil health [14]. The challenge is to identify soil management practices that promote SOC sequestration while ensuring productivity and profitability for farmers.

Although most current bioenergy crop production is from annual row crops such as corn (*Zea mays* L.) and soybean [*Glycine max* (L.) Merr.], we do not discuss these here because they offer few or no environmental benefits over current agricultural practises [16]. Instead, we focus on second-generation bioenergy crops (Table 1) as identified by the US Department of Energy (DOE) [15, 16]: (1) the highly productive perennial warm-season grasses (WSGs) as emerging biofuel feedstocks for cellulosic ethanol production [e.g. switchgrass (*Panicum virgatum L.*), miscanthus (*Miscanthus giganteus*), Napiergrass (elephantgrass) (*Pennisetum purpureum* Schum.), and energy cane (*Saccharum officinarum* L.)]; and (2) the short-rotation woody crops

In the US, average yields of selected WSGs ranged from 4.0 to 14.0 Mg ha-1 yr-l [15]. More recent results indicate the potential to reach 20.0 Mg ha-1 yr-1 in several locations. Energy cane and Napiergrass can produce yields from 12.2 to 32.4 Mg ha-1 yr-1. However, these species are very sensitive to frost and are limited to the southern regions of the US. Short rotation woody crop yields range from 10.0 to 17.0 Mg ha-1 yr-1. This is 2-3 times the yields normally achieved by traditional forest management. Observed positive attributes of these bioenergy crops include: being perennials, they require less tillage than conventional row crops (tillage is costly and is detrimental to maintenance of SOC stocks); their deep and extensive rooting patterns effec‐ tively exploit and extract water and soil nutrients deep-down the soil profile; they require little or no external nutrient inputs due to nutrient retention and efficient nutrient cycling between growing years or harvests; most are adapted to low-quality land (nutrient-depleted, compact‐ ed, poorly drained, and eroded); and they are tolerant to a host of abiotic stresses, such as drought, heat, cold, low pH, metal toxicity and salinity. The perennial grasses switchgrass and miscanthus have the additional advantage of the C4 pathway of photosynthesis, which is

To date, there are a few empirical studies of SOC sequestration under both WSGs and SRWCs cropping systems. Therefore, the following discussion is based on data from studies carried out mainly in the US and Europe. More data is needed to provide a comprehensive under‐

Growing WSGs offers promise to displace fossil fuels and reduce net CO2 emissions through SOC sequestration [17]. When compared to row crops, the greater SOC storage under WSGs is attributed to greater aboveground biomass and root biomass concentration [18, 19]. Switch‐ grass produces about 6.7 Mg ha−1 yr−1 of total root biomass and half of this is found below the 30-cm soil depth [19], while about 10% is found below the 60-cm soil depth. Warm-season grass root systems, which often extend to nearly 3 m in the soil profile, enhance translocation

**2.1. Bioenergy crop yields, impacts on soil organic carbon storage and greenhouse gas**

(SRWCs) [e.g. hybrid poplar (*Populus* spp.) and willow (*Salix* spp.)].

associated with high water use efficiency (WUE).

standing for these cropping systems.

*2.1.1. Impacts of Warm Season Grasses*

**mitigation potential**

254 CO2 Sequestration and Valorization

of SOC to deeper layers, thereby reducing decomposition of SOC and promoting long-term C storage [19]. Data from a regional study across 42 paired sites in Minnesota, North Dakota and South Dakota showed that switchgrass stored about 2.0 Mg ha−1 of SOC in the 0- to 5-cm depth, 7.7 Mg ha−1 in the 30- to 60-cm depth, and 4.4 Mg ha−1 in the 60- to 90-cm depth [20]. The fraction of SOC stored in deeper (>30 cm) soil profiles is critical for long-term SOC sequestration because this fraction has longer residence times and slower turnover due to reduced microbial processes and fluctuations in soil water content and temperature [21]

In central and northwest Europe, annual increases of SOC for the conversion of croplands to miscanthus differed widely among studies and among sites within one study, ranging from 6.9 to 7.7 Mg C ha-1 yr-1 [22 - 24]. Data from soil samples of miscanthus fields and adjacent reference croplands at four different locations showed SOC change rates ranging from 2.6 to 2.8 Mg C ha-1 yr-1 [22], while [23] found change rates of 0.8 to 2.2 Mg C ha-1 yr-1. Since all the above studies were conducted in central and northwest Europe, this scatter emphasizes the impact of site-specific factors on SOC sequestration. Carbon sequestration rates in miscanthus fields were determined at 0.1-0.2 Mg C ha-1yr-1 [25], while a comparative study of C seques‐ tration in miscanthus and willow in Poland reported sequestration rates of 0.6 and 0.3 Mg C ha-1 yr-1, respectively [26]. Observations by [27] show that the establishment of miscanthus on grasslands may lead to SOC losses and that the co-benefit of SOC sequestration mainly occurs under miscanthus established on former croplands.

The amount of SOC sequestered by WSGs is a function of site-specific factors; soil texture, management practices, initial SOC levels, and climate [28, 29]. The sequestration rate can be greater in soils with initial low SOC levels than in those that are about to reach their C saturation level. On average, perennial WSGs can sequester between 0 and 3.0 Mg C ha−1 yr−1 in the 0- to 5-cm soil depth [28]. In a clay loam in Iowa, switchgrass stored 0.8 Mg C ha−1 yr−1 in plant and litter biomass, which was higher than that in cool-season grass and row crops [30]. The SOC content under switchgrass was 45% more than under row crops planted in a sandy loam in Alabama [31]. A comparative evaluation of SOC sequestration in croplands, and under 5 to 8 yr-old WSGs across 10 sites in Indiana, showed only moderate increases in SOC storage: WSGs - 22.4 g C kg−1, compared to 19.8 g C kg−1under croplands [32]. A review of 115 studies of grassland SOC sequestration from 17 countries including Australia, Brazil, Canada, New Zealand, the UK, and the US reported SOC sequestration rates ranging from 0.1 to 3.0 Mg C ha−1 yr−1 [33]. Site-specific rates depended on biome type and climate.

Carbon sequestration rates during early grassland establishment on disturbed lands ranged from 0.3 to 0.5 Mg C ha−1 yr−1 [33, 34]. For the short term, perennial grasses may offer better economic returns due to the quicker establishment and the annual harvests. Data on SOC sequestration under Napiergrass and energy cane is scarce. One study on Napiergrass in the southeast US showed an average total soil C increase of 3.2 Mg ha−1 after four years of cultivation [35].

While growing WSGs increases SOC concentration by 0.3 to 0.5 Mg ha−1 yr−1 [33, 34], crop residue removal (e.g., corn stover) reduces SOC sequestration by 1 to 1.5 Mg ha−1 yr−1 in the 0 to 30-cm soil depth [12]. Growing WSGs thus increases SOC concentration while providing alternative feedstocks for cellulosic ethanol production, unlike crop residue removal, which results in large losses of SOC. The increase in SOC with WSGs not only plays a major role in mitigating GHG emissions, it also enhances ecosystem functions.

#### *2.1.2. Impacts of Short Rotation Woody Crops*

Information on SOC storage has not been widely reported for SRWCs. The limited data available so far comes from some few studies in the US, Canada and Europe. Two US studies were established in Minnesota to specifically address the issue of SOC sequestration under SRWCs [36, 37]. In the study by [36], soil C content was monitored in poplar plantations established on previously tilled agricultural prairie land, and compared to adjacent control grass and arable fields, while in the study by [37], poplar plantations were established on grass and arable lands. In both studies, there was a net loss of SOC over the first 6–12 years of initial establishment of the plantations. The SOC loss was largely from the surface 30 cm of soil, suggesting that it was due to enhanced decomposition. While [37] observed no increase in SOC under the poplar plantations, [36] reported an average rate of SOC increase of 1.6 Mg C ha-1 yr-1. Other research also indicates that SRWCs may result in high erosion and runoff rates during the first year(s) of establishment.

As with perennial WSGs, SRWCs have a huge potential to increase SOC due to the abundant above- and belowground biomass input. For example, one study by [38] showed that willow trees produce finer root biomass than corn (5.8 Mg ha-1 versus 0.9 Mg ha-1 for corn). On average, SRWCs can store between 0 and 1.6 Mg C ha−1 yr−1 in the 0- to 100-cm soil depth [39]. On relatively fertile soils in Canada, willow trees stored more SOC than either corn or switchgrass after 4 yr of establishment [40]. A determination of C sequestration in forests and willow cultivations showed sequestration rates at about 0.2 Mg C ha-1yr-1 [25]. In the US Midwest region, gains in SOC under SRWCs were greater than in croplands when the initial SOC concentration was low [41].

In the UK, estimates of SOC sequestration rates under SRWCs were based on willow and poplar plantations that were coppiced for varying time periods at three sites [42]. The average annual increase in SOC under SRWCs was estimated at between 1.2% and 2.2% [43]. In a study of three mixed poplar, aspen and willow plantation sites across Germany, substantial increases in SOC (0.1 to 0.6 Mg C ha-1 yr-1) in the upper surface soil were observed after seven to nine years [44]. For much of Europe, research investigating potential SOC sequestration under SRWCs is still in its infancy.

According to [41], SRWCs should be grown in marginal and degraded lands where there will be greater accumulation of SOC rather than in croplands or natural forests. However, [45] argue that perennial WSGs may be more suitable for C sequestration since SRWCs take time for canopy closure, making soil more prone to SOC loss. If a management strategy increases the above- and belowground biomass, then a net gain in SOC sequestration can result.

## *2.1.3. Greenhouse gas mitigation potential*

The amount of SOC sequestered by WSGs is a function of site-specific factors; soil texture, management practices, initial SOC levels, and climate [28, 29]. The sequestration rate can be greater in soils with initial low SOC levels than in those that are about to reach their C saturation level. On average, perennial WSGs can sequester between 0 and 3.0 Mg C ha−1 yr−1 in the 0- to 5-cm soil depth [28]. In a clay loam in Iowa, switchgrass stored 0.8 Mg C ha−1 yr−1 in plant and litter biomass, which was higher than that in cool-season grass and row crops [30]. The SOC content under switchgrass was 45% more than under row crops planted in a sandy loam in Alabama [31]. A comparative evaluation of SOC sequestration in croplands, and under 5 to 8 yr-old WSGs across 10 sites in Indiana, showed only moderate increases in SOC storage: WSGs - 22.4 g C kg−1, compared to 19.8 g C kg−1under croplands [32]. A review of 115 studies of grassland SOC sequestration from 17 countries including Australia, Brazil, Canada, New Zealand, the UK, and the US reported SOC sequestration rates ranging from 0.1 to 3.0 Mg C

Carbon sequestration rates during early grassland establishment on disturbed lands ranged from 0.3 to 0.5 Mg C ha−1 yr−1 [33, 34]. For the short term, perennial grasses may offer better economic returns due to the quicker establishment and the annual harvests. Data on SOC sequestration under Napiergrass and energy cane is scarce. One study on Napiergrass in the southeast US showed an average total soil C increase of 3.2 Mg ha−1 after four years of

While growing WSGs increases SOC concentration by 0.3 to 0.5 Mg ha−1 yr−1 [33, 34], crop residue removal (e.g., corn stover) reduces SOC sequestration by 1 to 1.5 Mg ha−1 yr−1 in the 0 to 30-cm soil depth [12]. Growing WSGs thus increases SOC concentration while providing alternative feedstocks for cellulosic ethanol production, unlike crop residue removal, which results in large losses of SOC. The increase in SOC with WSGs not only plays a major role in

Information on SOC storage has not been widely reported for SRWCs. The limited data available so far comes from some few studies in the US, Canada and Europe. Two US studies were established in Minnesota to specifically address the issue of SOC sequestration under SRWCs [36, 37]. In the study by [36], soil C content was monitored in poplar plantations established on previously tilled agricultural prairie land, and compared to adjacent control grass and arable fields, while in the study by [37], poplar plantations were established on grass and arable lands. In both studies, there was a net loss of SOC over the first 6–12 years of initial establishment of the plantations. The SOC loss was largely from the surface 30 cm of soil, suggesting that it was due to enhanced decomposition. While [37] observed no increase in SOC under the poplar plantations, [36] reported an average rate of SOC increase of 1.6 Mg C ha-1 yr-1. Other research also indicates that SRWCs may result in high erosion and runoff rates

As with perennial WSGs, SRWCs have a huge potential to increase SOC due to the abundant above- and belowground biomass input. For example, one study by [38] showed that willow trees produce finer root biomass than corn (5.8 Mg ha-1 versus 0.9 Mg ha-1 for corn). On average,

ha−1 yr−1 [33]. Site-specific rates depended on biome type and climate.

mitigating GHG emissions, it also enhances ecosystem functions.

*2.1.2. Impacts of Short Rotation Woody Crops*

during the first year(s) of establishment.

cultivation [35].

256 CO2 Sequestration and Valorization

In contrast to row crops, evidence shows that biofuels from perennial WSGs and SRWCs have a positive GHG mitigation potential. Carbon sequestration associated with bioenergybased cropping systems could potentially offset 1 to 2 Pg C yr−1 globally [46]. Three GHGs: CO2, N2O and CH4 contribute directly to the GHG emission balance of cropping systems [47]. Nitrous oxide is produced in soils mostly through nitrification and denitrification. The majority of agricultural soils are minor sinks or sources of CH4 [48]. Carbon dioxide emissions come from diesel fuel used in farm operations and during the manufacture of agricultural machinery, seed, and agrochemicals. Accounting for the GHG emission balance of bioenergy cropping systems entails evaluating their net impact on all emissions associat‐ ed directly and indirectly to bioenergy feedstock production [1].

In the US, 60 million (M) ha of land is potentially available for conversion to bioenergy plantations. The conversion to bioenergy crops could mitigate C emissions at a rate 5.4 Mg C ha−1 yr−1. On 60 M ha, that would represent 324 Tg C yr−1, which is equivalent to a 20% reduction from current fossil-fuel CO2 emissions [39]. In one study it was estimated that SRWCs sequester soil C at an average rate of 0.7 Mg C ha-1 yr-1 (Table 2). In addition, SRWCs could also generate substantial reductions in fertilizer and fuel use, reducing CO2 and N2O emissions, for a net GHG mitigation potential of 1.1 Mg C ha-1 yr-1. According‐ ly, as much as 40 M ha of highly eroded, degraded or mined lands in the US could be planted with SRWCs with limited negative impact on the production of key food and fiber crops [49].


Positive numbers depict removal of greenhouse gases from atmosphere or prevented emissions. Source: [49].

**Table 2.** Carbon sequestration and greenhouse gas (GHG) mitigation potential in short rotation woody crops (SRWCs) – based on 35 field-based observations.

Fluxes of N2O and CH4 of SRWCs remain unknown. Uncertainties regarding N2O emissions from nitrification and denitrification complicate the estimation of a system's GHG mitigation potential or GWP [1]. However, estimates of CH4 consumption by soils are small, generally found to have little impact on net GWP [48]. One study across Europe noted much lower N2O emissions in forested versus agricultural lands [50]. A review of the potential of SOC seques‐ tration and CO2 offset by SRWCs in the US [51] reported C sequestration rates ranging from 0 to 1.6 Mg C ha-1 yr-1, and a GHG mitigation potential of 5.4 Mg C ha-1 yr-1. In contrast, little difference was found in the N2O emissions of annual crops versus those of poplar plantations [52]. Other researchers have estimated the bioenergy displacement of fossil fuels from SRWCs being as much as 4.9 to 5.5 Mg C ha-1 yr-1 [49]. In the majority of studies, N2O emissions are the main source or driver of the cropping system's GWP. The assumption is that C sequestration by bioenergy crops of around 0.25 Mg C ha-1 yr-1 makes bioenergy neutral in terms of the GHG emission balance [53]. Quantifying these factors is necessary to determine the net effect of bioenergy cropping systems on soil C sequestration and GHG emissions.

#### **2.2. Monitoring soil organic carbon storage and greenhouse gas fluxes**

Soil organic carbon is the single most important indicator of soil quality [54, 55]. In addition to the importance of SOC as an indicator of soil quality, accurate measurements of soil C are also crucial in the assessment of the GHG emission balance of bioenergy systems. In general, methods for the determination of soil C involve direct and indirect approaches. Direct methods employ field and laboratory measurements of soil C stocks or CO2 fluxes over time, while indirect methods combine data collected in direct soil C measurements and process-based simulation models to predict soil C changes temporally and spatially.

#### *2.2.1. Field monitoring*

The International Climate Change Treaty (Kyoto Protocol) recognized removal of CO2 from the atmosphere by plants as a valid approach to mitigating climate change [56] and identified the need to conduct long-term monitoring of C stocks under various land uses. While the methods for analyzing the SOC concentration of a given soil sample are well established and easily carried out with high precision and negligible analytical error [57], field documentation of SOC changes faces many challenges because of the heterogeneity of soils, environmental conditions, land use history, sampling methods and analytical errors [58]. Here we highlight some of the challenges.

**GHG category Switch to short rotation woody crops**

Soil carbon **0.7** (-2.1 to 3.6) Nitrous oxide **0.2** Methane **no data** Process and upstream emissions **0.2** Sum of GHGs **1.2** (-1.6 to 4.0) Maximum U.S. applicable area, M ha **40**

bioenergy cropping systems on soil C sequestration and GHG emissions.

**2.2. Monitoring soil organic carbon storage and greenhouse gas fluxes**

simulation models to predict soil C changes temporally and spatially.

Positive numbers depict removal of greenhouse gases from atmosphere or prevented emissions. Source: [49].

**Table 2.** Carbon sequestration and greenhouse gas (GHG) mitigation potential in short rotation woody crops (SRWCs)

Fluxes of N2O and CH4 of SRWCs remain unknown. Uncertainties regarding N2O emissions from nitrification and denitrification complicate the estimation of a system's GHG mitigation potential or GWP [1]. However, estimates of CH4 consumption by soils are small, generally found to have little impact on net GWP [48]. One study across Europe noted much lower N2O emissions in forested versus agricultural lands [50]. A review of the potential of SOC seques‐ tration and CO2 offset by SRWCs in the US [51] reported C sequestration rates ranging from 0 to 1.6 Mg C ha-1 yr-1, and a GHG mitigation potential of 5.4 Mg C ha-1 yr-1. In contrast, little difference was found in the N2O emissions of annual crops versus those of poplar plantations [52]. Other researchers have estimated the bioenergy displacement of fossil fuels from SRWCs being as much as 4.9 to 5.5 Mg C ha-1 yr-1 [49]. In the majority of studies, N2O emissions are the main source or driver of the cropping system's GWP. The assumption is that C sequestration by bioenergy crops of around 0.25 Mg C ha-1 yr-1 makes bioenergy neutral in terms of the GHG emission balance [53]. Quantifying these factors is necessary to determine the net effect of

Soil organic carbon is the single most important indicator of soil quality [54, 55]. In addition to the importance of SOC as an indicator of soil quality, accurate measurements of soil C are also crucial in the assessment of the GHG emission balance of bioenergy systems. In general, methods for the determination of soil C involve direct and indirect approaches. Direct methods employ field and laboratory measurements of soil C stocks or CO2 fluxes over time, while indirect methods combine data collected in direct soil C measurements and process-based

The International Climate Change Treaty (Kyoto Protocol) recognized removal of CO2 from the atmosphere by plants as a valid approach to mitigating climate change [56] and identified

Mg C ha-1 yr-1

Number of observations **35**

– based on 35 field-based observations.

258 CO2 Sequestration and Valorization

*2.2.1. Field monitoring*

Variability in soil C estimates can be due to a number of factors or conditions: uncertainties in bulk density measurements as impacted by above and below-ground biomass inputs; rock and stone contents; compaction; sediment loss; and deposition by wind and water erosion. These conditions can drastically affect the accuracy of soil C sequestration estimates, implying that a relatively large number of soil samples will be required to obtain more reliable estimates. For example, [59] suggest that 20 to 30 soil samples (250 cm3 ) will be needed for estimates of mean SOC within ±10% at the 95% confidence level, while [60] estimate that more than 100 samples (2.54 cm cores) would be needed to detect a 2–3% change in soil C under switchgrass.

The small annual incremental changes of soil C are difficult if not impossible to detect by present analytical methods. According to [61], noticeable soil C changes require a longer time period (e.g., at least 20 years). Statistically significant differences in SOC were observed after 5 to 10 yr [34, 62 - 64]. Coleman et al. [41] collected soil samples from 27 study sites across north central US to compare the soil C of short rotation poplar plantations to adjacent agricultural crops and woodlots. Soil organic carbon ranged from 20 to more than 160 Mg ha-1 across the sampled sites. Woodlots consistently contained greater SOC than the other crops, especially at depth. Surprisingly, [41] observed little difference between paired poplar and switchgrass. Furthermore, there was no evidence of changes in poplar SOC relative to adjacent agricultural soils when considered for stand ages up to 12 years.

A review of the different approaches to measurement and monitoring of SOC stocks shows that the main challenge is in designing an efficient, cost-effective sampling and SOC stock estimation system [65]. Current field sampling methods rely on a set of measurements that are extrapolated to represent a given geographic area. However, such intensive sampling is time consuming and costly. Furthermore, quantifying SOC changes at national or regional scales requires even higher sampling densities [66]. Conant et al. [65] are however optimistic that alternative measurement techniques, such as the infrared probe used by [67] or gamma-ray spectroscopy [68], could allow field-wide measurements of SOC stocks, even though the technologies currently require extensive site-specific instrument calibration. In addition, eddy covariance methods have improved greatly over the last two decades, and if coupled with appropriate and thoroughly tested terrestrial biosphere models, could be an alternative method to identify portions of the net ecosystem exchange associated with SOC [69]. Until remotely sensed measurements are calibrated to subsurface information, remotely sensed SOC inventories should only be considered as minimum estimates [70].

Given the many benefits and functions of SOC, a number of soil monitoring networks (SMN) that target soil C and other related data have been implemented or are under development in a number of countries throughout the world. For example, the European SMN – ENVASSO (ENVironmental ASsessment of Soil for mOnitoring; [71] has partners from 25 countries and a total of 33 334 soil C monitoring sites. The density of the European networks is high, with a median coverage of 300 km² for each monitoring site [58]. Although most SMNs are designed primarily for documenting soil C stocks and related data, their objectives do vary from country to country. For a more detailed review of the concepts, design and objectives of some of the SMNs across the world, see [72].

In the US there has been considerable investment in soil C surveys in accordance with the Energy Independence and Security Act (EISA) [73]. Two such programs are (1) the U.S. Geological Survey (USGS) carbon sequestration program (www.usgs.gov/climate\_landuse/ carbon\_seq/). Two USGS national-scale assessments (subsurface and terrestrial) covering all 50 states and primary ecosystems are currently under way, and will provide the most com‐ prehensive accounting of C storage potential in the US. The second program (2) is the USDA Natural Resources Conservation Service (NRCS) Rapid Assessment of US Soil Carbon (RaCA) (http://soils.usda.gov/survey/raca/). In general, the two programs were initiated to develop statistically reliable quantitative estimates of amounts and distribution of C stocks for US soils under various lands, differing agricultural management; to provide data to support model simulations of soil C change related to land use change, agricultural management, conserva‐ tion practices, and climate change; and to provide a scientifically and statistically defensible inventory of soil C sequestration and GHG emissions for the US.

In Canada, a group of energy industries formed the GEMCo (Greenhouse Gas Emissions Management Consortium) in recognition of the potential of soil C sequestration as a possible GHG offset mechanism, together with the need to develop soil C measuring and verification protocols for C trading [74, 75]. The GEMCo, in partnership with Canadian soil conservation associations and government agencies, went on to initiate work aimed at monitoring and documenting soil C sequestration at the field level [69].

Some national [76, 77] and global [78, 79, 80] scale assessments of soil C stocks have been made. The World Bank-Global Environmental Fund's sustainable development project is supporting similar soil C monitoring initiatives in developing countries, e.g., Mexico [69]. In the US, the USDA NRCS soil databases; (1) State Soil Geographic (STATSGO) [81], provides generalized soil data for the entire US, and (2) the Soil Survey Geographic (SSURGO) provides detailed soil information for selected US counties [82]. Both these databases can be used to estimate SOC stocks [69]. Although at a lower resolution, the FAO's Soil Map of the World [79, 83] provides global estimates of SOC stocks.

Globally, a major issue with SMNs or any other soil C monitoring program is that there is no consensus on soil C measurement, monitoring, and verification methods. Post et al. [69] stress that there is an urgent need to develop robust, science-based, flexible, and agreed-upon practical protocols for monitoring and verifying temporal and spatial changes in soil C.

#### *2.2.2. Application of simulation models*

Computer simulation models can complement and extend the applicability of information collected in field trials [84] and can be applied to obtain realistic assessments of bioenergy feedstock production effects on SOC storage and the GHG emission balance. Combining measurement of SOC with models provides more flexibility and reliability because the model can be continually updated with new field sample data, and will be able to account for both C sequestration and GHG mitigation potential through LCA.

a total of 33 334 soil C monitoring sites. The density of the European networks is high, with a median coverage of 300 km² for each monitoring site [58]. Although most SMNs are designed primarily for documenting soil C stocks and related data, their objectives do vary from country to country. For a more detailed review of the concepts, design and objectives of some of the

In the US there has been considerable investment in soil C surveys in accordance with the Energy Independence and Security Act (EISA) [73]. Two such programs are (1) the U.S. Geological Survey (USGS) carbon sequestration program (www.usgs.gov/climate\_landuse/ carbon\_seq/). Two USGS national-scale assessments (subsurface and terrestrial) covering all 50 states and primary ecosystems are currently under way, and will provide the most com‐ prehensive accounting of C storage potential in the US. The second program (2) is the USDA Natural Resources Conservation Service (NRCS) Rapid Assessment of US Soil Carbon (RaCA) (http://soils.usda.gov/survey/raca/). In general, the two programs were initiated to develop statistically reliable quantitative estimates of amounts and distribution of C stocks for US soils under various lands, differing agricultural management; to provide data to support model simulations of soil C change related to land use change, agricultural management, conserva‐ tion practices, and climate change; and to provide a scientifically and statistically defensible

In Canada, a group of energy industries formed the GEMCo (Greenhouse Gas Emissions Management Consortium) in recognition of the potential of soil C sequestration as a possible GHG offset mechanism, together with the need to develop soil C measuring and verification protocols for C trading [74, 75]. The GEMCo, in partnership with Canadian soil conservation associations and government agencies, went on to initiate work aimed at monitoring and

Some national [76, 77] and global [78, 79, 80] scale assessments of soil C stocks have been made. The World Bank-Global Environmental Fund's sustainable development project is supporting similar soil C monitoring initiatives in developing countries, e.g., Mexico [69]. In the US, the USDA NRCS soil databases; (1) State Soil Geographic (STATSGO) [81], provides generalized soil data for the entire US, and (2) the Soil Survey Geographic (SSURGO) provides detailed soil information for selected US counties [82]. Both these databases can be used to estimate SOC stocks [69]. Although at a lower resolution, the FAO's Soil Map of the World [79, 83]

Globally, a major issue with SMNs or any other soil C monitoring program is that there is no consensus on soil C measurement, monitoring, and verification methods. Post et al. [69] stress that there is an urgent need to develop robust, science-based, flexible, and agreed-upon practical protocols for monitoring and verifying temporal and spatial changes in soil C.

Computer simulation models can complement and extend the applicability of information collected in field trials [84] and can be applied to obtain realistic assessments of bioenergy feedstock production effects on SOC storage and the GHG emission balance. Combining

inventory of soil C sequestration and GHG emissions for the US.

documenting soil C sequestration at the field level [69].

provides global estimates of SOC stocks.

*2.2.2. Application of simulation models*

SMNs across the world, see [72].

260 CO2 Sequestration and Valorization

Although several soil C and GHGs emission models have been developed for conventional agricultural and forest systems, most of them have not been fully parameterized and effectively tested for application on WSGs and SRWCs. Detailed documentation on the main SOC models, their main features and comparative performance is widely available in the literature [85, 86, 87, 88, 89]. We focus on the applicability of a select number of process-based models that have been, or are currently being applied and hence could be adopted or adapted to monitor SOC sequestration and the GHG emission balance of bioenergy crops.

The amount of SOC in a given soil can increase or decrease depending on numerous factors, including climate, vegetation type, nutrient availability, disturbance, land use, and manage‐ ment practice. Most common models that describe this dynamic behavior or turnover of SOC follow the approach by [90], who proposed that soil C dynamics can be represented by two pools comprising one C pool and one residue pool, and can be described as follows:

$$\mathbf{dC\_s/dt} = \mathbf{hC\_i} \cdot \mathbf{kC\_s} \tag{1}$$

where Cs is soil carbon (Mg ha−1), t is time (yr), Ci is carbon inputs from plant residues and roots (Mg ha−1), h is a humification factor and k is the apparent soil decomposition rate (yr−1). This model, while simple, is adequate for hypothesis development and testing, and for fieldlevel data interpretation [91]. However, as more information on SOC dynamics became available, new multi-compartment models accounting for variations in turnover rates for different SOC pools were developed [92]. For example, the most widely applied soil C simulation models, CENTURY [93] and RothC [94], divide SOM into different pools with varying turnover rates (Table 3). This partitioning is supposed to reflect the average biochem‐ ical and physical properties of SOM in the pools (e.g., metabolic, structural and recalcitrant) or the accessibility of the SOM to soil microbes or catalytic enzymes, or constraints imposed on decomposition by environmental conditions [95].

In the CENTURY model, SOC is divided into the active, slow and passive pools, with mean residence times (MRTs, i.e., the time required for half of the SOM to decompose) of 1.5, 25, and about 1000 yr, respectively. The organic inputs are separated into metabolic (readily decom‐ posable; MRT of 0.1–1 yr) and structural (difficult to decompose; MRT of 1–5 yr) pools as a function of the lignin:N ratio [96]. A more detailed discussion of the SOC pools is given in [95]. It is, however, important to note that the C pool compartments are only conceptual, and have not been verified experimentally [96]. Despite this drawback, compartmental models have proved to be reasonably good in simulating changes in SOC stocks.

Although the CENTURY model was originally developed for grasslands [97], versions of the model have successfully been developed to simulate crops [98] and forest systems [99]. Among the nine SOC models that were evaluated against long-term datasets at a NATO Advanced Research Workshop held at Rothamsted, UK, in 1995, the CENTURY model had the best potential for adaptation to bioenergy cropping systems because of its integrated plant-soil approach and the availability of specific forestry subroutines [88]. A similar comparative study of six simulation models; RothC, CENTURY, EPIC [100], DNDC (DeNitrification DeCompo‐ sition) [101], ECOSYS [102] and SOCRATES (Soil Organic Carbon Reserves And Transforma‐ tions in agro-EcoSystems) [103] was conducted by [87]. Based on the model simulated outputs, the SOCTRATES model was best able to describe temporal (soil C changes observed in longterm experiments) and spatial SOC dynamics, and hence could be applied to scale-up soil C sequestration at the regional scale [87]. In another study, the model successfully predicted SOC change at eighteen long-term crop, pasture and forestry trials from North America, Europe and Australia [104]. A major advantage of the SOCRATES model is that it requires minimal data inputs, and is specifically designed to examine the impact of land use and land use change on SOC storage.

Despite sharing several basic concepts with the CENTURY model, the RothC model's applic‐ ability in monitoring SOC stocks in bioenergy cropping systems is rather limited due to the fact that it is solely concerned with soil processes, and does not contain a sub-model for plant production. In addition, it does not calculate annual returns of plant C to the soil from aboveground yields [88].The model has however been successful applied to describe SOC evolution under several long-term experiments, across a variety of land use and climate types [105].


DPM, decomposable plant material; BIO, microbial biomass; RPM, resistant plant material; HUM, humified organic matter; IOM, inert organic matter; POM, particulate organic matter; OM, organic matter. Source: [95].

**Table 3.** The pools of SOM defined according to their mean residence times (MRTs) and corresponding compound classes.

Due to the general applicability of the CENTURY model, a number of process-based simulation models such as EPIC, CropSyst [106] and DAYCENT [107], either adopted or adapted the SOC modeling algorithms from the CENTURY model. EPIC has since evolved into a comprehensive agro-ecosystem model with a SOC module capable of simulating SOC dynamics in a wide range of plant species, including crops, native grasses and trees [100]. Recently, [108] success‐ fully applied the EPIC model to simulate the production and environmental effects of perennial biomass feedstocks grown on marginal lands. The comparative assessment of six alternative cropping systems over 20 years showed that once well established, perennial herbaceous grasses have a direct GHG emissions mitigation capacity (-2.3±0.1 Mg C ha-1 yr-1) which is higher than that of conventional crops. Meki et al. [84] applied the sister model, ALMANAC [109], to evaluate sustainable energy sorghum biomass harvest thresholds and tillage effects on SOC and bulk density. Model results showed that 75% of energy sorghum biomass can be removed from a continuous no till system without any detrimental effects on SOC storage. This level of biomass removal, however, significantly increased soil bulk density, a critical indicator of soil compaction.

potential for adaptation to bioenergy cropping systems because of its integrated plant-soil approach and the availability of specific forestry subroutines [88]. A similar comparative study of six simulation models; RothC, CENTURY, EPIC [100], DNDC (DeNitrification DeCompo‐ sition) [101], ECOSYS [102] and SOCRATES (Soil Organic Carbon Reserves And Transforma‐ tions in agro-EcoSystems) [103] was conducted by [87]. Based on the model simulated outputs, the SOCTRATES model was best able to describe temporal (soil C changes observed in longterm experiments) and spatial SOC dynamics, and hence could be applied to scale-up soil C sequestration at the regional scale [87]. In another study, the model successfully predicted SOC change at eighteen long-term crop, pasture and forestry trials from North America, Europe and Australia [104]. A major advantage of the SOCRATES model is that it requires minimal data inputs, and is specifically designed to examine the impact of land use and land use change

Despite sharing several basic concepts with the CENTURY model, the RothC model's applic‐ ability in monitoring SOC stocks in bioenergy cropping systems is rather limited due to the fact that it is solely concerned with soil processes, and does not contain a sub-model for plant production. In addition, it does not calculate annual returns of plant C to the soil from aboveground yields [88].The model has however been successful applied to describe SOC evolution under several long-term experiments, across a variety of land use and climate types [105].

Structural 2–5 100–200 Polysaccharides

DPM POM

Slow RPM 15–100 10–25 Lignified tissues

Passive HUM 500–5000 7–10 Humic substances

DPM, decomposable plant material; BIO, microbial biomass; RPM, resistant plant material; HUM, humified organic matter;

**Table 3.** The pools of SOM defined according to their mean residence times (MRTs) and corresponding compound

IOM, inert organic matter; POM, particulate organic matter; OM, organic matter. Source: [95].

IOM Clay: OM complexes

Amino acids Starch

Polysaccharides

Waxes Polyphenols

Biochar

**Residue type Century RothC Residence time (yr) C:N Compounds** Litter Metabolic DPM 0.1–0.5 10–25 Simple sugars

SOM Active BIO 1–2 15–30 Living biomass

on SOC storage.

262 CO2 Sequestration and Valorization

classes.

By integrating outputs from the REAP model [111] with environmental impact estimates derived from the EPIC model, [110] were able to evaluate the environmental and economic impacts of increased ethanol production from switchgrass. Highlights from this study indicated that commercial-scale switchgrass production could be associated with reduced erosion and reduced nutrient losses, but the projected increases in N fertilizer application and the associated N2O emissions could offset soil C sequestration benefits and result in substantial increases in total GHGs emissions from the agricultural sector.

In many European countries, soil C and GHG simulation relies on the Yasso07 soil C model [112]. This model describes litter decomposition and soil C cycle based on the chemical quality of organic matter and the climatic conditions. Total soil C stocks are estimated for the top-most 1 m soil layer. Since Yasso07 is a reasonably new model [113], not many studies on the model's applicability to GHG gas inventory purposes are available. The strategy when developing the Yasso07 model was to achieve a soil C model that had fewer input data requirements, thus ensuring wider applicability [114]. A recent study in southern Finland by [112] showed the Yasso07 model accurately predicting the long-term soil C stocks of a coniferous forest soil (67.0 Mg ha-1), which accumulated at an average rate of 0.12 ± 0.06 Mg ha-1 yr-1. Previously, [115] evaluated the performance of the Yasso07 model in predicting the soil C stock changes in afforestation and deforestation sites. Overall, there was a good agreement between simulated and measured soil C stocks at most sites.

Garten and Wullschleger [60] used a simple, parameterized two-compartment model for SOC dynamics to predict the potential for soil C sequestration under switchgrass in the southeastern US. Model simulations indicated a measurable and verifiable recovery of soil C (~12% increase) on degraded lands through one decade of switchgrass production. A sensitivity analysis of the model indicated that switchgrass potential to sequester C depends on initial soil C inventories, prevailing climate, soil type, and site management.

The Terrestrial Ecosystem Model (TEM) is a process-based model that describes C and N dynamics of plants and soils for terrestrial ecosystems (http://ecosystems.mbl.edu/tem/). The

TEM uses spatially referenced information on climate, elevation, soils, vegetation and water availability as well as soil- and vegetation-specific parameters to make monthly estimates of C and N fluxes and pool sizes of terrestrial ecosystems. Conversion of agricultural lands to bioenergy crops consequently affects the ecosystem C balance. Qin et al. [116] applied the TEM model to evaluate the impacts of this change on the C balance, assuming several land use change scenarios from corn, soybean, and wheat to bioenergy crops of switchgrass and miscanthus. Overall, SOC increased significantly when land use changed from corn, soybean or wheat, to either switchgrass or miscanthus, reaching an average of ~100 Mg C ha-1 in the two bioenergy crops.

Relatively simple statistical models can also be useful in estimating SOC losses associated with land conversion. Anderson–Teixeira, et al. [117] used reduced ANOVA and the forced- and free-intercept regression statistical models to estimate and compare SOC losses associated with land conversion and SOC sequestration rates under corn, sugarcane, switchgrass, miscanthus and mixed native grass. Their predicted SOC accumulation rates under switchgrass, miscan‐ thus and mixed native grass ranged from 0.1- 1.0 Mg C ha-1 yr-1 in the top 30 cm, while SOC losses under corn and sugarcane ranged from 3.0–8.0 Mg C ha-1 yr-1. These results are consistent with data from numerous field studies reported elsewhere in the literature [20, 28, 31, 40].

Life-cycle assessments are commonly used to determine whether a biomass biofuel system results in a net reduction in GHG emissions or an improved energy source. In LCAs, all input and output data in all phases of the product's life cycle including biomass production, feedstock storage, feedstock transportation, biofuel production, biofuel transportation and final use are required. Because field measurements of bioenergy LCA components are difficult to carry out, simulation models have been applied effectively to obtain realistic estimates of the bioenergy crop's GWP or GHG mitigation potential.

Standard and useful tools for LCA are multidimensional spreadsheet models such as the GREET [118] and GHGenius [119] software models, which were specifically designed to address full LCA of biofuels. These spreadsheet models have advantages in that they are userfriendly, publicly available, straightforward, and relatively transparent. In a study commis‐ sioned by the US Department of Energy, Wang [120] applied the GREET model to characterize the soil C sequestration and GHG mitigation potential for the bioenergy crops switchgrass, poplars and willows. Model estimates indicated that there is a 70-85% reduction in GHGs emission when agricultural lands are converted to grasslands and forests.

DAYCENT is the daily time-step version of the CENTURY biogeochemical model [121]. DAYCENT simulates fluxes of C and N among the atmosphere, vegetation and soil [107]. Adler et al. [122] routinely applied the DAYCENT model to assess and compare GHG fluxes and biomass yields for corn, soybean, alfalfa, hybrid poplar, reed canarygrass, and switchgrass as bioenergy crops in Pennsylvania, US. All simulated cropping systems provided net GHG sinks compared with the fossil fuel life cycle, even in the long term when there were no further increases in soil C sequestration due to SOC storage reaching a new steady state. Switchgrass and hybrid poplar (Populus spp.) provided the largest net GHG sinks, (> 0.5 Mg C ha–1 yr–1) for biomass conversion to ethanol, and > 1.1 Mg C ha–1 yr–1 for biomass gasification for electricity generation. Compared with the LCA of fossil fuels, switchgrass-derived biofuels reduced GHG emissions by more than 115%. Similar results were reported by [123], who also applied the DAYCENT model to calculate the net GHG fluxes, when converting from cotton and unmanaged grasslands to a switchgrass system in the southern U.S. Long-term simula‐ tions predicted a decrease in GHG emissions (0.3–1.0 Mg C ha-1 yr-1) for conversions from cotton to switchgrass at N application rates of 0–135 kg N ha-1. Conversely, conversion from unman‐ aged grasses to switchgrass resulted in annual increases of net GHG emissions (0.1–0.4 Mg C ha-1 yr-1) for switchgrass at no and low (45 kg N ha-1) fertilization rates. Using the DAYCENT model, Davis et al. [124] estimated the effects of replacing corn ethanol with switchgrass and miscanthus on GHG emissions in Central US. Overall, there was a 29 to 473% reduction in GHG emissions even after accounting for emissions associated with indirect land-use change. Conversion from a high-input annual crop to a low-input perennial crop for biofuel production can thus transition the central US from a net source to a net sink for GHGs. A simulation study with the DAYCENT model showed that cultivation of energy cane on former pasture on Spodosol soils in the southeast US has the potential for high biomass yield and a positive GHG emission balance [35].

The USGS General Ensemble Biogeochemical Modeling System (GEMS) [125] was developed to integrate well-established ecosystem biogeochemical models, such as CENTURY, with various spatial databases for the simulation of biogeochemical cycles over large areas. GEMS accounts for the impacts of atmospheric and climatic changes (e.g., changes in precipitation and temperature, CO2 enrichment, N deposition), land use change (e.g., biofuel crops, tillage, forest logging, reforestation, urbanization, fire fuel treatment), natural disturbances and extreme events (land fire, insect outbreak, drought, flooding, hurricane) on ecosystem C balances and GHG emissions.

Even with all the uncertainties associated with LCAs, most suggest that WSGs and SRWCs have higher net energy balances and positive GHG emission balances than conventional crops [122, 126, 127]. Sarker [128] concludes that the use of computational models to quantify GHG fluxes in LCAs of bioenergy crops gives a more accurate representation of the system than using constant IPCC [129] emission factors. The USDA-ARS GRACEnet (Greenhouse gas Reduction through Agricultural Carbon Enhancement network: www.ars.usda.gov/research/ GRACEnet) research initiative compares various GHG mitigation strategies based on field studies in several agroecoregions of the US. It is hoped that results from this environmental research program will be important to inform public policy on valuing the benefits of bioen‐ ergy crop production, and in addition, provide the much needed data for model GHG emission calibration and testing.

#### **2.3. Carbon Credit Markets and Enabling Policies**

TEM uses spatially referenced information on climate, elevation, soils, vegetation and water availability as well as soil- and vegetation-specific parameters to make monthly estimates of C and N fluxes and pool sizes of terrestrial ecosystems. Conversion of agricultural lands to bioenergy crops consequently affects the ecosystem C balance. Qin et al. [116] applied the TEM model to evaluate the impacts of this change on the C balance, assuming several land use change scenarios from corn, soybean, and wheat to bioenergy crops of switchgrass and miscanthus. Overall, SOC increased significantly when land use changed from corn, soybean or wheat, to either switchgrass or miscanthus, reaching an average of ~100 Mg C ha-1 in the

Relatively simple statistical models can also be useful in estimating SOC losses associated with land conversion. Anderson–Teixeira, et al. [117] used reduced ANOVA and the forced- and free-intercept regression statistical models to estimate and compare SOC losses associated with land conversion and SOC sequestration rates under corn, sugarcane, switchgrass, miscanthus and mixed native grass. Their predicted SOC accumulation rates under switchgrass, miscan‐ thus and mixed native grass ranged from 0.1- 1.0 Mg C ha-1 yr-1 in the top 30 cm, while SOC losses under corn and sugarcane ranged from 3.0–8.0 Mg C ha-1 yr-1. These results are consistent with data from numerous field studies reported elsewhere in the literature [20, 28, 31, 40]. Life-cycle assessments are commonly used to determine whether a biomass biofuel system results in a net reduction in GHG emissions or an improved energy source. In LCAs, all input and output data in all phases of the product's life cycle including biomass production, feedstock storage, feedstock transportation, biofuel production, biofuel transportation and final use are required. Because field measurements of bioenergy LCA components are difficult to carry out, simulation models have been applied effectively to obtain realistic estimates of

Standard and useful tools for LCA are multidimensional spreadsheet models such as the GREET [118] and GHGenius [119] software models, which were specifically designed to address full LCA of biofuels. These spreadsheet models have advantages in that they are userfriendly, publicly available, straightforward, and relatively transparent. In a study commis‐ sioned by the US Department of Energy, Wang [120] applied the GREET model to characterize the soil C sequestration and GHG mitigation potential for the bioenergy crops switchgrass, poplars and willows. Model estimates indicated that there is a 70-85% reduction in GHGs

DAYCENT is the daily time-step version of the CENTURY biogeochemical model [121]. DAYCENT simulates fluxes of C and N among the atmosphere, vegetation and soil [107]. Adler et al. [122] routinely applied the DAYCENT model to assess and compare GHG fluxes and biomass yields for corn, soybean, alfalfa, hybrid poplar, reed canarygrass, and switchgrass as bioenergy crops in Pennsylvania, US. All simulated cropping systems provided net GHG sinks compared with the fossil fuel life cycle, even in the long term when there were no further increases in soil C sequestration due to SOC storage reaching a new steady state. Switchgrass and hybrid poplar (Populus spp.) provided the largest net GHG sinks, (> 0.5 Mg C ha–1 yr–1) for biomass conversion to ethanol, and > 1.1 Mg C ha–1 yr–1 for biomass gasification for electricity generation. Compared with the LCA of fossil fuels, switchgrass-derived biofuels

emission when agricultural lands are converted to grasslands and forests.

two bioenergy crops.

264 CO2 Sequestration and Valorization

the bioenergy crop's GWP or GHG mitigation potential.

It has been demonstrated conclusively that bioenergy crops can sequester soil C and mitigate GHG emissions. Accreditation of bioenergy production to CCMs will provide the much needed incentives for land managers and bioenergy feedstock producers to expand produc‐ tion, both for GHG mitigation benefits and the resulting improvements to soil health. Pio‐ neered by the European Union's (EU) Emissions Trading Schemes (ETSs) launched in 2005, several CCMs are now operating across the world - see [130]. For agriculture, the schemes or programs allow producers and landowners to earn income by storing C in their soil through GHG capture projects. Soil C credits are then traded on CCMs, much like any other agricultural commodity. Large companies and other entities purchase C credits on the CCMs to offset their own GHG emissions into the atmosphere.

Carbon credits are allocated to an enterprise or project based on verifiable GHG emission reductions. For bioenergy feedstock producers to participate in CCMs there is need for lowcost, simple, but yet reliable methods for verifying soil C change in bioenergy cropping systems under different soil types, climates and management practices. As discussed earlier, field monitoring of soil C to document soil C changes is not only cost prohibitive but is generally not feasible due to the spatial variability in C stocks, and for most systems, the small annual incremental changes are not detectable by present analytical methods. According to [131], quantifying soil C may be even more costly than C stored in above-ground biomass, consid‐ ering that soil C cannot be observed in the ways that biomass can. It is apparent that fully parameterized process-based models could be adopted or adapted to provide the much needed tools to track soil C sequestration and GHG emissions for the CCMs. As pointed out by [132], a major reason for developing soil C or agroecosystem models is to be able to use them to predict SOC changes under climate-soil management combinations lacking in soil C measurements.

Models such as Century, RothC, EPIC, SOCRATES, and more recently, C-Farm [92] can be used to assess and verify C sequestration rates and GHG mitigation potential for the CCMs. For most models, applicability is often limited because they require numerous data inputs. More user-friendly tools with minimum calibration requirements are needed. The Comet 2.0 tool (http://www.cometvr.colostate.edu/) is an easy to use Carbon Management Tool and Greenhouse Gas Accounting system for farm-ranch-orchard operations. Growers in the US can apply the Comet tool to obtain farm level mean estimates and uncertainty of C sequestra‐ tion and GHG emissions from annual crops, hay, pasture and range, perennial woody crops (orchards, vineyards), agroforestry practices, and fossil fuel usage. The tool links a large set of databases containing information on soils, climate and management practices to dynamically run the CENTURY model as well as empirical models for soil N2O emissions and CO2 from fuel usage for field operations. Similarly, C-FARM is a user-friendly cropping systems model that allows calculating rates of C storage for the soil profile on a layer-by-layer basis. For simplicity, the model is built around a robust single-pool C model and can provide estimates of short- and long-term on-farm C storage rates. According to [92], C-FARM is a useful tool for growers, consultants, and state and federal agencies that need to develop C storage estimates for local conditions.

The energy crisis of 1973, and more recent concerns about global climate change and energy security have led to significant policy support for bioenergy throughout the world. The United States, European Union (EU), Australia, Canada and Switzerland spent at least \$11 billion on biofuel subsidies in 2006 [133]. Here we briefly examine policies that could stimulate increased C sequestration by promoting bioenergy feedstocks production. Some suggest that, policy changes must focus on finding new avenues to divert public support to producers and landowners in the form of stewardship or green payments as well as making progress in encouraging the success of private C trading systems [28].

Most producers are reluctant to take their land out of conventional crop production into bioenergy crops for a long period of time. Governments should introduce policies that include programs that help farmers make the transition to growing bioenergy crops. Although only partially implemented, the USDA Biomass Crop Assistance Program (BCAP) is meant to provide cost-share payments for establishing dedicated energy crops and annual payments to cover the foregone income from alternative uses of the land [28]. Such payments give producers the opportunity to preserve their land by increasing the SOC pool, while at the same time producing a crop that helps to meet future energy needs.

It is however important to note that while subsidies do create incentives for a transition to biofuels, they may also encourage increased fuel consumption by lowering the price of blended fuel, thus mitigating the benefits of land conversion to biofuels through increased GHG emissions. Many environmental economists recognize that a tax or fee on CO2 emission from fossil fuel sources is the most efficient system to reduce emissions and spread the burden equitably across all sources: industrial and personal [134]. A CO2 tax is currently under discussion in a number of countries, including the US and EU. In general, the tax would be applied to all fossil fuels before combustion according to the C content of the respective fuel.

Although CCMs keep developing throughout the world, their sustained future existence is uncertain as they face many challenges: the global economic crisis and its impacts on trading markets; failure by governments to pass cap-and-trade legislation, especially the US, one of the world's largest GHG emitters; great uncertainty over future investments in some devel‐ oping countries as CCM investors do not know where demand will come from; and finally, the international community's inability to agree on a post-2012 Kyoto framework in Copen‐ hagen (COP15) and Cancun (COP16), which has seriously damaged the confidence of the private sector in the long-term viability of existing CCMs [135].

## **3. Concluding remarks**

programs allow producers and landowners to earn income by storing C in their soil through GHG capture projects. Soil C credits are then traded on CCMs, much like any other agricultural commodity. Large companies and other entities purchase C credits on the CCMs to offset their

Carbon credits are allocated to an enterprise or project based on verifiable GHG emission reductions. For bioenergy feedstock producers to participate in CCMs there is need for lowcost, simple, but yet reliable methods for verifying soil C change in bioenergy cropping systems under different soil types, climates and management practices. As discussed earlier, field monitoring of soil C to document soil C changes is not only cost prohibitive but is generally not feasible due to the spatial variability in C stocks, and for most systems, the small annual incremental changes are not detectable by present analytical methods. According to [131], quantifying soil C may be even more costly than C stored in above-ground biomass, consid‐ ering that soil C cannot be observed in the ways that biomass can. It is apparent that fully parameterized process-based models could be adopted or adapted to provide the much needed tools to track soil C sequestration and GHG emissions for the CCMs. As pointed out by [132], a major reason for developing soil C or agroecosystem models is to be able to use them to predict SOC changes under climate-soil management combinations lacking in soil C

Models such as Century, RothC, EPIC, SOCRATES, and more recently, C-Farm [92] can be used to assess and verify C sequestration rates and GHG mitigation potential for the CCMs. For most models, applicability is often limited because they require numerous data inputs. More user-friendly tools with minimum calibration requirements are needed. The Comet 2.0 tool (http://www.cometvr.colostate.edu/) is an easy to use Carbon Management Tool and Greenhouse Gas Accounting system for farm-ranch-orchard operations. Growers in the US can apply the Comet tool to obtain farm level mean estimates and uncertainty of C sequestra‐ tion and GHG emissions from annual crops, hay, pasture and range, perennial woody crops (orchards, vineyards), agroforestry practices, and fossil fuel usage. The tool links a large set of databases containing information on soils, climate and management practices to dynamically run the CENTURY model as well as empirical models for soil N2O emissions and CO2 from fuel usage for field operations. Similarly, C-FARM is a user-friendly cropping systems model that allows calculating rates of C storage for the soil profile on a layer-by-layer basis. For simplicity, the model is built around a robust single-pool C model and can provide estimates of short- and long-term on-farm C storage rates. According to [92], C-FARM is a useful tool for growers, consultants, and state and federal agencies that need to develop C storage

The energy crisis of 1973, and more recent concerns about global climate change and energy security have led to significant policy support for bioenergy throughout the world. The United States, European Union (EU), Australia, Canada and Switzerland spent at least \$11 billion on biofuel subsidies in 2006 [133]. Here we briefly examine policies that could stimulate increased C sequestration by promoting bioenergy feedstocks production. Some suggest that, policy changes must focus on finding new avenues to divert public support to producers and

own GHG emissions into the atmosphere.

266 CO2 Sequestration and Valorization

measurements.

estimates for local conditions.

The potential for bioenergy landscapes under WSGs and SRWCs to sequester soil C and mitigate GHG emissions is indisputable. The functions and benefits of SOC in the maintenance of soil quality, productivity and sustainability of cropping systems are well documented. It is imperative to monitor the SOC integrity of bioenergy crops through the establishment of stringent SOC measurement and quantification methodologies. Available information suggests the best approach would be to combine field measurements with process-based simulation modeling. While models can capture complex C dynamics as impacted by sitespecific factors or conditions, field measurements are crucial for the verification of model outputs. Besides promoting participation in CCMs by producers, accurate monitoring of SOC stocks will also allow assessment of whether a given bioenergy cropping system meets government-stipulated life cycle GHG emission reductions. Whether or not countries imple‐ ment policies that stimulate increased C sequestration through bioenergy cropping, the importance of enhancing the soil resource base should be increasingly highlighted.

## **Acknowledgements**

Preparation of this book chapter was supported by Texas A&M AgriLife Research and the USDA-ARS, Grassland, Soil and Water Research Laboratory, Temple, Texas, through Specific Cooperative Agreement: 58-6206-1-053, and partially funded by the U.S. Navy, Office of Naval Research. USDA is an equal opportunity provider and employer.

## **Author details**

Manyowa N. Meki1\*, James R. Kiniry2 , Kathrine D. Behrman2 , Meghan N. Pawlowski3 and Susan E. Crow3

\*Address all correspondence to: nmeki@brc.tamus.edu, jim.kiniry@ars.usda.gov

1 Texas A&M AgriLife Blackland Research and Extension Center, Temple, US.

2 USDA, Agricultural Research Service, Grassland Soil and Water Research Laboratory, Temple, US.

3 Dept. of Natural Resources and Environmental Management, University of Hawaii at Manoa, Honolulu, US.

## **References**


ment policies that stimulate increased C sequestration through bioenergy cropping, the

Preparation of this book chapter was supported by Texas A&M AgriLife Research and the USDA-ARS, Grassland, Soil and Water Research Laboratory, Temple, Texas, through Specific Cooperative Agreement: 58-6206-1-053, and partially funded by the U.S. Navy, Office of Naval

, Kathrine D. Behrman2

2 USDA, Agricultural Research Service, Grassland Soil and Water Research Laboratory,

3 Dept. of Natural Resources and Environmental Management, University of Hawaii at

[1] Meki NM., Kemanian AR., Potter SR., Blumenthal JM., Williams JR., Gerik, TJ. Crop‐ ping systems effects on sorghum grain and biomass yield, soil organic carbon and global warming potential in central and south Texas. Journal Agricultural Systems

[2] Miller SA., Landis AE., Theis TL. Environmental tradeoffs of biobased production.

[3] US EPA (Environmental Protection Agency). Renewable Fuel Standard Program: (RFS2) Regulatory Impact Analysis EPA-420-R-10-006, February 2010. http://

[4] Meki NM., Kiniry RJ. A Dynamic Tool. International Innovation: The Global Fore‐

www.epa.gov/otaq/renewablefuels/index.htm (accessed 15 August 2013).

cast, October 2013. Research Media, UK, p118-120, ISSN 2051-8544. 2013.

Environmental Science and Technology 2007;41:5176–5182.

\*Address all correspondence to: nmeki@brc.tamus.edu, jim.kiniry@ars.usda.gov

1 Texas A&M AgriLife Blackland Research and Extension Center, Temple, US.

, Meghan N. Pawlowski3

and

importance of enhancing the soil resource base should be increasingly highlighted.

Research. USDA is an equal opportunity provider and employer.

**Acknowledgements**

268 CO2 Sequestration and Valorization

**Author details**

Susan E. Crow3

Temple, US.

**References**

Manoa, Honolulu, US.

2013;117: 19-29.

Manyowa N. Meki1\*, James R. Kiniry2


house gas-altered climate in the Central United States: A Simulation study. Agriculture, Ecosystems and Environment 2000;78: 31–47.


[31] Ma Z., Wood CW., Bransby DJ. Soil management impacts on soil carbon sequestra‐ tion by switchgrass. Biomass Bioenergy 2000;18: 469–477.

house gas-altered climate in the Central United States: A Simulation study.

[19] Frank AB., Berdahl JD., Hanson JD., Liebig MA., Johnson HA. Biomass and carbon

[20] Liebig MA., Johnson HA., Hanson JD., Frank AB. Soil carbon under switchgrass

[21] Blanco-Canqui H. Energy crops and their implications on soil and environment.

[22] Kahle P., Boelcke B., Zacharias S. Effects of *Miscanthus x giganteus* cultivation on chemical and physical soil properties. Journal of Plant Nutrition and Soil Science

[23] Hansen EM., Christensen BT., Jensen LS., Kristensen K. Carbon sequestration in soil beneath long-term Miscanthus plantations as determined by C-13 abundance. Bio‐

[24] Felten D., Emmerling C. Accumulation of Miscanthus-derived carbon in soils in rela‐ tion to soil depth and duration of land use under commercial farming conditions.

[25] Bradley RI., King JA. 2004. A review of farm management techniques that have im‐ plications for carbon sequestration – validating an indicator. OECD Expert Meet. Farm Management Indicators and the Environment. 8-12 March, 2004, Palmerston North. http://webdomino1.oecd.org/comnet/agr/farmind.nsf/22afaeb‐

ba539ba74c1256a3b004d5175/b3b8d25f219f4ae3c1256bd5004874f1/\$FILE/Brad‐

[26] Borzêcka-Walker M., Faber A., Borek R. Evaluation of carbon sequestration in ener‐ getic crops (*Miscanthus* and coppice willow). International Agrophysics 2008;22:

[27] Don A., Osborne B., Hastings A. et al. Land-use change to bioenergy production in Europe: implications for the greenhouse gas balance and soil carbon. Global Change

[28] Lemus R., Lal R. Bioenergy crops and carbon sequestration. Critical Reviews in Plant

[29] Blanco-Canqui H., Lal R., Lemus R. Soil aggregate properties and organic carbon for switchgrass and traditional agricultural systems in the southeastern United States.

[30] Tufekcioglu A., Raich JW., Isenhart TM., Schultz RC. Biomass, carbon and nitrogen dynamics of multi-species riparian buffers within an agricultural watershed in Iowa,

Agriculture, Ecosystems and Environment 2000;78: 31–47.

Agronomy Journal 2010;102: 403-419.

mass and Bioenergy 2004;26: 97–105.

ley1.pdf (accessed 25 August 2013).

Biology Bioenergy 2012;4: 372–391.

Soil Science 2005;170: 998–1012.

US. Agroforestry Systems 2003;57: 187–198.

Sciences 2005;24: 1–21.

1999;162: 27–32.

270 CO2 Sequestration and Valorization

185-190.

partitioning in switchgrass. Crop Science 2004; 44: 1391–1396.

Journal of Plant Nutrition and Soil Science 2012;175: 661–670.

stands and cultivated cropland. Biomass Bioenergy 2005;28: 347–354.


[58] Saby NPA., Bellamy PH., Morvan X., Arrouays D., Jones RJA., Verheijen FGA., Kib‐ blewhite MG., Verdoot AY., Üveges JB., Freudenschuß A., Simota C. Will European soilmonitoring networks be able to detect changes in topsoil organic carbon? Global Change Biology 2008;14: 1-11

[45] Harmon ME., Ferrell W., Franlkin JF. Effects of carbon storage of conversion of old-

[46] Cannell MGR. Carbon sequestration and biomass energy offset: Theoretical, poten‐ tial, and achievable capacities globally, in Europe and the UK. Biomass and Bioener‐

[47] Bronson KF., Mosier AR. Nitrous oxide emissions and methane consumption in wheat and corn cropped systems. In: Harper LA., Mosier AR., Duxbury JM., Rolston DE. (eds.) Agriculture Ecosystem Effects on Trace Gases and Global Climate Change.

[48] Robertson GP., Paul EA., Harwood RR. Greenhouse gases in intensive agriculture: contributions of individual gases to radiative forcing of the atmosphere. Science

[49] Tuskan GA., Walsh ME. Short-rotation woody crop systems, atmospheric carbon di‐ oxide and carbon management: A US case study. Forestry Chronicle 2001;77: 259–64.

[50] Machefert SE., Dise NB., Goulding KT., Whitehead PG. Nitrous oxide emission from a range of land uses across Europe. Hydrological Earth System Science 2002;6: 325–

[51] Sartori F., Lal R., Ebinger, HE., Parrish DJ. Potential Soil Carbon Sequestration and CO2 Offset by Dedicated Energy Crops in the USA. Critical Reviews in Plant Scien‐

[52] Scheer C., Wassmann R., Kienzler K., Ibragimov N., et al. Methane and nitrous oxide fluxes in annual and perennial land-use systems of the irrigated areas in the Aral Sea

[53] Volk TA., Verwijst T., Tharakan PJ., Abrahamson LP., White EH. Growing fuel: a sustainability assessment of willow biomass crops. Frontiers in Ecology and the En‐

[54] Reeves DW. The role of soil organic matter in maintaining soil quality in continuous

[55] Karlen DL., Varvel GE., Johnson JM., Baker JM., Osborne SL., Novak JM., Adler PR., Roth GW., Birrell SJ. Monitoring soil quality to assess the sustainability of harvesting

[56] Marland G., Schlamadinger B. The Kyoto Protocol could make a difference for the optimal forest-based CO2 mitigation strategy: some results from GORCAM. Environ‐

[57] Nelson DW., Sommers LE. Total carbon, organic carbon, and organic matter. In: Sparks DL. (ed.) Methods of soil analysis. Part 3: Chemical methods. Madison, WI:

growth forests to young forests. Science 1990;247: 699–702.

ASA-CSSA-SSSA, Madison, WI. 1993. p133–144.

Basin. Global Change Biology 2008;14: 2454–2468.

cropping systems. Soil Tillage Research 1997;43: 131–167.

corn stover. Agronomy Journal 2011;103: 288–295.

mental Science and Policy 1999;2: 111-124.

Soil Science Society of America. 1996.

gy 2003; 24: 97–116.

272 CO2 Sequestration and Valorization

2000;289 : 1922–1925.

ces 2006;25: 441-472.

vironment 2004;2(8): 411- 418.

337.


soils.usda.gov/survey/geography/ssurgo/description.html. (accessed 9 September 2013).

[83] Sombroek WG., Nachtergaele FO., Hebel A. Amounts, dynamics and sequestration of carbon in tropical and subtropical soils. Ambio 1993;22: 417–426.

[71] Kibblewhite M., Jones RJA., Baritz R., Huber S., Arrouays D., Micheli E., Dufour MJD. ENVASSO. Environmental assessment of soil for monitoring. In: EC desertifi‐ cation meeting, Brussels. 2005. http://eusoils.jrc.ec.europa.eu/projects/envasso/ (ac‐

[72] Wesemael Bas van, Paustian K., Andrén O., Cerri EPC., et al. How can soil monitor‐ ing networks be used to improve predictions of organic carbon pool dynamics and

[73] EISA (Energy Independence and Security Act) http://www.afdc.energy.gov/laws/

[74] Ellert, BH., Janzen HH., McConkey B. A Method for measuring soil carbon change on the Canadian Prairies. In: Lal R., Kimble J., Follett R. (eds.) An International Workshop on Assessment Methods for Soil C Pools. Lewis Publishers. CRC Press,

[75] Izaurralde RC., McGill WB., Bryden A., Graham S., Ward M., Dickey P. Scientific challenges in developing a plan to predict and verify carbon storage in Canadian Prairie soils. In: Lal R., Kimble J., Follett R. (eds.) An International Workshop on As‐ sessment Methods for Soil C Pools. Lewis Publishers. CRC Press, Inc., Boca Raton,

[76] Kern JS., Turner DP., Dodson RF. Spatial patterns in soil organic carbon pool size in the Northwestern United States. In: Lal R., Follett R., Stewart BA. (eds.) Soil Process‐ es and the Carbon Cycle. Advanced Soil Science. CRC Press, Inc., Boca Raton, Fla.,

[77] Tarnocai C., Ballard M. Organic carbon in Canadian soils. In: Lal R., Kimble J., Levine E. (eds.) Soil Processes and Greenhouse Effect. USDA, Soil Conservation Service, Na‐

[78] Post WM., Pastor J., Zinke PJ., Stangenberger AG. Global soil nitrogen. Nature

[79] Eswaran H., Van den Berg E., Reich P., Kimble J. Global Soil Carbon Resources. In: Lal R., Kimble J., Levine E., Stewart BA. (eds.), Advances in Soil Science: Soils and

[80] Batjes NH. 1996. Total carbon and nitrogen in the soils of the world. European Jour‐

[81] NRCS (Natural Resources Conservation Service) 1994, State Soil Geographic (STATS‐ GO) Data Base - Data Use Information. US Department of Agriculture. http:// soils.usda.gov/survey/geography/ssurgo/description\_statsgo2.html. (accessed 9 Sep‐

[82] NRCS (Natural Resources Conservation Service). 1995 Soil Survey Geographic (SSURGO) Data Base – Data Use Information. US Department of Agriculture. http://

Global Change, Lewis Publishers, CRC Press. Boca Raton, FL. 1995. p. 27–43.

tional Soil Survey Center, Lincoln, NE, USA. 1994. p31–45.

CO2 fluxes in agricultural soils. Plant Soil 2011;338: 247–259.

cessed on 11 September 2013)

274 CO2 Sequestration and Valorization

Inc., Boca Raton, FL. 2000.

FL. 1998. p433–446.

1998. p29-44.

1985;317: 613–616.

tember 2013).

nal Soil Science 1996;47: 151–163.

eisa. 2007. (accessed 11 September 2013).


Powlson DS., Smith P., Smith JU. (eds.) NATO ASI Series I. Springer, Berlin; 1996. p237–246.


[108] Gelfand I., Sahajpal R., Zhang X., Izaurralde RC., Gross KL., Robertson GP. Sustaina‐ ble bioenergy production from marginal lands in the US Midwest. 2013 doi:10.1038/ nature11811.

Powlson DS., Smith P., Smith JU. (eds.) NATO ASI Series I. Springer, Berlin; 1996.

[95] Dungait AJ., Hopkins DW., Gregory AS., Whitmore AP. Soil organic matter turnover is governed by accessibility not recalcitrance. Global Change Biology 2012;18: 1781–

[96] Six J., Jastrow JD. Organic matter turnover. Encyclopedia of Soil Science. 2002.

[97] Parton WJ., Schimel DS., Cole CV., Ojima DS. Analysis of factors controlling soil or‐ ganic matter levels in Great Plains grasslands. Soil Science Society America Journal

[98] Parton WJ., Rasmussen PE. Long-term effects of crop management in wheat/fallow: II. CENTURY model simulations. Soil Science Society America Journal 1994;58:

[99] Johnson K., Scatena FN., Pan Y. Short- and long-term responses of total soil organic carbon to harvesting in a northern hardwood forest. Forest Ecology and Manage‐

[100] Izaurralde RC., Williams JR., McGill WB., Rosenberg NJ., Jakas MC. Simulating soil C dynamics with EPIC: Model description and testing against long-term data. Eco‐

[101] Li C., Frolking S., Harriss R. Modeling carbon biogeochemistry in agricultural soils.

[102] Grant RF., Juma NG., McGill WB. Simulation of carbon and nitrogen transformations in soil: mineralization. Soil Biology and Biochemistry 1993;25: 1317–1329.

[103] Grace PR., Ladd JN. SOCRATES v2.00 User Manual. Cooperative Research Centre for Soil and Land Management, PMB 2 Glen Osmond 5064, South Australia. 1995.

[104] Grace PR.,, Ladd JN., Robertson GP., Gage SH. SOCRATES - A simple model for pre‐ dicting long-term changes in soil organic carbon in terrestrial ecosystems. Soil Biolo‐

[105] Kelly RH., Parton WJ., Crocker GJ., Grace PR., Klir J., Korschens M., Poulton, PR., Richter DD. Simulating trends in soil organic carbon in long-term experiments using

[106] Stockle CO., Donatelli M., Nelson R. CropSyst, a cropping systems simulation model.

[107] Del Grosso SJ., Mosier AR., Parton WJ., Ojima DS. DAYCENT model analysis of past and contemporary soil N2O and net greenhouse gas flux for major crops in the US.

p237–246.

276 CO2 Sequestration and Valorization

1796.

p936-942.

530-536.

1987;51: 1173-1179.

ment 2010;259: 1262–1267.

logical Modeling 2006;192: 362-384.

Global Biogeochemical Cycles 1994;8: 237–254.

gy and Biochemistry 2006;38: 1172–1176.

the century model. Geoderma 1997;81: 75-90.

Soil Tillage Research 2005;83: 9-24.

European Journal of Agronomy 2003;18: 289–307.


[121] Parton WJ., Hartman MD., Ojima DS., Schimel DS. DAYCENT: Itsland surface sub‐ model: description and testing. Global Planetary Changes 1998;19: 35-48.

[122] Adler PR., Del Grosso SJ., Parton WJ. Life cycle assessment of net greenhouse gas flux for bioenergy cropping systems. Ecological Applications 2007;17: 675–691. [123] Chamberlain JF., Miller S., Frederick J. Using DAYCENT to quantify on-farm GHG emissions and N dynamics of land use conversion to N-managed switchgrass in the

Southern U.S. Agriculture, Ecosystems and Environment 2011;141: 332-341.

10.1890/110003.

278 CO2 Sequestration and Valorization

2008;96: 224-236.

1608-1619.

Press, Cambridge, UK. 1997.

Energy Policy 2011;39: 6040–6054.

[124] Davis SC., Parton WJ., Del Grosso SJ., Keough C., Marx E., AdlerPR., DeLucia EH. Impact of second-generation biofuel agriculture on greenhouse-gas emissions in the corn-growing regions of the US. Frontiers in Ecology and the Environment 2011. doi:

[125] Liu S., Tan Z., Chen M., Liu J., Wein A., Li Z., Huang J., Oeding J., Young CJ., Verma SB., Suyker AE., Faulkner S., McCarty GW. The General Ensemble Biogeochemical Modeling System (GEMS) and its applications to agricultural systems in the United States. In: Liebig MA., Franzluebbers AJ., Follett RF. (eds.) Managing agricultural greenhouse gases-coordinated agricultural research through GRACEnet to address

[126] Boehmel C., Lewandowsk I., Claupein W. Comparing annual and perennial energy cropping systems with different management intensities. Agricultural Systems

[127] Tilman D, Hill J., Lehman, C. Carbon-negative biofuels from low-input high-diversi‐

[128] Sarkar S., Miller SA., Frederick JR., Chamberlain JF. 2011. Modeling nitrogen loss from switchgrass agricultural systems. Biomass and Bioenergy 2011;35: 4381-4389.

[129] IPCC (Intergovernmental Panel on Climate Change). Greenhouse Gas Inventory Ref‐ erence Manual. Intergovernmental Panel on Climate Change. Cambridge University

[130] Perdan S. Azapagic A. Carbon trading: Current schemes and future developments.

[131] Antle JM., Capalbo SM., Mooney M., Elliott ET., Paustian KH. 2003. Spatial heteroge‐ neity, contract design, and the efficiency of carbon sequestration policies for agricul‐ ture. Journal of Environmental Economics and Management 2003;46(2): 231–50. [132] He X., Izaurralde RC., Vanotti MB., Williams JR., Thomson AM. Simulating longterm and residual effects of nitrogen fertilization on corn yields, soil carbon seques‐ tration, and soil nitrogen dynamics. Journal of Environmental Quality 2006;35:

our changing climate: San Diego, Calif., Academic Press. 2012. p309-323.

ty grassland biomass. Science 2006;314: 1598-600.

[135] IETA (International Emissions Trading Association):www.ieta.org. (accessed 15 Au‐ gust 2013).

## **Pre-Injection Phase: Site Selection and Characterization**

B. Llamas, M. Arribas, E. Hernandez and L.F. Mazadiego

Additional information is available at the end of the chapter

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

## **1. Introduction**

Depleted oil and gas fields, deep saline aquifers and non-mineable coal beds fulfil the special requirements for CO2 storage and have become global references to develop this industry [2]. Depleted oil or gas fields have been well researched due to the associated industrial and economic value of these emplacements. While deep saline aquifers are among the most promising emplacements, since theoretically these structures offer the highest value in terms of capacity [2],[3], the risk associated with the exploration of these potential sites is greater than that for already investigated depleted oil or gas fields. In such cases, it is likely that multicriteria algorithms can facilitate evaluation to find the best option under consideration, so the results of this process will help the decision-maker to decrease the risk associated with the exploration of these potential emplacements.

The site selection phase comprises the identification, characterization and selection of em‐ placements that could be suitable for CCS among a list of candidates [4]. This phase is generally completed by the definition of the qualification criteria and the provision of the evidence concerning the reliable functioning of the emplacement according to these criteria. The selection of a suitable site also depends on the scale of the assessment. In this regard, every scale will be related to a different resolution and detail of information. It is possible to differentiate at least three levels, namely, basin level (identifying and quantifying large potential storage areas), regional level (increased level of detail, identifying areas to prospect) and local level (very detailed structures, pre-engineering site selection).

While the main mechanisms for CO2 storage have been identified, and a series of criteria for site assessment and selection have been developed [5],[6],[7], a multicriteria algorithm to quantify and rank the potential areas under consideration has, to our knowledge, never been applied in an absolute mode, meaning that any alternative will be compared against a pattern.

© 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.

Thus, the objective of this study was to develop a methodology based on multicriteria algorithms for assessing the best emplacement for CCS within a range of alternatives.

In general, most of the areas which could be suitable for storing CO2 are not well explored geologically. As a result, further exploration of the subsurface must be carried out, which implies higher cost and risks. In order to reduce the risk of failure, it is necessary to define a previous phase, which should be based on (1) data collection, and (2) the definition of a criterion and multicriteria decision tool.

The collection of data is independent of each country or region. It involves the study and processing of existing data, such as from geophysical surveys along with existing wells (shallow and deep) developed for oil and gas exploration, water resources and exploration, mining activities (exploration and extraction), exploration for nuclear waste deposits and underground natural gas storage activities. The data from these industrial activities can be complemented by other academic evaluations, including postgraduate theses and/or peer reviewed scientific articles.

**Figure 1.** Work flow proposed for basin screening (Definition phase)

Figure 1 represents a proposed work flow [4]. The screening phase could be differentiated by the Data Recompilation task and the Multicriteria Decision Tool. It is integrated as a prelimi‐ nary phase, and it is connected with a second phase called *characterization phase,* which corresponds to site maturation and testing.

## **2. Site criteria definition**

Thus, the objective of this study was to develop a methodology based on multicriteria

In general, most of the areas which could be suitable for storing CO2 are not well explored geologically. As a result, further exploration of the subsurface must be carried out, which implies higher cost and risks. In order to reduce the risk of failure, it is necessary to define a previous phase, which should be based on (1) data collection, and (2) the definition of a

The collection of data is independent of each country or region. It involves the study and processing of existing data, such as from geophysical surveys along with existing wells (shallow and deep) developed for oil and gas exploration, water resources and exploration, mining activities (exploration and extraction), exploration for nuclear waste deposits and underground natural gas storage activities. The data from these industrial activities can be complemented by other academic evaluations, including postgraduate theses and/or peer

Figure 1 represents a proposed work flow [4]. The screening phase could be differentiated by the Data Recompilation task and the Multicriteria Decision Tool. It is integrated as a prelimi‐ nary phase, and it is connected with a second phase called *characterization phase,* which

algorithms for assessing the best emplacement for CCS within a range of alternatives.

criterion and multicriteria decision tool.

**Figure 1.** Work flow proposed for basin screening (Definition phase)

corresponds to site maturation and testing.

reviewed scientific articles.

282 CO2 Sequestration and Valorization

It is not only necessary to evaluate specific sites with a technical point of view. Sometimes problems can involve economic aspects or social acceptance, which could make a CO2 storage site much more difficult (and costly). Therefore, it is necessary to establish a high level differentiation between technical and socioeconomic criteria (later we call this differentiation the 1st criteria level of the Analytical Hierarchy Process).

There is no standardization in this aspect, so the selection of the criteria should be as careful as possible, and should include all aspects that can make an area suitable or not. The criteria proposed in this chapter are based on the direct experience of the research group involved in this publication and on the evaluation of several publications and projects that focus on site selection methodology [10],[4].

Even though the criteria may be as described in the next sub-chapters, it is convenient to consider another type of classification, which is based on time scale and the possibility of modifying each criterion along that scale. Moreover, we can describe:


#### **2.1. Technical criteria**

These criteria relate to scientific aspects or parameters to provide confidence in the findings about the subsurface structure. Deep subsurface exploration implies higher risk because the exploration techniques available are expensive and the probability of success is not high. In order to reduce the technical risk, it is necessary to define those criteria relevant for considering every critical issue.

**Tectonics.** This parameter comprises aspects like the structural definition of the trap at a basin and local level. The former considers whether the sedimentary basin is convergent or diver‐ gent, as well as its neotectonic activities. The latter relates to the type of the structure and trap, whether it is an anticline, syncline or seal fault [10]. Furthermore, the geo-mechanic evaluation should be taken into account in order to evaluate the maximum CO2 injection pressure in the storage formation. This injectivity criterion defines the maximum capacity per unit of time.

**Geology.** This criterion evaluates the storage and the cap-rock formations. The common formations considered suitable for storing CO2 are sandstones and limestones. These geolog‐ ical formations tend to have high porosity. While sandstone has a primary porosity, which is much more homogeneous, limestone presents secondary porosity, which is created by processes of diagenesis (e.g., leeching of minerals or the generation of a fracture system). Other parameters to be defined in storage formations are permeability, and thickness of the forma‐ tion. The plasticity, porosity and thickness of the cap-rock formations should also be evaluated and considered. In relation to the plasticity of the cap-rock formation, it is more desirable to have a ductile rather than a brittle behavior.

**Hydrogeology.** This parameter should consider both the dynamics and the fluid quality in the reservoir formation. Hydrodynamic criteria describe the natural dynamic flow system and hence the potential for hydrodynamic trapping within the basin under assessment. Shallow, short flow systems therefore do not meet the geological requirements for maintaining super‐ critical CO2, in terms of depth, pressure and temperature, and do not have sufficient residence time to immobilize the injected CO2 by one of the other trapping mechanisms, such as residual trapping, solution trapping or mineral trapping. Flow rate is controlled by the driving forces of the fluid, including its buoyancy and hydraulic gradient, and by the permeability and porosity characteristics of the reservoir rock through which the fluid is moving.

An additional sub-criterion that should be considered is the quality of the fluid under consideration. In this regard, it is mandatory to consider the principle of sustainability, whereby present actions cannot compromise future generations. Fluid quality is measured in Total Dissolved Solids (TDS), which define drinkable water (TDS < 3000 ppm). Thus, when drinkable water is present, the aquifer is not suitable for CO2 storage. On the other hand, TDS higher than 10,000 ppm should be considered suitable for CO2 storage operations, as the quality of this water makes it not suitable for any other activity.

**Subsurface conditions of the CO2.** The most efficient way to store CO2 underground is to store it under supercritical conditions, [2]. This special state of the CO2 provides similar densities to a liquid (i.e., increasing the capacity per volumetric unit), while the viscosity is similar to a gas (i.e., increasing the capacity per unit of time). Supercritical conditions can be reached in CO2 geology storage when the depth of the storage formation is higher than 800 m and the geothermal gradient is low (below 25 ◦C/km) [5,[6].

**Capacity.** The current methods for estimating CO2 storage potential and capacity are based on widely accepted assumptions about geological trapping mechanisms, storage media and operating timeframes reviewed previously by other authors [2][6]. Using the concept of resources and reserves, the CSLF Task Force proposed a Techno-Economic Resource-Reserve Pyramid for CO2 Storage Capacity, [11]. The various capacities are nested within the resourcereserves pyramid, and defined as Theoretical, Effective, Practical and Matched Storage Capacity.

This parameter is considered a technical criterion, but it is also relevant for economic reasons: it is necessary to define the CO2 emitter in order to relate CO2 emissions and capacity.

**Other geological formations.** CO2 storage implies a long time period. Considering that the whole life cycle of the project includes characterization (pre-injection), injection, closure and post-closure as the most general phases, in addition to the storage itself, the project can last up to 100 years. Moreover, if there are any future modifications to the original conditions, such as mining activities or the development of oil fields, the confinement mechanisms may not be guaranteed. For this reason, it is necessary to consider every shallow geological formation that could be of potential use in the future. This will fulfill the sustainability principle previously described.

#### **2.2. Socioeconomic criteria**

parameters to be defined in storage formations are permeability, and thickness of the forma‐ tion. The plasticity, porosity and thickness of the cap-rock formations should also be evaluated and considered. In relation to the plasticity of the cap-rock formation, it is more desirable to

**Hydrogeology.** This parameter should consider both the dynamics and the fluid quality in the reservoir formation. Hydrodynamic criteria describe the natural dynamic flow system and hence the potential for hydrodynamic trapping within the basin under assessment. Shallow, short flow systems therefore do not meet the geological requirements for maintaining super‐ critical CO2, in terms of depth, pressure and temperature, and do not have sufficient residence time to immobilize the injected CO2 by one of the other trapping mechanisms, such as residual trapping, solution trapping or mineral trapping. Flow rate is controlled by the driving forces of the fluid, including its buoyancy and hydraulic gradient, and by the permeability and

An additional sub-criterion that should be considered is the quality of the fluid under consideration. In this regard, it is mandatory to consider the principle of sustainability, whereby present actions cannot compromise future generations. Fluid quality is measured in Total Dissolved Solids (TDS), which define drinkable water (TDS < 3000 ppm). Thus, when drinkable water is present, the aquifer is not suitable for CO2 storage. On the other hand, TDS higher than 10,000 ppm should be considered suitable for CO2 storage operations, as the quality

**Subsurface conditions of the CO2.** The most efficient way to store CO2 underground is to store it under supercritical conditions, [2]. This special state of the CO2 provides similar densities to a liquid (i.e., increasing the capacity per volumetric unit), while the viscosity is similar to a gas (i.e., increasing the capacity per unit of time). Supercritical conditions can be reached in CO2 geology storage when the depth of the storage formation is higher than 800 m and the

**Capacity.** The current methods for estimating CO2 storage potential and capacity are based on widely accepted assumptions about geological trapping mechanisms, storage media and operating timeframes reviewed previously by other authors [2][6]. Using the concept of resources and reserves, the CSLF Task Force proposed a Techno-Economic Resource-Reserve Pyramid for CO2 Storage Capacity, [11]. The various capacities are nested within the resourcereserves pyramid, and defined as Theoretical, Effective, Practical and Matched Storage

This parameter is considered a technical criterion, but it is also relevant for economic reasons:

**Other geological formations.** CO2 storage implies a long time period. Considering that the whole life cycle of the project includes characterization (pre-injection), injection, closure and post-closure as the most general phases, in addition to the storage itself, the project can last up to 100 years. Moreover, if there are any future modifications to the original conditions, such as mining activities or the development of oil fields, the confinement mechanisms may not be guaranteed. For this reason, it is necessary to consider every shallow geological formation that

it is necessary to define the CO2 emitter in order to relate CO2 emissions and capacity.

porosity characteristics of the reservoir rock through which the fluid is moving.

of this water makes it not suitable for any other activity.

geothermal gradient is low (below 25 ◦C/km) [5,[6].

Capacity.

have a ductile rather than a brittle behavior.

284 CO2 Sequestration and Valorization

These criteria include both economic aspects and parameters related to the social acceptance of the emplacement and its activity.

**The quantity and quality of the geological data,** although considered in many cases as an economic aspect, is one of the main criteria, as this information is used to determine and quantify the technical criteria referred to in the above text.

The more information is available, the fewer characterization methods will have to be applied. In this case, the characterization program (geophysics and wells) will be less expensive, and the risk of failure of the geological exploration will be reduced.

**CO2 sources.** This criterion relates to the location details of the major stationary sources, and the distance between them and the areas of interest. An additional parameter to be considered is the flue gas quality (CO2 quality). Certain gases (NOX or SOX) can increase the interaction between the storage formation and the fluid injected [12]. These chemical interactions between injected gases and the storage or caprock formation materials can create precipitates or dissolutions, which can modify storage specifications of the area or structure.

**Regional location.** Whether an area or structure is onshore or offshore has an important economic consideration, since generally it is likely to be cheaper and technically easier to implement a CO2 injection site onshore rather than offshore. On the other hand, public perception and land use issues may dictate that offshore sites are preferential for many CO2 storage projects.

**Maturity of the Area.** This criterion considers those aspects which define, at a local level, the location of the areas under consideration. It includes key aspects such as the climate, existing infrastructure that may be affected by the geological exploration or CO2 storage, and infra‐ structure that is required to develop the exploration program and engineering activities to develop the emplacement.

**Areas of interest: population, environmental and cultural resources.** This criterion refers to aspects that can affect acceptance of the emplacement by the community. Protest against the storage activity will be very demanding and obtaining legal approvals and permits may be delayed. For this reason, some special areas including major cities and those areas protected by the Natura 2000 network should not be considered as optimal areas for CO2 storage. Furthermore, other areas with relevant monuments or archaeological sites should be taken into account in the model for assessing a potential emplacement.

## **3. Multicriteria decision making (MCDM)**

Multiple criteria decision making (MCDM) is a methodology developed for making decisions in the presence of multiple, usually conflicting, criteria. Evaluation methods and multicriteria decisions include the selection of a set of feasible alternatives, the simultaneous optimization of several objective functions, and a decision-making process and evaluation procedures that must be rational and consistent. The application of a mathematical model of decision-making will help to find the best solution, establishing the mechanisms to facilitate the management of information generated by the various disciplines of knowledge.

Those problems in which decision alternatives are finite are called Discrete Multicriteria Decision problems. Such problems are most common in reality and this case scenario will be applied in solving the problem of *site selection for storing CO2*. Discrete MCDM is used to assess and decide on issues that by nature or design support a finite number of alternative solutions. Recently, Multicriteria Decision Analysis has been applied to hierarchy policy incentives for CCS [15] or to assess the role of CCS [16].

Assessment methods and criteria decision include selection among a set of feasible alterna‐ tives, optimization with various objective functions simultaneously, a decision-maker and rational and consistent procedures for assessment. Its principles are derived from matrix theory, graph theory, organizational theory, measurement theory, theory of collective deci‐ sions, operations research and economics.

The main evaluation methods are: linear weighting (scoring), multi-attribute utility (MAUT), overcoming relationships and hierarchical analysis (AHP).

Some of the advantages of AHP over other methods of Multicriteria Decision are:


#### **3.1. Analytical Hierarchy Process (AHP)**

AHP is one of the most extensively used and powerful MCDM. Nowadays it is used by many companies in solving various multicriteria problems, ranking these in the following categories: selection, prioritization and assessment, provision of resources against a standard assessment, management and quality management and strategic planning. For example, AHP has been applied in the analysis of location, resource allocation, outsourcing, evaluation, manufactur‐ ing, marketing, supplier selection, finance, energy, education and risk analysis, [17]. This widespread use shows the suitability of AHP in solving various types of business decisionmaking problems.

The AHP overcomes the problems with a scoring approach by structuring complexity as a hierarchy and by deriving ratio scale measures through pairwise relative comparisons. Pairwise comparisons are basic to the AHP method. Hence, when comparing a pair of criteria, sub-criteria or alternatives, a ratio of relative importance can be established. The pairwise comparison process can be performed using words, numbers, or graphical bars.

**Figure 2.** AHP Components: Four steps to build a hierarchy or network structure.

Once the model is built, pairwise comparisons are made with all individual elements (criteria, sub-criteria and alternatives). This process allows giving numerical values to the judgments provided by people, which is also able to measure how each element contributes to each level of the hierarchy. Furthermore, the process is based on a well-defined structure consisting of arrays, and the ability of the eigenvalues to generate values or to approximate weights of each criterion [18], [19], [20]. The problem of finding a nonzero solution to this set of equations is very common in engineering and physics and is known as an *eigenvalue problem.*

In order to carry out these comparisons, the AHP uses a fundamental scale of numbers that have proven absolute in practice and that have been experimentally validated for physical problems and decisions. This scale assigns mathematical values with respect to quantitative or qualitative attributes, homogenizing each valuable criterion.

Figure 4 illustrates the process followed for every criterion. As an example, Original Fluid Quality should be evaluated considering the Water Quality for different uses (agricultural sector, human use) and the European Directive for CO2 storage; it is possible to establish different mathematical values for each measurable criterion.

#### **3.2. Construction of the decision tree**

decisions include the selection of a set of feasible alternatives, the simultaneous optimization of several objective functions, and a decision-making process and evaluation procedures that must be rational and consistent. The application of a mathematical model of decision-making will help to find the best solution, establishing the mechanisms to facilitate the management

Those problems in which decision alternatives are finite are called Discrete Multicriteria Decision problems. Such problems are most common in reality and this case scenario will be applied in solving the problem of *site selection for storing CO2*. Discrete MCDM is used to assess and decide on issues that by nature or design support a finite number of alternative solutions. Recently, Multicriteria Decision Analysis has been applied to hierarchy policy incentives for

Assessment methods and criteria decision include selection among a set of feasible alterna‐ tives, optimization with various objective functions simultaneously, a decision-maker and rational and consistent procedures for assessment. Its principles are derived from matrix theory, graph theory, organizational theory, measurement theory, theory of collective deci‐

The main evaluation methods are: linear weighting (scoring), multi-attribute utility (MAUT),

Some of the advantages of AHP over other methods of Multicriteria Decision are:

**3.** It allows measuring quantitative and qualitative criteria using a common scale.

**5.** It enables checking the consistency index and making corrections, if applicable.

**6.** It generates a synthesis and provides the ability to perform sensitivity analyses.

**4.** It includes participation of different people or groups of interest to build consensus.

**7.** It is easy to use and allows the solution to be complemented with mathematical optimi‐

AHP is one of the most extensively used and powerful MCDM. Nowadays it is used by many companies in solving various multicriteria problems, ranking these in the following categories: selection, prioritization and assessment, provision of resources against a standard assessment, management and quality management and strategic planning. For example, AHP has been applied in the analysis of location, resource allocation, outsourcing, evaluation, manufactur‐ ing, marketing, supplier selection, finance, energy, education and risk analysis, [17]. This widespread use shows the suitability of AHP in solving various types of business decision-

of information generated by the various disciplines of knowledge.

CCS [15] or to assess the role of CCS [16].

286 CO2 Sequestration and Valorization

sions, operations research and economics.

**1.** It has a mathematical basis.

zation methods.

making problems.

**3.1. Analytical Hierarchy Process (AHP)**

overcoming relationships and hierarchical analysis (AHP).

**2.** It enables breaking down and analyzing a problem in parts.

As a major conclusion, a decision tree has been proposed. This model considers all the criteria described above, and they have been classified.

Weight assessment has been defined considering the AHP method: each level of the criteria and sub-criteria has been compared, using a comparison matrix, which should be constructed considering the consistency principle (it should fulfill the transitivity and reciprocity rules).

**Figure 3.** Scientific scale and translation into an AHP scale. This translation facilitates the evaluation of individual crite‐ ria (quantitative or qualitative) using a homogeneous scale.

**Figure 4.** Analytical Hierarchy Process: criterion tree proposed to identify the best site for storing CO2.

In order to "recover" or find the vector of weights, [w1, w2, w3,..., wn] given to these ratios, the matrix product of matrix A with the vector w can be calculated and considered in an equation, which is described as the eigenvalue matrix equation. The problem of obtaining a nonzero solution to this set of equations is very common in engineering and physics and is known as an eigenvalue problem.

## *A* · *W* =*n* · *W*

Where *[A]* is the pairwise comparison matrix – where *n* Is the dimension – and *[W]* is the weight matrix (eigenvalues) for every criterion.

Site (Sn) assessment is evaluated using the formula:

$$\mathbf{S}\_{\boldsymbol{n}} = \sum\_{i=1}^{i=n} \mathbf{W}\_i \cdot \mathbf{V}\_i$$

In order to "recover" or find the vector of weights, [w1, w2, w3,..., wn] given to these ratios, the matrix product of matrix A with the vector w can be calculated and considered in an equation, which is described as the eigenvalue matrix equation. The problem of obtaining a

**Figure 4.** Analytical Hierarchy Process: criterion tree proposed to identify the best site for storing CO2.

**Figure 3.** Scientific scale and translation into an AHP scale. This translation facilitates the evaluation of individual crite‐

ria (quantitative or qualitative) using a homogeneous scale.

288 CO2 Sequestration and Valorization

Where *Wi* are the weights of each criterion, and *Vi* are the values assigned for the specific conditions of each site.

#### AHP as an absolute mode

AHP is a multicriteria methodology which has been developed for use in two different ways: relative and absolute mode. In the first case, all the alternatives are compared between each other, but no more than seven alternatives are recommended for evaluation at the same time.

There are two reasons to justify this limitation: (1) consistency principle and (2) neurons. Pairwise comparisons errors increase due to inconsistent judgments. It is possible to distribute this inconsistency among all the alternatives under evaluation. If the number of alternatives/ elements is low, the priorities will be less affected by this inconsistency. The neuronal explan‐ ation has its limits in the brain's ability to identify simultaneous events: the more criteria exist for pairwise comparison, the greater the risk of inconsistent judgments will be.

For this study, we consider the AHP algorithm in absolute mode. It requires a standard with which to compare alternatives. The process leads to absolute preservation in the rank of the alternatives no matter how many are introduced. In this case, it is possible to define a standard considering the best values for each criterion (see Table 1and Table 2).


**Table 1.** Technical criteria and their best values for CO2 storage alternatives.


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

## **4. CO2SITE ASSESS: Informatics tool**

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 which allows representing the results and *CO2 site storage assessment*.

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 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: bernardo.llamas@upm.es

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. The horizontal extent of the plume given by TOUGH2 simulations is

 100.75 m 95.58 m 0.054090814 140.49 m 135.17 m 0.039357846 168.37 m 165.55 m 0.017034129 191.34 m 191.16 m 0.00094162 211.27 m 213.72 m -0.011463597 229.23 m 234.12 m -0.020886725 235.34 m 252.88 m -0.069360962 260.81 m 270.34 m -0.035251905 275.70 m 286.74 m -0.038501779 289.36 m 302.25 m -0.042646816

**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 injection lasts for 25 years before it is shut down, and the total simulation time is 50 years. 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 injection lasts for 25 years before it is shut down, and the total simulation time is 50 years. 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. Simultaneously, CO2 concentration builds up beneath the shale layer and finally breaks through the capillary barrier upon sustaining sufficient CO2 column height. 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. Simultaneously, CO2 concentration builds up beneath the shale layer and finally breaks through the capillary barrier upon sustaining sufficient CO2 column height. 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/

[2] Pruess K. TOUGH2: A General Numerical Simulator for Multiphase Fluid and Heat Flow, LawrenceBerkeley Laboratory Report LBL-29400, California: Berkeley; 1999. [3] Pruess K., Oldenburg C., and Moridis G. TOUGH2 User's Guide, Version 2.0 (re‐ vised), Lawrence Berkeley Laboratory Report LBL-43134, California, Berkeley; 2011.

[4] Celia M.A. and Nordbotten J.M. Practical Modeling Approaches for Geological Stor‐

[5] Celia M.A. and Nordbotten J.M. How Simple Can We Make Models for CO2 Injec‐

[6] PetraSim User Manual, Thunderhead Engineering website, http://www.thunderhea‐ deng.com/downloads/petrasim/PetraSim-4-manual.pdf(accessed 24 May 2013).

[7] PetraSim Example Manual, Thunderhead Engineering website, http://www.thunder‐ headeng.com/downloads/petrasim/4/examples/PetraSimExamples.pdf(accessed 24

[8] Workshop on Numerical Models for Carbon Dioxide Storage in Geological Forma‐ tions, University of Stuttgart website, http://www.iws.uni-stuttgart.de/CO2-work‐

[9] Ebigbo A., Nordbotten J.M., and Class H. Numerical Investigations of CO2 Sequestra‐ tion in Geological Formations: Problem-Oriented Benchmarks, Problem 1, CO2Plume Evolution and Leakage through an Abandoned Well. R&D-Program of GEOTECH‐ NOLOGIEN website, http://www.hydrosys.uni-stuttgart.de/co2-workshop/Prob‐

[10] Ebigbo A., Nordbotten J.M., and Class H. Numerical Investigations of CO2 Sequestra‐ tion in Geological Formations: Problem-Oriented Benchmarks, Problem 2, Enhanced CH4 Recovery in Combination with CO2Storage in Depleted Gas Reservoirs. R&D-Program of GEOTECHNOLOGIEN website, http://www.hydrosys.uni-stuttgart.de/

[11] Ebigbo A., Nordbotten J.M., Class H. and Eigestad G. Numerical Investigations of CO2 Sequestration in Geological Formations: Problem-Oriented Benchmarks, Prob‐ lem 3, Estimation of the CO2 Storage Capacity of a Geological Formation. R&D-Pro‐ gram of GEOTECHNOLOGIEN website, http://www.hydrosys.uni-stuttgart.de/co2-

age of Carbon Dioxide. Underground Water 2009; 47(5) 627-638.

tion, Migration, and Leakage? Energy Procedia 2001; 4 3857-3864.


**Provisional chapter**

## **Seismic Reflectivity in Carbon Dioxide Accumulations: A Review Seismic Reflectivity in Carbon Dioxide Accumulations: a Review**

Claudia L. Ravazzoli and Julián L. Gómez Claudia L. Ravazzoli\* and Julián L. Gómez

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

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

**Introduction**

[22] Doughty C., Freifeld B.M., and Trautz R.C. Site Characterization for CO2 Geologic Storage and Vice Versa: The Frio Brine Pilot, Texas, USA as A Case Study. Environ‐

[23] Bachu S. and Bennion B. Effects of In-situ Conditions on Relative Permeability Char‐ acteristics of CO2-Brine Systems. Environmental Geology 2007; 54 1635–1656.

[24] Chadwick R.A. and Noy D.J. History- Matching Flow Simulations and Time-Lapse Seismic Data from the Sleipner CO2 Plume. In: Proceedings of the 7th Petroleum Geol‐

[25] Audigane P., Gaus I., Czernichowski-Lauriol I., Pruess K., and Xu T. Two-Dimen‐ sional Reactive Transport Modeling of CO2 Injection in a Saline Aquifer at the Sleip‐

[26] Singh V., Cavanagh A., Hansen H., Nazarian B., Iding M., and Ringrose P. Reservoir Modeling of CO2 Plume Behavior Calibrated Against Monitoring Data from Sleipner,

[27] Zhu C. CO2-Water-Rock Interaction in Geological Carbon Sequestration, Seminar

[28] Zhu C. and Lu P. Personal Communication, Department of Geological Sciences, Uni‐

mental Geology 2008; 54 1635–1656.

342 CO2 Sequestration and Valorization

ogy Conference, 2010; 7 1171–1182

versity of Indiana, 2012.

ner Site. American Journal of Science 2007; 307 974-1008.

Norway. Society of Petroleum Engineers 2010; 134891-PP.

presentation at Washington University in St. Louis, 2011.

It is widely recognized that the continuous increase in the concentration of carbon dioxide (CO2) in our atmosphere is a major cause of global climate change [1]. Consequently, in accordance with the objectives of the Kyoto Protocol, many carbon dioxide capture and storage projects are being developed worldwide to reduce and stabilize the emission of this greenhouse gas into the atmosphere. The geological storage of CO2 in many cases becomes a feasible option to accomplish this goal, giving rise to the science of *carbon sequestration*, a challenging task for governments, scientists and engineers [2].

The main targets for geological disposal of CO2 are depleted hydrocarbon reservoirs and deep saline aquifers. While the latter are more numerous, their characterization requires detailed studies because they are typically not explored for prospecting. The geological sequestration of carbon dioxide requires careful prior study and subsequent monitoring to prevent this gas from leaking to the surface. In general, the significant contrast between the physical properties of natural reservoir fluids and those of carbon dioxide allows the utilization of time lapse seismic methods to monitor the evolution of the injected CO2 volume.

While it is accepted that time lapse surveys are able to monitor the presence or absence of CO2 within a geological formation, their ability to quantify the saturation, state and volume of this fluid within the reservoir still needs research efforts from the geophysical community. This makes it necessary to search for reliable seismic indicators. Theoretical and numerical modeling provides, in this sense, a tool to explore useful correlations between the seismic response and relevant parameters of the CO2 repository.

The significance of the *seismic reflectivity* has long been recognized by many authors, resulting in the publication of numerous works for different constitutive models and applications. In close connection with this, the analysis of seismic amplitude variation with angle, hereafter called *AVA* is frequently used in reservoir geophysics to obtain information about the rocks and pore fluids. The behavior of the reservoir reflectivity for different overall CO2 saturations, distribution types, layer thicknesses, frequencies and CO2 physical states

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. © 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.

©2012 Ravazzoli and Gøsmez, licensee InTech. This is an open access chapter distributed under the terms of

controls the amplitude of seismic wave reflections and strongly conditions the detectability of the injected CO2 volume, as reported by [3],[4] and [5] among others.

With this idea, using rock physics models, in this chapter we focus on the theoretical analysis of the properties and variations of the seismic reflectivity and related parameters in a porous rock partially saturated with CO2 and brine. In our tests we model the seismic effects of different CO2 accumulations formed below impermeable seals.

We calibrate our models using information on the Utsira formation, a shallow saline aquifer in the Sleipner field (offshore Norway). Millions of tonnes of CO2 have been separated since October 1996 from natural gas and re-injected into the aquifer, consisting of a highly porous and unconsolidated sandstone, with several thin intra-reservoir shale layers. These shale intervals act as temporary seals causing accumulations of CO2 beneath them [6],[7], whose saturations increase with time. A 50-100 m thick high velocity layer overlies the sand, forming the reservoir caprock, known as the Nørland formation. It is mainly composed of shales, siltstones and mudstones [8]. Another high-velocity shale layer of about 5 m thick is located about 10 m below the top of the sand.

As pointed out by [9] the *topmost* CO2 layer is of special interest for monitoring the Sleipner injection site, since its growth reveals the total upward flux and its changes over time. Thus, the analysis of the reflectivity associated with this kind of CO2 accumulation is an important topic in CO2 sequestration problems. The present book chapter focuses on this subject and is structured in three main sections considering different aspects of the problem.

First, in Section 1 we present a description of the models and rock physics tools used in the subsequent sections. Next, in Section 2 we model the compressional P-wave reflection coefficient at the interface between a caprock and a permeable porous layer partially saturated by a mixture of brine and CO2 under liquid, supercritical and gaseous conditions. For this analysis we consider a simple model consisting of two halfspaces, combining a fluid substitution procedure and a Gassmann-Hill formulation to take into account variable spatial distributions of the fluids. We perform a sensitivity analysis of the standard AVA coefficients (intercept, gradient and curvature) in the near offset range, to investigate their correlations with CO2 saturation and its thermodynamic state within the geologic formation. Also, we analyze the effect of modeling the CO2 by means of simple and complex equations of state.

The two-halfspace representation considered in Section 2 is valid when the accumulation thicknesses are larger than the involved seismic wavelengths. When this assumption is not valid, the interference between the multiple waves reflected at the different boundaries give rise to strong frequency dependence, *reflectivity dispersion* and *tuning* effects which must be taken into account. To do this, in Section 3 we compute the *generalized* P-wave reflection coefficients by means of a reflectivity method. This leads us to the study of frequency dependent AVA and *amplitude vs. frequency AVF*, which are topics of great interest due to the increasing use of spectral decomposition techniques on seismic data. This also allows us to analyze the pattern of maxima in the reflectivity amplitude and their associated *peak frequencies*, which are strongly related to the thickness of the CO2 layer.

Seismic monitoring of geological CO2 injection sites is mostly based on the seismic reflections coming from high saturation CO2 accumulations. This is due to their large seismic amplitudes and good signal to noise ratio. However, low-saturation zones with dispersed CO2, or *saturation transition zones* may have an important role in the propagation of waves within the reservoir, giving rise to amplitude and phase changes of the seismic signals. Transition zones have been studied by different authors in the geophysical literature considering linear velocity trends with depth and constant density [10], [11]. Therefore, using the parameters of the Sleipner field, in Section 3.3 we model the reflectivity response of a CO2 transition layer, defined by a given linear vertical CO2 saturation profile, which results in a non-linear velocity trend with depth.

In general, emphasis is placed on establishing correlations between the seismic reflectivity and related attributes (intercept, gradient, curvature, peak frequencies) and the overall CO2 saturation, its physical state and thickness of the accumulation. These results are intended to help in the understanding of the expectable variations in a seismic time lapse study. They can be extended to other CO2 repositories with proper calibration of the rock and fluid properties.

Many of the results and procedures presented in this chapter are a revision and extension of those published by us in [4], [12], [13], [14] and [15].

## **1. Theoretical framework and assumptions**

2 Carbon Sequestration

controls the amplitude of seismic wave reflections and strongly conditions the detectability

With this idea, using rock physics models, in this chapter we focus on the theoretical analysis of the properties and variations of the seismic reflectivity and related parameters in a porous rock partially saturated with CO2 and brine. In our tests we model the seismic effects of

We calibrate our models using information on the Utsira formation, a shallow saline aquifer in the Sleipner field (offshore Norway). Millions of tonnes of CO2 have been separated since October 1996 from natural gas and re-injected into the aquifer, consisting of a highly porous and unconsolidated sandstone, with several thin intra-reservoir shale layers. These shale intervals act as temporary seals causing accumulations of CO2 beneath them [6],[7], whose saturations increase with time. A 50-100 m thick high velocity layer overlies the sand, forming the reservoir caprock, known as the Nørland formation. It is mainly composed of shales, siltstones and mudstones [8]. Another high-velocity shale layer of about 5 m thick is

As pointed out by [9] the *topmost* CO2 layer is of special interest for monitoring the Sleipner injection site, since its growth reveals the total upward flux and its changes over time. Thus, the analysis of the reflectivity associated with this kind of CO2 accumulation is an important topic in CO2 sequestration problems. The present book chapter focuses on this subject and

First, in Section 1 we present a description of the models and rock physics tools used in the subsequent sections. Next, in Section 2 we model the compressional P-wave reflection coefficient at the interface between a caprock and a permeable porous layer partially saturated by a mixture of brine and CO2 under liquid, supercritical and gaseous conditions. For this analysis we consider a simple model consisting of two halfspaces, combining a fluid substitution procedure and a Gassmann-Hill formulation to take into account variable spatial distributions of the fluids. We perform a sensitivity analysis of the standard AVA coefficients (intercept, gradient and curvature) in the near offset range, to investigate their correlations with CO2 saturation and its thermodynamic state within the geologic formation. Also, we analyze the effect of modeling the CO2 by means of simple and complex equations of state. The two-halfspace representation considered in Section 2 is valid when the accumulation thicknesses are larger than the involved seismic wavelengths. When this assumption is not valid, the interference between the multiple waves reflected at the different boundaries give rise to strong frequency dependence, *reflectivity dispersion* and *tuning* effects which must be taken into account. To do this, in Section 3 we compute the *generalized* P-wave reflection coefficients by means of a reflectivity method. This leads us to the study of frequency dependent AVA and *amplitude vs. frequency AVF*, which are topics of great interest due to the increasing use of spectral decomposition techniques on seismic data. This also allows us to analyze the pattern of maxima in the reflectivity amplitude and their associated *peak*

is structured in three main sections considering different aspects of the problem.

*frequencies*, which are strongly related to the thickness of the CO2 layer.

Seismic monitoring of geological CO2 injection sites is mostly based on the seismic reflections coming from high saturation CO2 accumulations. This is due to their large seismic amplitudes and good signal to noise ratio. However, low-saturation zones with dispersed CO2, or *saturation transition zones* may have an important role in the propagation of

of the injected CO2 volume, as reported by [3],[4] and [5] among others.

different CO2 accumulations formed below impermeable seals.

located about 10 m below the top of the sand.

In this section we summarize the modeling tools used later for the applications. First we explain the procedure for the calculation of the elastic properties of the partially saturated rocks for variable CO2 saturation. Next we describe the laws for the computation of the physical properties of the CO2 and brine for the different temperature and pressure states. Finally, we outline the reflectivity method for the computation of the generalized seismic reflectivity of the layered model

### **1.1. Elastic properties of CO**2 **bearing rocks**

To begin with the description of the model we consider that after the injection, a volume of carbon dioxide occupies part of the pore volume of a geologic reservoir which at the pre-injection state was fully saturated with brine. For simplicity we assume that CO2 displaces, without dissolution, part of the *in-situ* brine giving rise to a two-phase fluid saturation. From now on the corresponding saturations of free CO2 and brine are denoted as *S<sup>g</sup>* and *Sbr*, respectively, so that *Sbr* + *S<sup>g</sup>* = 1.

For the study of partially saturated rocks it is necessary to compute the bulk density and elastic coefficients (bulk and shear modulus) of the fluid saturated medium. The bulk density *ρ* is given by

$$
\rho = (1 - \phi)\rho\_s + \phi \left( S\_{br} \ \rho\_{br} + S\_g \ \rho\_g \right), \tag{1}
$$

where *φ* is the rock porosity, *ρ<sup>s</sup>* is the mineral grain density and *ρg*, *ρbr* are CO<sup>2</sup> and brine densities, respectively.

The formulation and solution of the energy and amplitude splitting problem when a monochromatic plane compressional wave strikes obliquely at a plane interface between two porous saturated media has been discussed by various authors, such as [16], [17] and [18]. In those papers the mechanical behavior of the porous rock was described using the classic constitutive relations and equations of motion of Biot [19],[20]. In the following sections our aim is to study the reflection of elastic waves and related parameters, for frequencies *f*

within the common seismic range *f* 120*Hz*. When a low frequency seismic wave propagates through a fluid saturated rock, due to fluid viscosity the solid and the fluid move in phase, so the medium behaves as an *effective* medium.

One of the most significant models to estimate the effective bulk modulus and the seismic velocity of a porous fluid saturated rock is to use Gassmann's relations [21] which can be written in the form:

$$K^G = K\_s \left(\frac{K\_m + Q}{K\_s + Q}\right), \quad \text{with } Q = \frac{K\_f(K\_s - K\_m)}{\phi(K\_s - K\_f)}.\tag{2}$$

In this equation the mechanical behavior of the fluid saturated porous medium is assumed to be elastic and isotropic. The coefficients *Km*, *K<sup>s</sup>* and *K<sup>f</sup>* are the bulk modulus of the dry matrix, the mineral grains and the pore fluid, respectively. The physical properties of the two-phase brine-CO2 fluid are computed using an *effective fluid* whose density is given by their weighted average and its compressibility by the isostress Reuss average of individual fluid bulk moduli [22]:

$$\frac{1}{K\_f} = \frac{S\_{br}}{K\_{br}} + \frac{S\_g}{K\_g},\tag{3}$$

where *Kbr* and *K<sup>g</sup>* are the bulk moduli of brine and CO2. In these equations, it is assumed that the mixture of CO2 and brine at the pore scale can be treated as a viscous single phase effective fluid. However, as pointed out in [23], this approach is strictly valid only when the pore fluids are uniformly mixed at very small scales so that the different wave-induced pressure increments in each fluid have time to diffuse and equilibrate during a seismic period. This critical distance is the so-called *diffusion length*. The presence of fluid heterogeneities over scales greater than this length (in which wave induced pore pressure gradients cannot equilibrate quickly), gives rise to a *patchy saturation*. It has been shown that for low frequencies such as those used in seismic exploration the effective modulus of a rock with patches of brine and CO2 of arbitrary geometry is given by [24],[23]

$$\frac{1}{\left(K^{P} + \frac{4}{3}\mu\right)} = \frac{S\_{br}}{K\_{br}^{G} + \frac{4}{3}\mu} + \frac{S\_{g}}{K\_{g}^{G} + \frac{4}{3}\mu};\tag{4}$$

where *K<sup>G</sup> br* and *<sup>K</sup><sup>G</sup> g* are the Gassmann's moduli (2) of the rock fully saturated with brine and CO2. The coefficient *µ* denotes the shear modulus of the rock, equal to that of the rock matrix, since the rigidity of a rock does not change due to the saturant fluid.

The existence of heterogeneous fluid distribution may give rise to mesoscopic wave attenuation and velocity dispersion phenomena in the seismic frequency range. There are many studies about numerical modeling of these effects, e.g. [25],[26],[27], showing that they are strongly dependent on the shapes and characteristic lengths of the patches. Given that these parameters are rarely known, to avoid dealing with this uncertainty, in our analysis we assume elastic layers where no attenuation-dispersion phenomena, associated with the patchy fluid distribution take place. However, these effects can be included in our formulation using appropriate viscoelastic models in the frequency domain. Because of the classic *correspondence principle*, this implies the replacement of the real constant elastic moduli by complex frequency dependent moduli [20].

The patchy (Hill) and uniform (Gassmann) fluid distributions are respectively upper and lower bounds for the seismic wave velocities [23]. Taking into account the uncertainties in the knowledge of the *in-situ* CO2 - brine distribution, in our models we assume that for each saturation it is reasonable to compute the average between both velocities.

It is important to mention the assumption of no chemical interactions between the pore fluids and the frame, which allows us to employ the fluid substitution procedure to consider that the pore space is saturated by brine and CO2 in variable proportions. Also, since the amount of CO2 dissolved in brine is a negligible fraction [28],[29], we do not take this effect into account. Nevertheless, when CO2 is injected in oil reservoirs this effect can not be neglected and should be included in the physical models, as shown in [4] and [28].

As pointed out in [8], no systematic variations in fluid pressure are observed in the Sleipner field. Thus, variations of rock elastic properties with effective pressure (related to the difference between confining and pore pressure), are not taken into account in the model. However, using appropriate effective pressure laws these effects can be included in the computations, as explained in [18].

### **1.2. Physical properties of carbon dioxide and brine**

4 Carbon Sequestration

written in the form:

fluid bulk moduli [22]:

where *K<sup>G</sup>*

*br* and *<sup>K</sup><sup>G</sup>*

so the medium behaves as an *effective* medium.

*K<sup>G</sup>* = *K<sup>s</sup>*

*K<sup>m</sup>* + *Q Ks* + *Q*

1 *K<sup>f</sup>*

patches of brine and CO2 of arbitrary geometry is given by [24],[23]

<sup>3</sup>*µ*) <sup>=</sup> *<sup>S</sup>br K<sup>G</sup> br* <sup>+</sup> <sup>4</sup> <sup>3</sup>*µ* +

matrix, since the rigidity of a rock does not change due to the saturant fluid.

and CO2. The coefficient *µ* denotes the shear modulus of the rock, equal to that of the rock

The existence of heterogeneous fluid distribution may give rise to mesoscopic wave attenuation and velocity dispersion phenomena in the seismic frequency range. There are many studies about numerical modeling of these effects, e.g. [25],[26],[27], showing that they are strongly dependent on the shapes and characteristic lengths of the patches. Given that these parameters are rarely known, to avoid dealing with this uncertainty, in our analysis we assume elastic layers where no attenuation-dispersion phenomena, associated with the patchy fluid distribution take place. However, these effects can be included in our formulation using appropriate viscoelastic models in the frequency domain. Because of the

1 (*K<sup>P</sup>* + <sup>4</sup>

within the common seismic range *f* 120*Hz*. When a low frequency seismic wave propagates through a fluid saturated rock, due to fluid viscosity the solid and the fluid move in phase,

One of the most significant models to estimate the effective bulk modulus and the seismic velocity of a porous fluid saturated rock is to use Gassmann's relations [21] which can be

In this equation the mechanical behavior of the fluid saturated porous medium is assumed to be elastic and isotropic. The coefficients *Km*, *K<sup>s</sup>* and *K<sup>f</sup>* are the bulk modulus of the dry matrix, the mineral grains and the pore fluid, respectively. The physical properties of the two-phase brine-CO2 fluid are computed using an *effective fluid* whose density is given by their weighted average and its compressibility by the isostress Reuss average of individual

> <sup>=</sup> *<sup>S</sup>br Kbr* + *Sg Kg*

where *Kbr* and *K<sup>g</sup>* are the bulk moduli of brine and CO2. In these equations, it is assumed that the mixture of CO2 and brine at the pore scale can be treated as a viscous single phase effective fluid. However, as pointed out in [23], this approach is strictly valid only when the pore fluids are uniformly mixed at very small scales so that the different wave-induced pressure increments in each fluid have time to diffuse and equilibrate during a seismic period. This critical distance is the so-called *diffusion length*. The presence of fluid heterogeneities over scales greater than this length (in which wave induced pore pressure gradients cannot equilibrate quickly), gives rise to a *patchy saturation*. It has been shown that for low frequencies such as those used in seismic exploration the effective modulus of a rock with

, with *<sup>Q</sup>* <sup>=</sup> *<sup>K</sup><sup>f</sup>* (*K<sup>s</sup>* <sup>−</sup> *<sup>K</sup>m*)

*Sg K<sup>G</sup> <sup>g</sup>* <sup>+</sup> <sup>4</sup> <sup>3</sup>*µ*

*g* are the Gassmann's moduli (2) of the rock fully saturated with brine

*<sup>φ</sup>*(*K<sup>s</sup>* <sup>−</sup> *<sup>K</sup><sup>f</sup>* ) . (2)

, (3)

; (4)

Depending on the *in-situ* pressure and temperature conditions the CO2 can exist in the subsurface in different phases. We recall that the *critical point* for CO2 occurs at a temperature *<sup>T</sup><sup>c</sup>* = 31.1 ◦C and a pressure *<sup>P</sup><sup>c</sup>* = 7.39 MPa. For temperatures higher than *T<sup>c</sup>* and pressures higher than *P<sup>c</sup>* the carbon dioxide is in a *supercritical* state, where it is compressible like a gas but with the density of a liquid. This characteristic of CO2 is particularly relevant for its underground storage since supercritical CO2 can fill the available pore volume with minimum buoyancy effects [30]. Temperatures and pressures near the critical point commonly occur in applications involving CO2, such as enhanced oil recovery techniques and sequestration projects [31]. However, as pointed out by [30], the depth at which CO2 supercritical conditions are present is highly variable and strongly dependent on surface temperature and geothermal gradients, even within a single basin. In addition, the pressure regime of the basin (normal or abnormal), is also very important and is related to its geologic history, existence of sealing faults, permeability barriers and the occurrence of overpressure generation mechanisms [32].

For the examples, the density and bulk moduli of brine for given *in-situ* temperature and pressure conditions are computed using the semi-empirical relations proposed by Batzle and Wang (BW) [29]. For the computations we consider a typical brine salinity of 50000 ppm. The corresponding properties of CO2 can be computed using some of the many equations of state (EoS) developed for real gases such as the two-parameter van der Waals (vW) equation [35] and the Peng and Robinson (PR) equation [34], or some more specific one such as Duan, Moller and Weare (DMW) [33], which involves fifteen parameters. To our knowledge, there is no full agreement in the geophysical literature about which EoS is the most appropriate to represent CO2 properties under the thermodynamic conditions found in geologic reservoirs. Therefore, we use different models to analyze the effect of the EoS in the seismic magnitudes of interest. With that purpose we selected three different temperature


**Table 1.** Physical properties of CO2 and brine under different thermodynamic conditions. The density and bulk moduli of CO2 are obtained using three equations of state: DMW [33], PR [34] and vW [35]. The CO2 bulk moduli estimations are corrected for adiabatic conditions. The properties of brine are computed using BW laws [29] for a salinity of 50000 ppm.

and pressure (*T*, *<sup>P</sup>*) values: (40◦C, 6 MPa) corresponding to gaseous state, (36◦C, 10 MPa) corresponding to supercritical state and (20◦C, 10 MPa) corresponding to liquid state. The resulting estimations are shown in Table 1 and will be used in the following sections. The temperature and pressure of the supercritical state are in agreement with the conditions reported in [8] for Utsira sandstone.

It must be pointed out that the different CO2 bulk moduli estimations in Table 1, to be used for the computations, were corrected to represent adiabatic conditions. The reason is that the fluid compression due to the passage of a wave is fast, so that the process is not isothermal [29], [37].

#### **1.3. Calibration of the rock physics model**

To calibrate our elastic model for the Utsira sandstone, we used the measured bulk density and the compressional and shear wave velocities in the preinjection state (i.e. under full brine saturation), given in [8] and [38]: *<sup>ρ</sup><sup>b</sup>* <sup>=</sup> 2050 kg/m3, *<sup>V</sup><sup>p</sup>* =2050 m/s and *<sup>V</sup><sup>s</sup>* <sup>=</sup> <sup>640</sup> m/s. Using the parameters given for the pore water *<sup>ρ</sup>br* <sup>=</sup> 1040 kg/m3, *<sup>K</sup>br* <sup>=</sup> 2.305 GPa and the average sandstone porosity *φ* = 0.37, the mineral grain density results in *ρs* = 2643 kg/m3. We also performed a Gassmann-inverse calculation to obtain: *<sup>K</sup><sup>m</sup>* <sup>=</sup> 2.569 GPa and *µ* = 0.84 GPa. For the shale layers we obtained the elastic parameters using the reported values: *<sup>V</sup><sup>p</sup>* <sup>=</sup> 2270 m/s, *<sup>V</sup><sup>s</sup>* <sup>=</sup> 850 m/s and *<sup>ρ</sup><sup>b</sup>* <sup>=</sup> 2100 kg/m3.

#### **1.4. The generalized reflectivity of a layered medium**

The computation of the reflection and transmission coefficients of elastic waves propagating through layered media is a subject of interest in different fields such as seismology, seismic prospecting and underwater acoustics. The first numerical procedures date back to the pioneering works of Thomson [39] and Haskell [40], who proposed matrix methods that transfer stresses and displacements through successive layers. Since then, different approaches and applications have been presented.

In this work we perform an implementation of the *reflectivity method* [41], [42], which is based on the continuity of particle displacements and stress components at the interfaces between sets of plane layers embedded between two halfspaces. This method builds up the reflection and transmission matrices iteratively by starting at the top of a lower bounding halfspace and adding one layer per iteration until the total stack response is constructed. This procedure is intuitively simple and exact [43], taking into account all internal reverberations. The recursion algorithm is

6 Carbon Sequestration

*Gaseous*

*Liquid*

*Supercritical*

reported in [8] for Utsira sandstone.

**1.3. Calibration of the rock physics model**

values: *<sup>V</sup><sup>p</sup>* <sup>=</sup> 2270 m/s, *<sup>V</sup><sup>s</sup>* <sup>=</sup> 850 m/s and *<sup>ρ</sup><sup>b</sup>* <sup>=</sup> 2100 kg/m3.

approaches and applications have been presented.

**1.4. The generalized reflectivity of a layered medium**

isothermal [29], [37].

**Physical state CO**2 **CO**2 **CO**2 **Brine EoS DMW PR vW BW**

*<sup>T</sup>* <sup>=</sup> <sup>40</sup>◦ <sup>C</sup> *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.0049 *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.0089 *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.0050 *<sup>K</sup>br* <sup>=</sup> 2.5986 GPa *<sup>P</sup>* <sup>=</sup> 6 MPa *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 149.8 *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 153.2 *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 142.1 *<sup>ρ</sup>br* <sup>=</sup> 1028.7 kg/m<sup>3</sup>

*<sup>T</sup>* <sup>=</sup> <sup>36</sup>◦ <sup>C</sup> *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.037245 *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.026120 *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.01999 *<sup>K</sup>br* <sup>=</sup> 2.6234 GPa *<sup>P</sup>* <sup>=</sup> 10 MPa *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 637.6 *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 706.4 *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 490.7 *<sup>ρ</sup>br* <sup>=</sup> 1030.4 kg/m<sup>3</sup>

*<sup>T</sup>* <sup>=</sup> <sup>20</sup>◦ <sup>C</sup> *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.0931 *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.1376 *<sup>K</sup><sup>g</sup>* <sup>=</sup> 0.0464 *<sup>K</sup>br* <sup>=</sup> 2.5009 GPa *<sup>P</sup>* <sup>=</sup> 10 MPa *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 853.5 *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 833.1 *<sup>ρ</sup><sup>g</sup>* <sup>=</sup> 565.2 *<sup>ρ</sup>br* <sup>=</sup> 1036.07 kg/m<sup>3</sup>

**Table 1.** Physical properties of CO2 and brine under different thermodynamic conditions. The density and bulk moduli of CO2 are obtained using three equations of state: DMW [33], PR [34] and vW [35]. The CO2 bulk moduli estimations are corrected

and pressure (*T*, *<sup>P</sup>*) values: (40◦C, 6 MPa) corresponding to gaseous state, (36◦C, 10 MPa) corresponding to supercritical state and (20◦C, 10 MPa) corresponding to liquid state. The resulting estimations are shown in Table 1 and will be used in the following sections. The temperature and pressure of the supercritical state are in agreement with the conditions

It must be pointed out that the different CO2 bulk moduli estimations in Table 1, to be used for the computations, were corrected to represent adiabatic conditions. The reason is that the fluid compression due to the passage of a wave is fast, so that the process is not

To calibrate our elastic model for the Utsira sandstone, we used the measured bulk density and the compressional and shear wave velocities in the preinjection state (i.e. under full brine saturation), given in [8] and [38]: *<sup>ρ</sup><sup>b</sup>* <sup>=</sup> 2050 kg/m3, *<sup>V</sup><sup>p</sup>* =2050 m/s and *<sup>V</sup><sup>s</sup>* <sup>=</sup> <sup>640</sup> m/s. Using the parameters given for the pore water *<sup>ρ</sup>br* <sup>=</sup> 1040 kg/m3, *<sup>K</sup>br* <sup>=</sup> 2.305 GPa and the average sandstone porosity *φ* = 0.37, the mineral grain density results in *ρs* = 2643 kg/m3. We also performed a Gassmann-inverse calculation to obtain: *<sup>K</sup><sup>m</sup>* <sup>=</sup> 2.569 GPa and *µ* = 0.84 GPa. For the shale layers we obtained the elastic parameters using the reported

The computation of the reflection and transmission coefficients of elastic waves propagating through layered media is a subject of interest in different fields such as seismology, seismic prospecting and underwater acoustics. The first numerical procedures date back to the pioneering works of Thomson [39] and Haskell [40], who proposed matrix methods that transfer stresses and displacements through successive layers. Since then, different

In this work we perform an implementation of the *reflectivity method* [41], [42], which is based on the continuity of particle displacements and stress components at the interfaces between sets of plane layers embedded between two halfspaces. This method builds up the reflection

for adiabatic conditions. The properties of brine are computed using BW laws [29] for a salinity of 50000 ppm.

$$\begin{aligned} \mathcal{R}\_{j-1} &= R\_{j-1}^d + T\_{j-1}^u \bar{\mathcal{R}}\_j (I - R\_{j-1}^u \bar{\mathcal{R}}\_j)^{-1} T\_{j-1}^d, \\ \mathcal{T}\_{j-1} &= \bar{\mathcal{T}}\_j (I - R\_{j-1}^u \bar{\mathcal{R}}\_j)^{-1} T\_{j-1}^d; \\ \bar{\mathcal{R}}\_j &= E\_j \mathcal{R}\_j E\_j, \\ \bar{\mathcal{T}}\_j &= \mathcal{T}\_j E\_j, \end{aligned} \tag{5}$$

were R*<sup>j</sup>* and T*<sup>j</sup>* are the total-reflection and total-transmission matrices of interface *j*. *<sup>R</sup>u*,*<sup>d</sup> <sup>j</sup>* and *<sup>T</sup> <sup>u</sup>*,*<sup>d</sup> <sup>j</sup>* are the upward and downward layer reflection and transmission matrices respectively. These matrices depend on the bulk density, P- and SV-wave seismic velocities and ray angles from interface *j*. *E<sup>j</sup>* is the layer phase-shift diagonal matrix and (*I* − *R<sup>u</sup> <sup>j</sup>*−1R¯ *<sup>j</sup>* ) <sup>−</sup><sup>1</sup> is known as the reverberation operator [44], with *<sup>I</sup>* being the identity matrix.

The recursion starts at the base of the layering *j* = *n* with R*<sup>n</sup>* = *R<sup>d</sup> <sup>n</sup>* = 0 and T*<sup>n</sup>* = *T<sup>d</sup> <sup>n</sup>* = *I* and continues to *j* = 1. When the top of the layering is reached, the generalized P-wave reflection coefficient for scaled displacements is obtained from *<sup>R</sup>pp*(*f*, *<sup>θ</sup>*) = <sup>R</sup>0[1,1]. This coefficient is a complex function of frequency *f* and incidence angle *θ*.

For a complete description of the matrices involved in this scheme the reader is referred to [44].

## **2. AVA analysis at the top of a CO**2 **accumulation**

In this section we model the seismic compressional wave reflection coefficient and related AVA parameters at the interface between Utsira sandstone, partially saturated by mixtures of CO2 and brine, overlain by a shale caprock, with the properties described in Section 1.3.

Following the classic approach, for the present applications both media are considered as elastic halfspaces, an assumption valid for layers of thicknesses larger than the wavelengths of the incident waves. Under these conditions the reflection coefficient are real and frequency independent. In this particular configuration, our reflectivity method reproduces the classic results of Zoeppritz [45] for two halfspaces.

In the near offset domain *θ <* <sup>60</sup>◦ or below the critical angle, we can assume that the *<sup>R</sup>pp* reflection coefficient as a function of the incidence angle *θ* can be approximated using the expression of Shuey [46]:

$$R\_{pp}(\theta) \approx A + B\sin^2\theta + C\left(\tan^2\theta - \sin^2\theta\right);\tag{6}$$

where *A*, *B* and *C* are known as the *AVA coefficients*. *A* is called the *intercept*, *B* the *gradient* and *C* the *curvature* [47],[48]. The intercept is equal to the normal incidence reflection coefficient and is controlled by the contrast in acoustic impedance between both media. The gradient is related to the contrast in density and in compressional P and shear SV wave

**Figure 1.** Compressional elastic wave velocity for variable CO<sup>2</sup> saturation *S<sup>g</sup>* obtained using Gassmann-Hill average under different thermodynamic conditions.

**Figure 2.** P-wave reflection coefficient vs. ray angle for different fixed partial CO2 saturations under supercritical conditions.

velocities [49]. The curvature parameter is important at far offsets and near critical angles [44]. Equation (6) is used in this section to carry out a parametric analysis of the AVA coefficients or to study their sensitivity by implementing a standard fitting procedure on the results obtained for *Rpp*(*θ*).

In Figure 1 we plot the compressional elastic wave velocity for the partially saturated CO2-brine Utsira sandstone in the low frequency range. The parameters of the CO2 at the liquid, supercritical and gaseous states correspond to the PR EoS and are listed in Table 1, along with those of the brine. As mentioned in Section 1, our velocity model for the sandstone corresponds to the Gassmann-Hill average for each saturation state. It is important to note the decreasing tendency of the wave velocity for increasing CO2 saturation, which becomes less marked as saturation increases, a feature that will also be observed in the AVA coefficients. The velocity model for the supercritical CO2 is in close agreement with that presented in [8] and [38] for the Utsira sandstone. According to this model, the seismic velocity can discriminate the liquid from the other two physical states, particularly for saturations lower than about 40%.

Next, in Figure 2 a)-d) we show the compressional reflection coefficient versus incidence angle from 0◦ to 50◦, for fixed values of CO<sup>2</sup> saturation in the range 0%–100%. The temperature is 36◦C and pressure 10 MPa, corresponding to the supercritical conditions given in Table 1, for the PR equation. As expected, an abrupt change in the reflectivity from the pre-injection to the post-injection state is observed, due to the high contrast between the physical properties of CO2 and brine in this state. As expected, the CO2 saturation decreases the impedance of the lower medium with respect to that of the upper one, resulting in a negative reflection coefficient for normal incidence, becoming more negative for higher angles and higher saturations. These curves are consistent with those published in [50]. The curves corresponding to liquid and gaseous states show very similar behavior, and consequently they are not shown. This means that the discrimination of the physical state from the AVA behavior may be difficult. This subject will be analyzed in more detail in terms of the AVA coefficients. According to the well known AVA classification [48], [51], the curves for the

8 Carbon Sequestration

different thermodynamic conditions.

results obtained for *Rpp*(*θ*).

for saturations lower than about 40%.

**Figure 1.** Compressional elastic wave velocity for variable CO<sup>2</sup> saturation *S<sup>g</sup>* obtained using Gassmann-Hill average under

**Figure 2.** P-wave reflection coefficient vs. ray angle for different fixed partial CO2 saturations under supercritical conditions. velocities [49]. The curvature parameter is important at far offsets and near critical angles [44]. Equation (6) is used in this section to carry out a parametric analysis of the AVA coefficients or to study their sensitivity by implementing a standard fitting procedure on the

In Figure 1 we plot the compressional elastic wave velocity for the partially saturated CO2-brine Utsira sandstone in the low frequency range. The parameters of the CO2 at the liquid, supercritical and gaseous states correspond to the PR EoS and are listed in Table 1, along with those of the brine. As mentioned in Section 1, our velocity model for the sandstone corresponds to the Gassmann-Hill average for each saturation state. It is important to note the decreasing tendency of the wave velocity for increasing CO2 saturation, which becomes less marked as saturation increases, a feature that will also be observed in the AVA coefficients. The velocity model for the supercritical CO2 is in close agreement with that presented in [8] and [38] for the Utsira sandstone. According to this model, the seismic velocity can discriminate the liquid from the other two physical states, particularly

**Figure 3.** Behavior of AVA parameters intercept *A*, gradient *B* and curvature *C* versus CO<sup>2</sup> saturation in different physical states.

preinjection case and for *S<sup>g</sup>* = 10% (in Figure 2 b) correspond to the AVA class IV. The set of AVA curves for higher saturations shows a very flat tendency for near angles (i.e. with very low gradients *B*), which can be observed in more detail in Figures 2 c) and d). Its classification will be discussed later in Figure 4.

To analyze the sensitivity of the AVA parameters we plot in Figure 3 the predicted variations of *A*, *B* and *C* with increasing saturation for the liquid, supercritical and gaseous conditions. We observe that the intercept *A* is negative, decreasing and monotone for the three physical states under consideration, being a sensitive attribute throughout the saturation range. The gradient *B* shows a trend similar to that of the intercept except for the gaseous state, in which stabilization is observed. This can be useful to distinguish this particular state from the other two.

The curvature attribute *C*, is the least sensitive of the three, both to physical state and saturation level. We only note appreciable changes for liquid CO2 and for saturations in the range from 0 to about 30%. However, the correct determination of such changes may be strongly limited by the seismic resolution. In connection with this, Brown et al. [3] pointed out that AVA variations on the order of 5% are seismically detectable. Taking that numerical threshold into account we can state that, for this model, variations in the parameter *C* do not contain much information with respect to variations in the CO<sup>2</sup> content. Consequently, an AVO time-lapse analysis can be based solely on the coefficients *A* and *B*,

where *A* is the most sensitive for the determination of the CO<sup>2</sup> saturation state. To inspect

**Figure 4.** Intercept (A) vs. gradient (B) crossplot for variable saturation and physical states.

the behavior of the reflectivity in terms of the AVA classification, in Figure 4 we present typical intercept-gradient crossplots (*B* vs. *A*) to observe their trends for variable saturations and physical states. According to our model, while for gaseous CO2 the AVA class is IV in all the saturation range, under liquid and supercritical conditions a change to class III can be expected for very high saturations. We also observe that since from intermediate to high gas saturations, the *B* − *A* trajectories of each CO<sup>2</sup> state are different, a time-lapse crossplot analysis may bring useful information about the physical state of the underground CO2. In order to analyze the influence of the state equations in our model estimates, in Figure 5 we display the coefficients *A* and *B* as a function of CO<sup>2</sup> saturation, for van der Waals [35], Peng & Robinson [34] and Duan et al. [33] EoS for the three physical states in Table 1. For gaseous and supercritical states the coefficients are very similar in all the saturation range. The most important discrepancies are observed when CO2 is at liquid state. From these experiments we conclude that for seismic modeling applications the choice of the EoS will not significantly affect the results.

The sensitivity and robustness of the AVA coefficients to the presence of low fractions of methane (CH4), which is a frequent impurity in CO2 injection sites, was considered in [15].

#### **3. Finite thickness CO**2 **accumulations: frequency dependent AVA**

So far we have modeled the CO2 reflectivity response in terms of the simple single layer representation. However, due to the spatial variations in petrophysical properties and natural permeability barriers, the CO2 tends to form accumulations of finite thickness.

In exploration geophysics, it is generally accepted that when the thickness of a layer is smaller than the predominant seismic wavelengths, the layer is said to be *thin*. Also, a common measure for the vertical seismic resolution is given by a quarter of the dominant wavelength [53]. As explained in [6],[8] and [54], the CO2 accumulations within the Utsira sandstone in the Sleipner field can be regarded as thin layers. The interference effects between the

<sup>352</sup> CO2 Sequestration and Valorization Seismic Reflectivity in Carbon Dioxide Accumulations: a Review 11 10.5772/57087 Seismic Reflectivity in Carbon Dioxide Accumulations: A Review http://dx.doi.org/10.5772/57087 353

**Figure 5.** Intercept and gradient model estimates for variable saturation and physical states using different equations of state: DMW in blue, PR in red and vW in green lines.

multiple waves reflected at the different boundaries give rise to strong frequency dependence in the reflectivity of the medium, a phenomenon known as *reflectivity dispersion* [11],[55]. In this case the dispersive character has a geometric origin and is due to layer reverberation; however reflectivity dispersion can also take place at interfaces between highly dissipative media. The constructive interference between some frequency components of the up-going and down-going reflections at the physical boundaries gives rise to *tuning* effects. In the context of CO2 sequestration these effects have been analyzed in several papers, such as [5], [8],[9], and [38].

As pointed out in [11] and [55], the modern use of time-frequency analysis and spectral decomposition tools [56] has shown that reflection events in practice are always frequency dependent. The spectral decomposition techniques allow obtaining and analyzing the seismic response of a medium at a given frequency. This makes it possible to investigate the reflected amplitude variation with angle at different frequencies (*spectral AVA*) and the amplitude variation with frequency (*AVF*) at a fixed angle.

At the time of this writing and to our knowledge, the AVF characteristics of CO2 layers have not been modeled yet. With this motivation, in this section we investigate the existing correlations between reflectivity, frequency, saturation, physical state and thickness for a simple accumulation model. From these experiments our aim is to assess whether time lapse changes in the reflected amplitudes can be correlated with changes in such parameters, which are of utmost importance in CO2 monitoring.

#### **3.1. Description of the model**

10 Carbon Sequestration

where *A* is the most sensitive for the determination of the CO<sup>2</sup> saturation state. To inspect

the behavior of the reflectivity in terms of the AVA classification, in Figure 4 we present typical intercept-gradient crossplots (*B* vs. *A*) to observe their trends for variable saturations and physical states. According to our model, while for gaseous CO2 the AVA class is IV in all the saturation range, under liquid and supercritical conditions a change to class III can be expected for very high saturations. We also observe that since from intermediate to high gas saturations, the *B* − *A* trajectories of each CO<sup>2</sup> state are different, a time-lapse crossplot analysis may bring useful information about the physical state of the underground CO2. In order to analyze the influence of the state equations in our model estimates, in Figure 5 we display the coefficients *A* and *B* as a function of CO<sup>2</sup> saturation, for van der Waals [35], Peng & Robinson [34] and Duan et al. [33] EoS for the three physical states in Table 1. For gaseous and supercritical states the coefficients are very similar in all the saturation range. The most important discrepancies are observed when CO2 is at liquid state. From these experiments we conclude that for seismic modeling applications the choice of the EoS will

The sensitivity and robustness of the AVA coefficients to the presence of low fractions of methane (CH4), which is a frequent impurity in CO2 injection sites, was considered in [15].

**3. Finite thickness CO**2 **accumulations: frequency dependent AVA** So far we have modeled the CO2 reflectivity response in terms of the simple single layer representation. However, due to the spatial variations in petrophysical properties and natural

In exploration geophysics, it is generally accepted that when the thickness of a layer is smaller than the predominant seismic wavelengths, the layer is said to be *thin*. Also, a common measure for the vertical seismic resolution is given by a quarter of the dominant wavelength [53]. As explained in [6],[8] and [54], the CO2 accumulations within the Utsira sandstone in the Sleipner field can be regarded as thin layers. The interference effects between the

permeability barriers, the CO2 tends to form accumulations of finite thickness.

**Figure 4.** Intercept (A) vs. gradient (B) crossplot for variable saturation and physical states.

not significantly affect the results.

Here we present another simple but useful approach to model the reflectivity of the topmost CO2 layer in the Sleipner field. For the computations we use the sandstone and caprock

**Figure 6.** Generalized complex reflection coefficient *Rpp* vs. angle for f= 50 Hz, *S<sup>g</sup>* = 100% and thicknesses *h* = 0, 5, 10 m and *h* → ∞.

parameters described in Section 1.3. To simulate the supercritical conditions of the CO2 injected into Utsira sandstone, carbon dioxide density and bulk moduli are calculated using the PR EoS at a temperature of 36◦C and pressure of 10 MPa (Table 1).

According to the scheme used for the reflectivity method (described in Section 1.4), the upper halfspace plays the role of the caprock and the lower halfspace is the fully brine saturated Utsira sandstone. Below the caprock, the injected carbon dioxide is assumed to form a layer in the sandstone of thickness *h* with spatially constant CO2 saturation *Sg*. The case of vertical variations in the CO<sup>2</sup> saturation or *saturation transitions* will be considered in Section 3.3. Thus for this model the generalized reflection coefficient is a complex, frequency-dependent function of the incidence angle, the CO2 saturation and thickness *h* in the CO2 bearing layer.

#### **3.2. Behavior of reflectivity with saturation, frequency and thickness**

In Figure 6 we display the generalized reflection coefficient *Rpp* vs. incidence angle for a frequency of *f* = 50 Hz, a saturation of *Sg* = 100 % and thicknesses of *h* = 0, 5, 10 m, considering also the limit *h* → ∞. For this CO<sup>2</sup> saturation the P-wave velocity (see Figure 1) is *V<sup>p</sup>* ≈ 1410 m/s, so that the associated wavelength *λ* = *Vp*/*f* is about 28 m and the vertical resolution limit, *λ*/4, at this frequency results in about 7 m. So, in these examples the 5 m layer is ultra-thin and the 10 m layer is thin.

When *h* = 0 and *h* = ∞, the model reduces to the classic two halfspace model, with a real reflection coefficient (shown in Figures 6 a) and b)), which can be analyzed in terms of the AVA classification. As can be seen in Figures 6 b) and c), the spectral AVA response for finite thickness layers differs considerably from that of the halfspaces, as expected. For 0 *<h<* ∞ the AVA character is highly dependent on frequency, a dependence that can be clearly observed in Figure 7. In particular, the curves for the thicknesses *h* = 5 m and *h* = 10 m show a periodic behavior; meaning that the different frequency components of a seismic signal are reflected with different amplitudes. To understand this periodicity on a theoretical basis, it is useful to consider the analytical result for the reflectivity of an acoustic layer embedded between two halfspaces, treated in [52] and [53]. From that solution it can

10.5772/57087

12 Carbon Sequestration

and *h* → ∞.

thickness *h* in the CO2 bearing layer.

the 5 m layer is ultra-thin and the 10 m layer is thin.

**Figure 6.** Generalized complex reflection coefficient *Rpp* vs. angle for f= 50 Hz, *S<sup>g</sup>* = 100% and thicknesses *h* = 0, 5, 10 m

parameters described in Section 1.3. To simulate the supercritical conditions of the CO2 injected into Utsira sandstone, carbon dioxide density and bulk moduli are calculated using

According to the scheme used for the reflectivity method (described in Section 1.4), the upper halfspace plays the role of the caprock and the lower halfspace is the fully brine saturated Utsira sandstone. Below the caprock, the injected carbon dioxide is assumed to form a layer in the sandstone of thickness *h* with spatially constant CO2 saturation *Sg*. The case of vertical variations in the CO<sup>2</sup> saturation or *saturation transitions* will be considered in Section 3.3. Thus for this model the generalized reflection coefficient is a complex, frequency-dependent function of the incidence angle, the CO2 saturation and

**3.2. Behavior of reflectivity with saturation, frequency and thickness**

In Figure 6 we display the generalized reflection coefficient *Rpp* vs. incidence angle for a frequency of *f* = 50 Hz, a saturation of *Sg* = 100 % and thicknesses of *h* = 0, 5, 10 m, considering also the limit *h* → ∞. For this CO<sup>2</sup> saturation the P-wave velocity (see Figure 1) is *V<sup>p</sup>* ≈ 1410 m/s, so that the associated wavelength *λ* = *Vp*/*f* is about 28 m and the vertical resolution limit, *λ*/4, at this frequency results in about 7 m. So, in these examples

When *h* = 0 and *h* = ∞, the model reduces to the classic two halfspace model, with a real reflection coefficient (shown in Figures 6 a) and b)), which can be analyzed in terms of the AVA classification. As can be seen in Figures 6 b) and c), the spectral AVA response for finite thickness layers differs considerably from that of the halfspaces, as expected. For 0 *<h<* ∞ the AVA character is highly dependent on frequency, a dependence that can be clearly observed in Figure 7. In particular, the curves for the thicknesses *h* = 5 m and *h* = 10 m show a periodic behavior; meaning that the different frequency components of a seismic signal are reflected with different amplitudes. To understand this periodicity on a theoretical basis, it is useful to consider the analytical result for the reflectivity of an acoustic layer embedded between two halfspaces, treated in [52] and [53]. From that solution it can

the PR EoS at a temperature of 36◦C and pressure of 10 MPa (Table 1).

**Figure 7.** Modulus of the generalized reflection coefficient *Rpp* vs. frequency for normal incidence, *S<sup>g</sup>* = 100 % and thicknesses *h* = 0, 5, 10 m and *h* → ∞.

be shown that the function *Rpp*(*f*, *θ*) has multiple maxima or *peaks* at frequencies *f<sup>n</sup> max* and minima or *notches* located at frequencies *f<sup>n</sup> min*, given respectively by [15]

$$f\_{max}^{n} = \left(\frac{1+2n}{2}\right)\frac{V\_p}{2h\cos\theta\_L}, \quad f\_{min}^{n} = n\frac{V\_p}{2h\cos\theta\_L}, \quad \text{for} \quad n = 0, 1, 2, \cdots \quad \text{and} \quad \binom{n}{2}$$

The occurrences of these extrema are directly proportional to the compressional velocity within the layer *Vp* (shown in Figure 1) and inversely proportional to its thickness *h*, where *θ<sup>L</sup>* is the ray angle inside the layer. The frequencies at which the reflected amplitude has maxima are known as *tuning frequencies* [58] or *peak frequencies*. The first peak frequency *f*0 *max* for normal incidence, hereafter denoted as *fp*, results in

$$f\_p = \frac{V\_p}{4h}.\tag{8}$$

Equation (8) is a very useful relation for thickness determination by means of the *fp* attribute, as shown in [38] in the context of CO2 sequestration.

For an arbitrary number of layers with different physical properties and variable thicknesses the generalized reflection coefficient also displays a pattern of maxima and minima. However, for such cases there are no simple analytical expressions and consequently a numerical procedure such as the one presented in this work enables predicting the spectral characteristics of their reflectivity.

For the particular cases *h* = 0 and *h* → ∞, plotted in Figures 7 a) and d), the reflectivity for normal incidence is real, constant and related to the impedance contrast between the caprock and the fully brine saturated sandstone and the fully CO2 saturated layer. In Figures 7 b) and c), the location of the peak close to 70 Hz for *h* = 5 m and the peaks close to 35 Hz and 105 Hz for *h* = 10 m can be verified. Note also that for the continuous component *f* = 0 Hz, these AVF curves converge to the response in Figure 7 a), given that for such a long wavelength the layers *h* = 5 m and *h* = 10 m are invisible to the seismic waves.

**Figure 8.** a) Modulus of *Rpp* vs. frequency and thickness for normal incidence and saturation *S<sup>g</sup>* = 100%. b) Relation between the first peak frequency *f<sup>p</sup>* and layer thickness *h*.

**Figure 9.** Modulus of *Rpp* vs. frequency and CO<sup>2</sup> saturation *S<sup>g</sup>* = for thickness *h* = 10 m and three incidence angles *θ* = 0◦, 30◦ and 60◦.

The modulus of *Rpp* versus frequency and thickness for normal incidence is shown in Figure 8 a), in which the reciprocal relation between peak frequencies and thicknesses is clearly displayed. As seen, at each thickness we have a different peak frequency. To study this dependency in more detail, in Figure 8 b) we display the first peak frequency versus thickness, in which it is instructive to note that for a frequency of 35 Hz the corresponding thickness *h* is about 10 m. This is the case studied in [38], where spectral decomposition is used at Sleipner to map the topmost CO2 accumulation. For very small thicknesses Figure 8 b) also indicates that the peak frequencies are well above the standard seismic frequency range, so high-resolution seismic acquisition is essential for these layers.

To analyze the sensitivity of the previous results to CO2 saturation, next in Figure 9, we plot the modulus of *Rpp* versus frequency and saturation for thickness *h* = 10 m and incidence angles 0◦, 30◦ and 60◦. The peaks of the reflectivity are located in the high CO<sup>2</sup> saturation range, since in that case the impedance contrast between the caprock and the sandstone is maximum. The sensitivity to *S<sup>g</sup>* depends on frequency showing bands of very weak reflectivity (in blue scale) and higher amplitudes (in red scale). The high reflectivity zones, which at the same time are sensitive to saturation, are located around the first and second peak frequencies for this example (*h* = 10 m), i.e. near 35 Hz and 105 Hz. These conclusions

**Figure 10.** a) The carbon dioxide saturation profile and its effects on the Utsira sandstone P-wave velocity b) and bulk density c) for *h*<sup>0</sup> = 10 m and *h<sup>t</sup>* = 5 m.

are valid for normal incidence *<sup>θ</sup>* = <sup>0</sup>◦ and also for larger angles, although a slight shift of the peak frequencies to higher values is observed, as expected from equation (7).

#### **3.3. The influence of saturation transitions**

14 Carbon Sequestration

0◦, 30◦ and 60◦.

between the first peak frequency *f<sup>p</sup>* and layer thickness *h*.

**Figure 8.** a) Modulus of *Rpp* vs. frequency and thickness for normal incidence and saturation *S<sup>g</sup>* = 100%. b) Relation

**Figure 9.** Modulus of *Rpp* vs. frequency and CO<sup>2</sup> saturation *S<sup>g</sup>* = for thickness *h* = 10 m and three incidence angles *θ* =

The modulus of *Rpp* versus frequency and thickness for normal incidence is shown in Figure 8 a), in which the reciprocal relation between peak frequencies and thicknesses is clearly displayed. As seen, at each thickness we have a different peak frequency. To study this dependency in more detail, in Figure 8 b) we display the first peak frequency versus thickness, in which it is instructive to note that for a frequency of 35 Hz the corresponding thickness *h* is about 10 m. This is the case studied in [38], where spectral decomposition is used at Sleipner to map the topmost CO2 accumulation. For very small thicknesses Figure 8 b) also indicates that the peak frequencies are well above the standard seismic frequency range, so

To analyze the sensitivity of the previous results to CO2 saturation, next in Figure 9, we plot the modulus of *Rpp* versus frequency and saturation for thickness *h* = 10 m and incidence angles 0◦, 30◦ and 60◦. The peaks of the reflectivity are located in the high CO<sup>2</sup> saturation range, since in that case the impedance contrast between the caprock and the sandstone is maximum. The sensitivity to *S<sup>g</sup>* depends on frequency showing bands of very weak reflectivity (in blue scale) and higher amplitudes (in red scale). The high reflectivity zones, which at the same time are sensitive to saturation, are located around the first and second peak frequencies for this example (*h* = 10 m), i.e. near 35 Hz and 105 Hz. These conclusions

high-resolution seismic acquisition is essential for these layers.

The seismic monitoring of carbon dioxide in geologic reservoirs is mostly focused on the characterization of accumulations of high saturation, due to their large seismic amplitudes. Nevertheless, low-saturation zones with dispersed CO2, or *saturation transitions* may have an important role in the propagation of waves within the reservoir, giving rise to amplitude and phase changes of the seismic signals. Our aim is to analyze whether the presence of a saturation transition below the main accumulation substantially alters the previous modeling results. At the same time we investigate the influence of the vertical extent of the transition zone on the reflectivity of the reservoir.

To study these effects, we consider a modification of the model used in Section 3.2, by introducing a variable CO2 saturation-depth profile. This model supposes that the migration of the injected CO2 is controlled by buoyancy and that the reservoir permeability is homogeneous and isotropic, resulting in a laterally constant gas saturation field. The saturation-depth profile is inspired by the axisymmetric flow simulations presented in [9].

For the following examples we represent the saturation *Sg*(*z*) with a piecewise linear function of the vertical depth *z*. The total thickness of the CO2 bearing layer is given by

$$h = h\_0 + h\_t,\tag{9}$$

where *h*0 is the thickness of a constant saturation zone below the caprock, with *Sg* = 100 % and *ht* being the thickness of the saturation transition. The vertical variation in saturation *Sg*(*z*) in turn produces a linear density profile *ρ*(*z*) and a non-linear velocity with depth *Vp*(*z*), as shown in Figure 10. The physical conditions of the CO2 are those described in Section 3.1. To compute the generalized reflection coefficient using the reflectivity method,

**Figure 11.** Normal incidence AVF for transition thicknesses *h<sup>t</sup>* = 0, 5, 10 and 20 m.

the upper halfspace of the model is the caprock, followed by the constant saturation layer with *h*<sup>0</sup> = 10 m, the saturation transition with variable *h<sup>t</sup>* and a fully brine saturated lower halfspace. The transition is discretized in a finite number of layers of 0.5 m with decreasing CO2 saturation with depth.

In Figure 11 a)-d) we compare the normal incidence AVF response for transitions of different thicknesses given by *ht*/*h*<sup>0</sup> = 0, 0.5, 1 and 2, that is *h<sup>t</sup>* = 0, 5, 10 and 20 m. The overall bulk saturation within the transition, given by

$$
\bar{S} = \frac{1}{h\_t} \int\_0^{h\_t} S\_{\mathcal{G}}\left(z\right) dz,\tag{10}
$$

is kept constant at *<sup>S</sup>*¯ = 50 % in the four cases. Due to the layered geometry and variable saturation of the model, the reflectivity vs. frequency shows a pattern of maxima and minima, with decreasing peak amplitudes with frequency. The presence of the saturation transition shifts the first peak to lower frequencies within the seismic band. In general for higher thicknesses *h<sup>t</sup>* we observe lower values of the first peak frequency, which indicates that the inverse proportionality between these parameters still holds. Although the theoretical relation in equation (7) is not applicable to this multilayer geometry, the distinct periodicities can help to discriminate thin and thick transition zones.

The effect of the saturation transition can be noted in more detail in Figure 12, where we plot the first peak frequency *f<sup>p</sup>* vs. thickness *h* for the case with and without saturation transition below the main accumulation. The discrepancies between the curves indicate that the CO2 saturation transition introduces some uncertainty in the estimation of the accumulation thickness using the *f<sup>p</sup>* attribute. In other words, if the existence of a transition is not taken into account, for a given value of *f<sup>p</sup>* the thickness of the accumulation, as given by equation (8), is underestimated leading also to an incorrect estimation of its associated CO2 volume.

A more extended and insightful study about the reflectivity and characterization of CO2 saturation transitions was published by us in [13] and [14]. In those works, a comparative

**Figure 12.** First peak frequency vs. thickness for the accumulation with and without saturation transition. The transition has *h*<sup>0</sup> = 10 m and *h<sup>t</sup>* from 0 to 10 m.

analysis of the results obtained with the layered saturation model and those obtained using the classic Wolf ramp [10], [11] was also presented.

## **4. Conclusions**

16 Carbon Sequestration

CO2 saturation with depth.

saturation within the transition, given by

**Figure 11.** Normal incidence AVF for transition thicknesses *h<sup>t</sup>* = 0, 5, 10 and 20 m.

can help to discriminate thin and thick transition zones.

the upper halfspace of the model is the caprock, followed by the constant saturation layer with *h*<sup>0</sup> = 10 m, the saturation transition with variable *h<sup>t</sup>* and a fully brine saturated lower halfspace. The transition is discretized in a finite number of layers of 0.5 m with decreasing

In Figure 11 a)-d) we compare the normal incidence AVF response for transitions of different thicknesses given by *ht*/*h*<sup>0</sup> = 0, 0.5, 1 and 2, that is *h<sup>t</sup>* = 0, 5, 10 and 20 m. The overall bulk

> *<sup>h</sup><sup>t</sup>* 0

is kept constant at *<sup>S</sup>*¯ = 50 % in the four cases. Due to the layered geometry and variable saturation of the model, the reflectivity vs. frequency shows a pattern of maxima and minima, with decreasing peak amplitudes with frequency. The presence of the saturation transition shifts the first peak to lower frequencies within the seismic band. In general for higher thicknesses *h<sup>t</sup>* we observe lower values of the first peak frequency, which indicates that the inverse proportionality between these parameters still holds. Although the theoretical relation in equation (7) is not applicable to this multilayer geometry, the distinct periodicities

The effect of the saturation transition can be noted in more detail in Figure 12, where we plot the first peak frequency *f<sup>p</sup>* vs. thickness *h* for the case with and without saturation transition below the main accumulation. The discrepancies between the curves indicate that the CO2 saturation transition introduces some uncertainty in the estimation of the accumulation thickness using the *f<sup>p</sup>* attribute. In other words, if the existence of a transition is not taken into account, for a given value of *f<sup>p</sup>* the thickness of the accumulation, as given by equation (8), is underestimated leading also to an incorrect estimation of its associated CO2 volume.

A more extended and insightful study about the reflectivity and characterization of CO2 saturation transitions was published by us in [13] and [14]. In those works, a comparative

*S<sup>g</sup>* (*z*) *dz*, (10)

*<sup>S</sup>*¯ <sup>=</sup> <sup>1</sup> *ht* In this chapter we presented a review on the theoretical seismic reflectivity function for carbon dioxide accumulations. For monitoring purposes, emphasis was placed on the study of its correlation with parameters of interest, such as accumulation thickness and overall saturation. This study is intended to predict quantitatively the changes that can be expected in seismic time lapse data analysis. Although for the examples we calibrated our rock physics models with information of the well known Sleipner field, similar analysis can be implemented for any other geological site, including the case of CO2–oil, or three phase fluid saturation with a proper adaptation of the physical models.

Given that the analysis of seismic AVA is a widely used tool in reservoir geophysics, using the classic Shuey's approximation, we modeled the behavior of intercept, gradient and curvature attributes under variable saturation and physical conditions. From our results we conclude that a monotonic increase in magnitude for the intercept can be expected in time lapse surveys, with significant variations with respect to the pre-injection state. This parameter shows clearly the decrease in the acoustic impedance in the sandstone for increasing CO2 saturation.

The gradient is expected to decay strongly for low CO2 saturations showing a different trend for saturations in the intermediate range, particularly for gaseous conditions. The curvature parameter is not expected to bring much more information about saturation, since in most cases its variations are smaller than those observed in the intercept and gradient. We have also shown that these AVA coefficients do not depend strongly on the choice of the equation of state to model the physical properties of CO<sup>2</sup> under *in-situ* reservoir conditions.

A loss of sensitivity on the AVA parameters, due to the stabilization of the seismic velocities, is characteristic of middle to high CO2 saturations. This imposes a physical limitation on the monitoring of the CO2 saturation. Our results suggest that a combined seismic time lapse analysis of the intercept and gradient AVA coefficients may be useful, not only to constrain the saturation state, but also to distinguish the physical state of the CO2 accumulated below the caprock.

Taking into account that the vertical scale lengths in a layered medium provide physical constraints on the behavior of the propagating waves, we introduced this parameter in the computations of the reflectivity for a better characterization of this attribute. This is an issue of particular importance when the thicknesses of the CO2 accumulations are small compared to the seismic wavelengths, a common situation in geologic repositories. The interference between multiple reflections gives rise to reflectivity dispersion effects and periodicities which can bring useful information. In this context, the use of modern spectral decomposition techniques plays an important role for the study of the amplitudes of reflected waves at different frequencies giving rise to amplitude versus frequency analysis, AVF. As shown in the examples, the periodicity of the reflectivity and the location of the first peak frequency can bring useful information about the thicknesses of the CO2 bearing layers. In this sense, the utilization of the theoretical models presented here can be very useful for interpreting the observations and predicting peak frequencies. Our model estimates also show that the sensitivity of the reflectivity to saturation is frequency dependent, showing that low and high reflectivity bands can be expected. The high reflectivity zones show at the same time significant sensitivity to CO2 saturation, a correlation that can be exploited for monitoring purposes.

We also analyzed the existing correlation between peak frequencies and thickness, to illustrate the utility of this spectral attribute in CO2 accumulations. Our modeling results also indicate that the presence of saturation transitions below a main accumulation is a source of some inaccuracy in thickness and volume determination.

Regarding the potential utilization of these models and attributes (AVA coefficients, AVF, peak frequencies) using real seismic data to monitor changes in the CO2 repository, some precautions must be taken to obtain meaningful estimations:


In this way the theoretical results presented in this work can be tested for real case situations and may become useful tools for carbon dioxide monitoring problems.

## **Acknowledgements**

This work was partially supported by CONICET, Argentina with grants PIP 00952, PIP 00777 and Universidad Nacional de La Plata, Argentina. Some of these results were obtained during the participation of the authors in CO2ReMoVe Project (http://www.co2remove.eu).

### **Author details**

18 Carbon Sequestration

the caprock.

purposes.

analysis of the intercept and gradient AVA coefficients may be useful, not only to constrain the saturation state, but also to distinguish the physical state of the CO2 accumulated below

Taking into account that the vertical scale lengths in a layered medium provide physical constraints on the behavior of the propagating waves, we introduced this parameter in the computations of the reflectivity for a better characterization of this attribute. This is an issue of particular importance when the thicknesses of the CO2 accumulations are small compared to the seismic wavelengths, a common situation in geologic repositories. The interference between multiple reflections gives rise to reflectivity dispersion effects and periodicities which can bring useful information. In this context, the use of modern spectral decomposition techniques plays an important role for the study of the amplitudes of reflected waves at different frequencies giving rise to amplitude versus frequency analysis, AVF. As shown in the examples, the periodicity of the reflectivity and the location of the first peak frequency can bring useful information about the thicknesses of the CO2 bearing layers. In this sense, the utilization of the theoretical models presented here can be very useful for interpreting the observations and predicting peak frequencies. Our model estimates also show that the sensitivity of the reflectivity to saturation is frequency dependent, showing that low and high reflectivity bands can be expected. The high reflectivity zones show at the same time significant sensitivity to CO2 saturation, a correlation that can be exploited for monitoring

We also analyzed the existing correlation between peak frequencies and thickness, to illustrate the utility of this spectral attribute in CO2 accumulations. Our modeling results also indicate that the presence of saturation transitions below a main accumulation is a source of some

Regarding the potential utilization of these models and attributes (AVA coefficients, AVF, peak frequencies) using real seismic data to monitor changes in the CO2 repository, some

• it is very important to use high quality time-lapse seismic data, with good repeatability

• given that the amplitude spectrum of the seismic data is a combination of the wavelet spectrum and the reflectivity of the layers, it is very important for the processing sequence

• correct identification and isolation of a window containing the seismic reflections

• the determination of peak frequencies and their time-lapse variations using an appropriate

In this way the theoretical results presented in this work can be tested for real case situations

This work was partially supported by CONICET, Argentina with grants PIP 00952, PIP 00777 and Universidad Nacional de La Plata, Argentina. Some of these results were obtained during the participation of the authors in CO2ReMoVe Project (http://www.co2remove.eu).

to include amplitude spectral balancing to remove wavelet overprint;

time-frequency decomposition on the selected seismic window is necessary.

associated with the top of CO2 accumulations is necessary;

and may become useful tools for carbon dioxide monitoring problems.

inaccuracy in thickness and volume determination.

and frequency content as wide as possible;

**Acknowledgements**

precautions must be taken to obtain meaningful estimations:

Claudia L. Ravazzoli\* and Julián L. Gómez

\*Address all correspondence to: claudia@fcaglp.unlp.edu.ar

CONICET and Faculty of Astronomical and Geophysical Sciences, National University of La Plata, Argentina

## **References**


[27] Rubino J.G., Ravazzoli C.L., and Santos J.E. Equivalent viscoelastic solids for heterogeneous fluid-saturated porous rocks. Geophysics 2009;74(1):N1-N13.

20 Carbon Sequestration

layer. Geophysics 2010;75(5):A31–A35.

layers. Geophysics 2012;77(3):D75-D83.

of Applied Physics 1962; 33: 1482-1498.

reservoir. Geophysics 1976; 41:882-894.

Phys. of Solids 1963;11:357-372.

1979;44:1777-1788.

Naturforschenden Gessellshaft in Zurich 1951; 96:1-23.

saturated rocks. Geophysics 1998; 63(3):918-924.

flow. Journal of Geophysical Research 2004;109:B01201.

Nacional de La Plata; 2013.

1983;48:148-162.

1992;91:1911-1923.

[11] Liner C.L, and Bodmann B. G.The Wolf ramp: reflection characteristics of a transition

[12] Gómez J. L., Ravazzoli C. L. and Carozzi F. E. Generalized reflectivity of CO2 partially saturated layers. SEG Technical Program Expanded Abstracts 2010;488-492.

[13] Gómez J. L, Ravazzoli C. L. Reflection characteristics of linear carbon dioxide transition

[14] Gómez J. L. and Ravazzoli C. L. Modelling the reflectivity of a carbon dioxide transition zone. Proceedings of the Third EAGE CO2 Geological Storage Workshop 2012;114-118.

[15] Gómez J. L. Seismic monitoring of CO2: reflectivity modeling and attribute analysis. PhD Thesis (in Spanish). Facultad de Ciencias Astronómicas y Geofísicas, Universidad

[16] Dutta N. and Odé H. Seismic reflections from a gas-water contact. Geophysics

[17] Santos J.E., Corberó J., Ravazzoli C. L. and Hensley J. Reflection and transmission coefficients in fluid saturated porous media. Journal of the Acoustical Society of America

[18] Ravazzoli C. L. Analysis of reflection and transmission coefficients in three-phase sandstone reservoirs. Journal of Computational Acoustics 2001;9(4):1437-1454.

[19] Biot M.A. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low frequency range. Journal of the Acoustical Society of America 1956;28:168-171.

[20] Biot M.A. Mechanics of deformation and acoustic propagation in porous media. Journal

[21] Gassmann F. Uber die elastizitat poroser medien. Vierteljahrschrift der

[22] Domenico S.N. Effect of brine-gas mixture on velocity in an unconsolidated sand

[23] Mavko G. and Mukerji T. Bounds on low-frequency seismic velocities in partially

[24] Hill, R. Elastic properties of reinforced solids: some theoretical principles. J. Mech.

[25] Dutta N.C, and Odé H. Attenuation and dispersion of compressional waves in fluid-filled porous rocks with partial gas saturation (White model). Part I: Biot theory Geophysics

[26] Pride S.R., Berryman J.G., and Harris J.M. Seismic attenuation due to wave-induced

