Open access peer-reviewed chapter

Numerical Simulation of CO2 Sequestration in Large Saline Aquifers

By Zheming Zhang and Ramesh K. Agarwal

Submitted: September 18th 2013Reviewed: September 19th 2013Published: March 12th 2014

DOI: 10.5772/57065

Downloaded: 1639

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 permeability 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 computations 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.

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) conventional 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 sequestration in saline aquifers considering the following facts.

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

  • Geological survey has confirmed vast geological distribution of deep saline formations possibly suitable for GCS in the US and Canada.

  • 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 of US overall GCS potential [1].

  • Since most of the saline formations are located deep underground, i.e., at least 800 m below sea level, they provide great potential for secure long-term sequestration.

  • Many surveys, research projects and commercial projects have already been conducted for SAGCS, making it attractive for further research and technical contributions.

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.

RCSPsLow-end capacity(billion metric tons)High-end capacity(billion metric tons)
BSCSP (Big Sky Caron Sequestration Partnership)981237
MGCS (Midwest Geological Sequestration Consortium)11158
MRCSP (Midwest Regional Carbon Sequestration Partnership)95123
PCOR (Plains CO2 Reduction Partnership)174511
SECARB (Southeast Regional Carbon Sequestration Partnership)137614089
SWP (Southwest Regional Partnership on Carbon Sequestration)2662801
WESTCARB (West Coast Regional Carbon Sequestration Partnership)821124

Table 1.

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

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 simulations using computation fluid dynamics (CFD) technology can be employed at industrial 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 developed 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 worldwide. 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:


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)=Qdtis 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:


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 hydro-geological 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

Geometry100 m in thickness; 3000 m in radius
PermeabilityIsotropic, 100 mDarcy
Pressure6.4 MPa
CO2 Density789.96 kg/m3
CO2 Viscosity0.0000712905 Pa∙s
Brine Density1029.69 kg/m3
Brine Viscosity0.001488427 Pa∙s
Relative PermeabilityLinear
Brine Residual Saturation0
CO2 Residual Saturation0
Capillary PressureNone
Injection Rate1 kg/s
Boundary ConditionOpen boundary
Domain Discretization300 × 20

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.

Figure 2.

CO2 plume after 1, 4, 7 and 10 years after injection

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

YearMaximum Migration based on Numerical Simulation (A)Maximum Migration based on Analytical Solution (B)Deviation((A-B)/B)
1100.75 m95.58 m0.054090814
2140.49 m135.17 m0.039357846
3168.37 m165.55 m0.017034129
4191.34 m191.16 m0.00094162
5211.27 m213.72 m-0.011463597
6229.23 m234.12 m-0.020886725
7235.34 m252.88 m-0.069360962
8260.81 m270.34 m-0.035251905
9275.70 m286.74 m-0.038501779
10289.36 m302.25 m-0.042646816

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

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

Domain Dimension1000 m × 1000 m × 160 m
Aquifer Depth2840 m ~ 3000 m
Aquifer Thickness30 m
Aquitard Thickness100 m
Distance between Injection Well and Leaky Well100 m
Injection & Leaky Well Geometry0.266 m × 0.266 m

Table 4.

Geometry parameters for benchmark problem #1

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.

Aquifer Permeability20 mDarcy
Leaky Path Permeability1 Darcy
Residual Brine Saturation0.2
Residual CO2 Saturation0.05
Relative PermeabilityLinear (kr = S)
Capillary PressureBrooks-Corey
Entry Pressure1.0 × 105 Pa
Brooks-Corey Parameter2.0

Table 5.

Hydrogeological parameters for benchmark problem #1

Other simulation parameters such as initial conditions and boundary conditions are summarized in Table 6.

Thermal ConditionIsothermal
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
Boundary ConditionsFixed-state on lateral boundaries
No mass flow on top and bottom boundaries
Initial CO2 Mass FractionXCO2 = 0
Initial Salt Mass FractionXsm = 0.20
Injection Rate8.87 kg/s
Simulation End Time1000 days

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 implemented as an initial condition for the subsequent simulation with CO2 injection. The equilibrium 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

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.

Max. LeakageTime at Max. LeakageLeakage at 1000th Day
TOUGH2 (WUSTL)0.225 %100 days0.115 %
TOUGH2 (BRGM)0.226 %93 days0.110 %

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 long-term 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.

Figure 9.

Schematic of the 3D five-spot pattern for benchmark problem #2

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

Figure 10.

CFD model and its mesh in the quarter of the 5-spot domain for benchmark problem #2

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

PermeabilityHorizontal: 50 mDarcy, Vertical: 5 mDarcy
Residual Brine Saturation0
Relative PermeabilityLiquid: immobile, Gas: linear (kr = S)
Capillary PressureNone

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.

Thermal ConditionIsothermal
Initial Condition (Temperature)66.7 oC
Initial Condition (Pressure)3.55 MPa
Boundary ConditionsNo mass flow at all boundaries; Constant pressure at CH4 production well
Initial CO2 Mass FractionXCO2 = 0
Initial CH4 Mass FractionXsm = 1
Injection Rate0.1 kg/s until shut-down

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.

Figure 11.

History of enhanced CH4 recovery for benchmark problem #2

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 comparisons for recovery factor and production well shut-down time with other investigators’ simulations [12].

Recovery FactorProduction Well Shut-down Time
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 heterogeneous 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

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.

Residual Brine Saturation0.2
Residual CO2 Saturation0.05
Relative PermeabilityBrooks-Corey
Capillary PressureBrooks-Corey
Entry Pressure1.0×104 Pa
Brooks-Corey Parameter2.0

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.

