Numerical Simulation of CO 2 Sequestration in Large Saline Aquifers

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.


Introduction
With heightened concerns over CO 2 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 CO 2 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 CO 2 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.

Saline aquifer geological carbon sequestration (SAGCS)
Studies of GCS have suggested that various geological structures can serve as potential CO 2 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 CO 2 (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 CO 2 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 CO 2 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 CO 2 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.

The approach for numerical simulations of GCS
The large spatial extent (of the order of kilometers) and time duration (centuries) for CO 2 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 CO 2 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 CO 2 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".

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  Table 1. Estimation of saline aquifer storage capacity of GCS for different regional carbon sequestration partnerships (RCSPs) [1] 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.

Simulation of in situ CO 2 migration and comparison with analytical solution
As a first step towards code validation, simulations for CO 2 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, b ag (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 CO 2 and brine respectively, and V (t) = ∫ Qdt is the volume of the injected CO 2 within time t. For a horizontal reservoir, setting b ag (r,t) = 0 yields Eq. (2), which gives a quick evaluation of the maximum CO 2 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 hydrogeological properties are used. CO 2 injection rate is set at 1 kg per year for ten years. Detailed model parameters used in the simulations are summarized in Table 1. CO 2 plume migration in each year is computed by the simulation and compared with the analytical solution given by Eq. (2).  The simulation time is ten years and the CO 2 migration within the aquifer is computed for each of the ten years. In Figure 2, the CO 2 plume after 1, 4, 7 and 10 years since the inception of injection is shown.
As seen in Figure 2, the injected CO 2 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 CO 2 migration reaches about 100 m from the injection well. In the following nine years, in situ CO 2 keeps migrating outwards and spreads to 300 m after the tenth year of injection. Physically, such a large radial migration of in situ CO 2 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 CO 2 migration underneath the caprock given by the analytical solution and TOUGH2 simulation As shown in Table 3, TOUGH2 simulations successfully predict the maximum CO 2 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 CO 2 dissolution is accounted for in TOUGH2 while it is neglected in the derivation of the analytical solution.
Since CO 2 dissolution is governed by the contact area between CO 2 and the ambient brine, the rate of CO 2 dissolution into brine should gradually increase over time as larger contact area becomes available as the CO 2 plume spreads. Nevertheless, Table 3 validates TOUGH2 as an accurate simulation tool for predicting the migration of in situ CO 2 .
A schematic of the CO 2 plume flow is shown in Figure 3. Although based on a simplified analytical model, Figure 3 shows in situ CO 2 migration due to the combined pressure-driven Darcy flow and the buoyancy-drive CO 2 transport.With the injection well located on the left side of Figure 3, the CO 2 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, CO 2 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 CO 2 and Darcy flow occurs, causing CO 2 to migrate more radially through the aquifer. The second region is marked as region (2) in Figure 3, where the CO 2 plume fully develops. In this region, buoyancy due to the density difference between CO 2 and brine becomes dominant and drives the upward movement of CO 2 along with lateral migration. In this region, the vertical movement of CO 2 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.

Simulation of benchmark problem #1 -CO 2 plume evolution and leakage through an abandoned well
A three-layered formation is modeled for the first benchmark problem [9]. CO 2 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 CO 2 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. 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 CO 2 , 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 CO 2 injection well. It can be either a crack in the formation or a physical abandoned well, which served as a pathway for upward CO 2 migration. Supercritical CO 2 is injected only into the lower aquifer through the injection well. Being less dense than brine, injected supercritical CO 2 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.   To accurately model the small dimensions of the wells and accurately capture the CO 2 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 CO 2 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. Other simulation parameters such as initial conditions and boundary conditions are summarized in Table 6.

Thermal Condition Isothermal
Initial 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 CO 2 to achieve equilibrium condition under gravity. The equilibrium state is then implemented as an initial condition for the subsequent simulation with CO 2 injection. The equilibrium simulation is critical to provide the simulation with CO 2 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 CO 2 saturation distribution throughout the aquifer after 80 days of CO 2 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 CO 2 leakage rate to CO 2 injection rate. Detailed comparisons using various simulation codes are shown in Figure 6 and are summarized in Table 7.  Table 7. Simulation results and comparisons with other codes for benchmark problem #1 As additional comparisons, the pressure perturbation and CO 2 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. As seen in Figure 8, the CO 2 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 CO 2 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".

Simulation of benchmark problem #2 -enhanced CH 4 recovery in combination with CO 2 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 CO 2 into the mature reservoir. With void space being occupied by the injected CO 2 , remaining oil/gas is pushed out of the reservoir. Meanwhile, the depleted reservoir becomes an ideal carbon sink for longterm storage. EOR/EGR with CO 2 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 CO 2 is injected in the middle of the domain at the bottom-most part. This is a direct result of CO 2 being heavier than CH 4 under the reservoir conditions. Injection of CO 2 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 CH 4 production to the original remaining CH 4 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 CO 2 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.  Initial conditions and boundary conditions and some other parameters of the domain are given in Table 9.

Thermal Condition Isothermal
Initial  Table 9. Simulation parameters for benchmark problem #2 The termination of the simulation depends solely on the mass fraction of CO 2 in the reservoir. It takes about 30 minutes of CPU time to run 2,000 days of simulation before major CO 2 contamination occurs. The recovery factor, production shut-down time, pressure and CO 2 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 comparisons for recovery factor and production well shut-down time with other investigators' simulations [12]. To visualize how the displacement process of CH 4 by CO 2 works, the density and CO 2 mass fraction profiles at production shut-down are examined in Figure 12.  by upward movement of CO 2 occurs at the production well despite the gravity segregation. This is due to the strong convective flow near the production well and mass diffusion.

Simulation of benchmark problem #3 -CO 2 injection in a heterogeneous geological formation
Accurate estimation of in situ CO 2 dissolution into the ambient brine is another important aspect of SAGCS simulation, since CO 2 becomes securely sequestrated once dissolved.
Overestimation of CO 2 dissolution could lead to underestimating the potential leakage; on the other hand, underestimation of CO 2 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 CO 2 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. CO 2 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 CO 2 , the amount of CO 2 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 CO 2 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. The porosity and permeability of the modeled Johanson formation areshown in Figure 15.  Table 11. Both gaseous and aqueous CO 2 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.  A comparison of the CO 2 migration after 50 years is given in Figure 17. 3) CO 2 dissolution into the ambient porewater is a process that takes place very slowly. 4) The greater slope of aqueous CO 2 during the first 25 years (when injection continues) implies enhanced carbon dissolution due to convection during the first 25 years.

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

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

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 CO 2 [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. 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 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 CO 2 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 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. Due to the relatively high porosity and permeability, CO 2 injection is assigned at the bottom Arkosic unit (bottom three sublayers). The injection rate is assigned to be 5 million tons per  [20] year and injection lasts for 50 years. CO 2 footprint at 5 th year, 25 th year and 50 th year since the beginning of injection is examined. As seen in Figure 20, the CO 2 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 CO 2 is observed. In the overlying sublayers, however, a strong secondary sealing effect that retards the vertical migration of gaseous CO 2 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 CO 2 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 CO 2 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 CO 2 , and thus a stronger vertical pressure gradient is required for mobile CO 2 to penetrate sublayer #20. When the capillary pressure is greater than the phase pressure of CO 2 , sublayer #20 appears to be "impermeable" to the underlying CO 2 plume. Consequently, gaseous CO 2 accumulates under this layer and continues spreading out laterally, finally reaching a maximum extent of approximately 3000 m. Meanwhile, the increased CO 2 column under sublayer #20 brings up its phase pressure. Once the phase pressure of CO 2 exceeds the entry pressure of sublayer #20, mobile CO 2 breaks the capillary barrier of its overlying layer and starts to penetrate it. Such accumulation-penetration-breakthrough behavior of gaseous CO 2 occurs each time the upward migrating CO 2 encounters an overlying sublayer with lower permeability. Because the high capillary entry pressure of the overlying layer temporarily prevents CO 2 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 CO 2 . Its contribution makes gaseous CO 2 barely reach the Eau Clair shale even after 50 years of injection.

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 CO 2 migration and dissolution for fullterm SAGCS simulation. A full-term simulation refers to a simulation that investigates the fate of in situ CO 2 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 CO 2 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°f rom 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. 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 CO 2 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.
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 CO 2 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 CO 2 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.
As seen from Figure 23, the highly asymmetric CO 2 plume suggests a strong tendency of movement upwards in the direction of the geological slope. Unlike the case of Mt. Simon SAGCS, the CO 2 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 determining the in situ CO 2 migration. Considering the relatively short-term injection (10 days) for the Frio SAGCS project, this implies that in situ CO 2 migrates mostly convectively. Furthermore, it demonstrates that the poor permeability caprock layer above the injection serves quite well as the CO 2 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 CO 2 saturation, etc. (2) Discrepancy still exists at the detailed simulation level. The results show that in our simulations, CO 2 saturation at the injection well reaches a maximum of 0.8 by the 10 th 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 CO 2 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 CO 2 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 CO 2 -brine interface. This also explains the slight overshoot in the plume migration to the north in our simulation.

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 CO 2 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. CO 2 is captured from the gas mixture produced from a nearby deeper natural gas reservoir. To date, approximately 1 million tons of supercritical CO 2 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 CO 2 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 secondary sealing effect will occur. Figure 25 shows a 2-D seismic image taken in 2008 revealing the CO 2 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 CO 2 . 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 km 2 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 CO 2 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 CO 2 migration, while avoiding intensive computational effort associated with full 3D modeling and simulation of the entire Utsira formation.

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 CO 2 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 37 o 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. CO 2 injection of 30 kg/s is assigned as a point source at the middle of the bottom-most sand layer. The simulation time is set at 15 years and the CO 2 plume profile is examined for each year. Figure 27 shows the cross-sectional view of gaseous CO 2 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 CO 2 . Similar to the case of Mt. Simon SAGCS, the injected CO 2 first migrates upwards driven by buoyancy until it reaches the first shale layer. Due to the low permeability and high capillary entry pressure, CO 2 is confined by this shale layer and is forced to migrate radially. Simultaneously, CO 2 concentration builds up beneath the shale layer and finally breaks through the capillary barrier upon sustaining sufficient CO 2 column height. The accumulation-penetration-breakthrough takes place each time the CO 2 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 CO 2 has very limited contact with the caprock of the Utsira formation by the third year of injection. Additionally, ten-year CO 2 flux analysis has been made for the topmost sandstone layer (Layer #9) since it is critical to identify the accumulation of CO 2 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 CO 2 in the topmost sandstone layer tends to increase until it becomes stabilized. This fact can be explained by mechanism of secondary sealing effect.

Model #2 -Detailed 3D model of the Utsira Layer #9 formation
In situ CO 2 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 CO 2 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 CO 2 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 CO 2 and has direct contact with the overlying caprock formation. Seismic surveys have shown striking growth of CO 2 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 CO 2 accumulations appeared after about three years of injection (recall that injection began in 1996), indicating CO 2 began to accumulate under the caprock. However, CO 2 migration in Layer #9 was not symmetric due to the topography of the caprock. The northward migration of initially impacted CO 2 , 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 CO 2 along a north-tending topographic ridge. CO 2 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, CO 2 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 CO 2 , the quantification of CO 2 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 CO 2 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 CO 2 , 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.   It should be noted that in the 3D Layer #9 model, the source of CO 2 is identified as "feeder" but not "injector" to emphasize that CO 2 is supplied from the lower aquifer through leakage pathways rather than by direct injection. Since the actual CO 2 injector is located at about 200 m under Layer #9, information on the injection rate recorded at the injector is not applicable for the CO 2 feeders in the Layer #9 model. To determine the CO 2 feeding rate to Layer #9, seismic surveys of CO 2 distribution are used to obtain its volume under in situ conditions, and then converted to mass flow rate. Information of CO 2 accumulative mass provided by Zhu and Lu [28] is summarized inTable 17.  Table 17 gives the CO 2 accumulation in Layer #9 from 1999 to 2008. It can be seen that CO 2 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 CO 2 phase pressure at lower aquifer and the capillary pressure of the overlying shale layer that determines the breakthrough of CO 2 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 CO 2 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 CO 2 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) CO 2 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 CO 2 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. CO 2 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. 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 CO 2 . Horizontal anisotropy of 2:10 is needed to capture the northern spill of CO 2 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 CO 2 breakthrough from the lower aquifer structure. Third, the fact that the injection gas is a CO 2 -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 CO 2 be modeled. This is consistent with the behavior of CO 2 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 CO 2 under the caprock (shale) without major leakage, implying that the caprock serves quite well as the non-permeable CO 2 barrier while exerting little resistance for the lateral flow of CO 2 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.

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.