Thermal ConditionIsothermal
Initial Temperature0.03 oC/m; 100 oC at 3000 m depth
Initial Pressure1075 Pa/m, 30.86 MPa at 3000 m depth
Boundary ConditionsFixed-state on lateral boundaries
No mass flow on fault, top and bottom boundaries
Initial CO2 Mass FractionXCO2 = 0
Initial Salt Mass FractionXsm = 0.1
Injection Rate15 kg/s (for 25 years), 0 kg/s thereafter
DiscretizationNumber of computational grids: 18804
Non-uniform for x-, y-, z-directions

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

Gashouse CO2 at 50th YearDissolved CO2 at 50th Year
TOUGH2 (WUSTL)87.9%12.1%
TOUGH2 (CO2/CRC)86.5%13.5%
IPARS (U. Texas)79.1%20.9%

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 porewater 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 convection during the first 25 years.

4.5. Conclusions of benchmark problems

The benchmark simulations presented in this chapter demonstrate that the TOUGH2 numerical 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 efficiency, 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 significantly.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. Low-permeable 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.

Figure 18.

Core injection area and elevation of Mt.Simon sandstone

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

Sub-UnitSublayerLayer Depth (m)Mean PorosityMean Permeability (mDarcy)Characteristic Capillary Pressure (bar)
Upper Unit12140 – 21500.06150.37
22150 – 21820.1093000.06
32182 – 21970.074100.28
42197 – 22030.0833.60.4875
52203 – 22300.1951100.1
62230 – 22320.0711.10.8
72232 – 22800.132100.083
Middle Unit82280 – 23220.0835.40.4125
92322 – 23310.241500.0875
102331 – 23400.08880.35
112340 – 23500.1568000.095
122350 – 23700.25800.125
132370 – 23780.1639000.095
142378 – 23850.1951050.1007
152385 – 23990.1638000.05
162399 – 24060.136720.1167
172406 – 24120.1567000.05
182412 – 24240.1291600.09
192424 – 24300.1618500.05
202430 – 24620.128600.15
Arkosic Unit212462 – 25000.20210000.05
222500 – 25020.141900.09
232502 – 25370.15110000.04

Table 14.

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

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

Figure 19.

(a) Permeability, (b) porosity, and computational mesh of the 24 sublayers of the Mt.Simon formation model 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 secondary 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 maximum 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 permeability. 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

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 full-term 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.

PermeabilityIsotropic, 1 mDarcy
Residual Brine Saturation0.15
Residual CO2 Saturation0.2
Relative Permeabilityvan Genuchten-Mualem with hysteresis
Capillary Pressurevan Genuchten-Mualem with hysteresis
Initial ConditionsP = 15.2 MPa, T = 59oC
Initial CO2 Mass FractionXCO2 = 0
Initial Salt Mass FractionXsm = 0.093

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.

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.

Figure 22.

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

The simulation time is set at 10 days to match the actual duration of injection. It takes approximately 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.

Figure 23.

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

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

CO2 saturation profiles given by the numerical simulations and the RST logs

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 pyramid-shaped plume indicate strong evidence of the dominant role of gravity segregation in determining 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. Furthermore, it demonstrates that the poor permeability caprock layer above the injection serves quite well as the CO2 barrier.

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 dry-out 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

The Sleipner project near the Norwegian coast on the North Sea is probably the most prestigious, 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 high-permeability sandstone and low-permeability shale. Therefore, it is expected that the secondary 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.

Figure 25.

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

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 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 37oC 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.

Figure 28.

Gaseous CO2 accumulation in the topmost sandstone layer

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 experience 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 accumulation of CO2 and its tendency for migration underneath the caprock. With such information 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 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].

Figure 29.

Amplitude maps of Layer #9 from 1999 to 2008, from Singh et al. [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 #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.

Temperature33 oC
Pressure8.6 MPa
Total Utsira Formation Area26100 km2
Total Utsira Formation Thickness50 m ~ 300 m
Layer#9 Area1600 m × 4900 m
Layer#9 Thickness3.5 m ~ 26.3 m
Shale PermeabilityHorizontal: 0.001 mDarcy, Vertical: 0.0001 mDarcy
Mudstone PermeabilityW-E: 2 mDarcy, N-S: 10 mDarcy, Vertical: 200 mDarcy
Utsira Porosity (shale/mudstone)35.7 %
Residual CO2 Saturation0.02
Residual Brine Saturation0.11
Relative Permeability TypeCorey/van Genuchten-Muller
Capillary PressureNone
Porewater Salinity3.3 %
North-south Geological Slope8.2 m/km, 5.8 m/km
CO2 Feeder LocationMain feeder: W-E: 516 m, N-S: 1210 m, bottom mudstone
Secondary feeder: W-E: 925 m, N-S: 2250 m, bottom mudstone

Table 16.

Hydrogeological properties of the Utsira Layer #9 model

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.

YearAccumulated Mass (kg)Yearly Feeding Mass (kg)Feeding Rate (kg/s)

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 prominent 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 history-matching 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 secondary-sealing 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

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.


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

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

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Zheming Zhang and Ramesh K. Agarwal (March 12th 2014). Numerical Simulation of CO2 Sequestration in Large Saline Aquifers, CO2 Sequestration and Valorization, Claudia do Rosario Vaz Morgado and Victor Paulo Pecanha Esteves, IntechOpen, DOI: 10.5772/57065. Available from:

chapter statistics

1639total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Seismic Reflectivity in Carbon Dioxide Accumulations: A Review

By Claudia L. Ravazzoli and Julián L. Gómez

Related Book

First chapter

Urban Air Pollution

By Bang Quoc Ho

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us