Grantee Research Project Results
2010 Progress Report: Understanding and Managing Risks Posed by Brines Containing Dissolved Carbon Dioxide
EPA Grant Number: R834383Title: Understanding and Managing Risks Posed by Brines Containing Dissolved Carbon Dioxide
Investigators: Falta, Ronald W. , Benson, Sally M. , Murdoch, Lawrence C.
Current Investigators: Falta, Ronald W. , Murdoch, Lawrence C. , Benson, Sally M.
Institution: Clemson University , Stanford University
EPA Project Officer: Aja, Hayley
Project Period: November 1, 2009 through October 31, 2012 (Extended to October 31, 2014)
Project Period Covered by this Report: January 1, 2010 through December 31,2010
Project Amount: $891,342
RFA: Integrated Design, Modeling, and Monitoring of Geologic Sequestration of Anthropogenic Carbon Dioxide to Safeguard Sources of Drinking Water (2009) RFA Text | Recipients Lists
Research Category: Targeted Research , Water
Objective:
Geologic disposal of supercritical carbon dioxide in saline aquifers and depleted oil and gas fields will cause large volumes of brine to become saturated with dissolved CO2 at concentrations of 50 g/l or more. As CO2 dissolves in brine, the brine density increases slightly. This property favors the long-term storage security of the CO2 because the denser brine is less likely to move upwards towards shallower depths. In fact, one proposed strategy for reducing risk from CO2 injection activities involves pre-dissolving the CO2 into brine at the surface, and injecting this brine into the disposal formation. While dissolved phase CO2 poses less of a threat to the security of shallower drinking water supplies, the risk is not zero. There are plausible mechanisms by which the CO2-laden brine could be transported to a shallower depth, where the CO2 would come out of solution (exsolve), forming a mobile CO2 gas phase. This significant mechanism for drinking water contamination has received little attention, and there are basic science and reservoir engineering questions that need to be addressed in order to reduce risks to underground drinking water supplies.
Six main activities were identified in the research proposal:
1) Laboratory Experiments. Laboratory core experiments that flood cores with CO2-saturated brines at reservoir pressure and temperature. These cores then are gradually depressurized, and imaged using a medical CT scanner to study the CO2 phase evolution and movement.
2) Pore-Scale and Core-Scale Modeling. Pore-scale multiphysics and multiphase continuum modeling of these experiments are used to develop a fundamental understanding of the exsolution and CO2 bubble coalescence phenomena as the CO2 starts to form a mobile phase.
3) CO2 Phase Relative Permeability Functions for Multiphase Flow Models. Core-scale multiphase continuum modeling to upscale the experimental results with a focus on developing effective relative permeability functions for use in field-scale modeling.
4) Regional-Scale Variable Density Groundwater Modeling. Simulations of regional scale behavior using a variable density groundwater flow model. The simulations are designed to evaluate the effects of upward hydraulic gradients, upward pressure driven brine flow (for CO2 brine injection projects) and most importantly, the effects of groundwater pumping from shallower aquifers.
5) Multiphase Flow Simulations of Field-Scale CO2 Injection. Local field-scale (hundreds to thousands of meters) multiphase simulations of the likely failure modes using realistic hydrogeologic and geologic conditions that are representative of CO2 storage locations.
6) Remediation Designs. These models will then be used to study remediation strategies and alternative storage methods for each CO2 release scenario.
Progress Summary:
During our first year on this project, we have made very good progress on Activities 1, 3, and 5. We have begun work on Activities 2 and 4, but have not yet started on Activity 6. A detailed summary of progress by activity area is provided in this report.
1. Laboratory Experiments
The laboratory apparatus is up and running, and several exsolution experiments have been conducted on Berea sandstones. Two different approaches to these measurements have been used. Both demonstrate the emergence of gaseous CO2 (exsolved CO2) when the pressure decreases. The second method also allows measurement of the relative permeability to CO2 and water in the presence of exsolved CO2.
Approach One—Rapid Depressurization: Core samples were initially filled with CO2 saturated water at 50° C and 1800 psi. After this, fluids were extracted as fast as possible (~100 psi/min) until the pressure dropped to a specified value. The core holder was then sealed and the system was monitored for 6 to 24 hours to observe equilibration of exsolved CO2 under the influence of gravity and capillary forces. This experiment was conducted using both homogeneous and heterogeneous samples of Berea sandstone.
Evolution of significant amounts of exsolved CO2 was observed in these experiments. CO2 saturations in core samples were 3%~5% at 1400 psi, 10%~15% at 1000 psi, 20%~27% at 800 psi, and 40%~46% at 400 psi. Exsolved CO2 had uniform distributions in homogeneous core samples but showed more complexity in heterogeneous core samples due to permeability variations and flow preference, shown in Figure 1.1. No statistically significant correlation between porosity and CO2 saturation was observed either before or after equilibration. No obvious upward migration of CO2 occurred during the pressure drop, nor during the equilibration period.
Multiphase flow simulation results using TOUGH2-ECO2N, based on measured core sample porosity, using a capillary curve and CO2/water relative permeability curve from a traditional core flooding experiment are shown in Figure 1.2. A reasonable, though imperfect match between the simulated and measured saturation distribution was observed after the pressure was dropped to 400 psi.
Although these experiments clearly demonstrated the exsolution of CO2, they provided little information about the mobility of exsolved CO2.
Figure 1.1. CO2 distributions compared with porosity distributions of homogeneous rock (left) and heterogeneous rock (right) at 400 psi. |
Figure 1.2. CO2 distribution in the vertical direction of the middle rock slice compared with numerical simulation result 6 hours after initial pressure drop. |
Approach Two—Slow Depressurization: These experiments were designed to provide more information about the mobility (relative permeability) of exsolved CO2 and water. Highly permeable Berea sandstone was used in these experiments. The core sample initially was filled with CO2 saturated water at 50° C and 1600 psi. The pressure was raised to 1800 psi. Following this, fluids were extracted from one end of the core holder at a sequence of constant withdrawal rates until pressure dropped to a specified value. With CT scanning and pressure drop measurements across the core, we were available to calculate relative permeabilities of CO2 and water based on mass balance. Equilibrations after pressure drops also were allowed to observe CO2 re-distribution in the core.
With much slower pressure drops compared to the first part, CO2 saturations in core samples were 3%~6% at 1200 psi, 6.5%~10% at 1000 psi, 13%~17% at 800 psi, and 41%~43% at 400 psi. The critical gas saturation, when the CO2 phase became mobile, was 11.7%~15.5% at core sample pressure 815 psi to 845 psi. This result has very significant implications for the security of dissolved CO2, because it shows that if the brine is depressurized, and the CO2 exsolves, that it does not become mobile until a high CO2 phase saturation is reached. This is much different than the behavior of separate phase CO2 movement when it is injected into a water filled rock. In that case, the CO2 phase can be mobile at much lower phase saturations because it preferentially seeks out the larger pore spaces in the porous media.
It is confirmed by repetitive scans that there is no correlation between CO2 saturation and porosity. Re-dissolution, rather than vertical re-distribution, appeared as equilibrium proceeded, shown in Figure 1.3.
Figure 1.3. Vertical CO2 distributions of the middle rock slice as equilibrium proceeded. |
The measured relative permeability curves (Figure 1.4) show the relative permeability of water was reduced to less than 0.1 once exsolved CO2 appeared above a few percent saturation. More importantly, the relative permeability of the CO2 phase was very low (10-5~10-3), even when exsolved CO2 saturation increased to more than 40%. This result implies that even if CO2 exsolves from a depressurized brine, the resulting CO2 phase will have extremely low mobility in the formation. The significant reduction in both water and CO2 mobility observed in these experiments may form a kind of permeability barrier, which prevents CO2’s migration or even blocks possible leakage paths.
Figure 1.4 shows two repetitive relative permeability curve measurements for Berea sandstone. Poor repetition exists between two experiments, especially for CO2 curves. The measurement errors and the artifacts close to critical gas saturations were introduced primarily by measuring small pressure differential at low flow rates in the highly permeable core sample (966 mD). Further experiments will be conducted with core samples of lower permeability and at higher flow rates to reduce errors in pressure measurements.
Figure 1.4. Relative permeability curves of CO2 and water. |
2. CO2 Phase Relative Permeability Functions for Multiphase Flow Models
Our work in this area so far has focused on the development and testing of new hysteretic relative permeability and capillary pressure functions for cases where supercritical or gas phase CO2 is moving through brine-filled porous media. As we have seen from the recent laboratory experiments, gas phase or supercritical CO2 that exsolves from a brine containing dissolved CO2 may behave much differently, and this will be a focus of our work during the coming year.
Characteristic curves that govern the interactions of immiscible phases in numerical flow simulators for environmental and petroleum applications often are used to model the flow of supercritical carbon dioxide during geologic carbon sequestration. These curves describe relative permeability and capillary pressure as functions of phase saturation in a multiphase simulation. Commonly used non-hysteretic characteristic curves follow only one path for wetting-phase drainage and imbibitions processes and are limited to using a single value for residual nonwetting phase saturation. Hysteretic characteristic curves take advantage of variable values of residual nonwetting phase saturation dependent on the history of the phase saturation in the respective grid-block. Hysteresis is critical to predicting the amount of trapped fluid phase when both wetting-phase drainage and imbibitions processes occur along a wide range of phase saturation values; such as when injected supercritical CO2 is subject to buoyant flow in a deep saline aquifer. We have developed a simple hysteretic model for characteristic curves that can be applied to environmental or carbon sequestration numerical modeling.
Multiphase flow numerical models, those used in the petroleum and environmental industries and those used for CO2 sequestration, use characteristic curves to describe the interactions of the separate phases. Characteristic curves are functions that describe capillary pressure–saturation and relative permeability–saturation relationships at the grid-block scale, and ultimately control the multiphase flow. Characteristic curves may differ between grid-blocks depending on the properties of the material each grid-block represents.
In a system consisting of two immiscible fluids in a porous medium, there exists an interface that distinguishes the two fluids as separate phases. This interface can shift through the pores as a response to pressure differences between the two phases, or capillary pressure (Pc). Increases in Pc can drive the nonwetting phase (the phase that does not preferentially adhere to solid surfaces) to advance into the pores, pushing the wetting-phase (the phase that preferentially adheres to solid surfaces) out while decreases in Pc reverse the process.
Figure 2.1 shows a typical non-hysteretic capillary pressure curve. As an increase in capillary pressure occurs, the model follows the branch upwards (wetting-phase drainage process) and wetting-phase saturation (Sw) can decrease until the irreducible wetting-phase saturation (Slr) is reached (signifying a complete wetting-phase drainage). A subsequent decrease in capillary pressure will cause the model to follow the branch back downwards (wetting-phase imbibition) as Sw increases. As long as the nonwetting phase saturation (Sn) reached at least the user-set value of Snr during wetting-phase drainage, the wetting-phase imbibition process can continue only until the Snr is reached. Otherwise, Sw can continue to increase to 1, leaving no residual nonwetting phase in the grid-block.
Figure 2.1. van Genuchten [1980] non-hysteretic capillary pressure curve.
A relative permeability curve is a function that controls the ease at which fluids may flow in a multiphase model. Relative permeability (kr) is a scaling factor (0–1) used to adapt the Darcy equation to multiphase flow (Equation 2.1). To clarify, a porous medium will have a specific conductivity value for a given fluid, e.g., water under single-phase conditions. The presence of a second fluid-phase, e.g., NAPL (non-aqueous phase liquid), however, in the pores of the medium will decrease the ease at which the water flows through the medium. Therefore, it is necessary to scale the fluid conductivity of the medium for each fluid phase to represent the competition of multiple phases for pore space fluid pathways.
Equation 2.1 shows the Darcy equation for multiphase flow where qβ is the Darcy flux for phase β, , krβ is the relative permeability for phase β, μβ is the viscosity of phase β, and Pβ represents the pressure gradient within phase β. Equation 2.2 shows the relative permeability of phase β being equal to the ratio of the effective permeability of phase β (keβ) to the single-phase permeability of phase β.
Figure 2.2 shows a typical set of non-hysteretic relative permeability curves. As wetting-phase drainage occurs, the relative permeability of the wetting-phase (krw) decreases while the relative permeability for the nonwetting phase (krn) increases. The nonwetting phase relative permeability curve in this figure has a residual saturation (Snr) equal to zero, while the wetting phase curve has a residual saturation (Slr) of 0.2. With this nonwetting phase relative permeability curve, there is no trapping of the nonwetting phase, regardless of the flow history. Similar nonwetting phase relative permeability curves have been devised where a single fixed value of the nonwetting residual saturation is used. However, these tend to underpredict the nonwetting phase relative permeability when the phase first invades a water filled system, and they can over predict the mobility if the nonwetting phase is being displaced by water.
|
Figure 2.2. Mualem [1976] non-hysteretic 2-phase relative permeability curve.
The characteristic curves in Figure 2.3 are hysteretic, meaning they describe the capillary pressure and relative permeability based not only on the current saturation but also on the history of saturation at the respective grid-block. These curves attempt to recreate the physical phenomenon known as hysteresis, where the amount of trapped fluid in the pores of a medium is dependent on the maximum amount of that fluid that initially entered the pores of the medium. Note the multiple wetting-phase imbibition branches for different
|
|
Figure 2.3. Hysteretic characteristic curves from Kaluarachchi (1992). Notice from the multiple imbibition branches the value of Snr depends on the saturation at the turning point.
turning-point saturations (the saturation value at the point where the process changes from wetting-phase drainage to wetting-phase imbibition). The most important aspect of hysteretic characteristic curves is the calculation of Snr—represented as the value of Sn at which a wetting-phase imbibition branch has led to a Pc or krn value of 0—which varies depending on the maximum value of Sn reached during the wetting-phase drainage process. Greater values of Sn reached during wetting-phase drainage correspond to greater values of Snr.
The amount of fluid trapped at the pore scale affects the overall mobility of the fluid in the subsurface. When there exists both drainage and imbibition processes, e.g., during buoyancy-driven flow of supercritical CO2 through geologic media, the use of hysteretic values for residual saturation is critical to the behavior of the model. Hysteresis has been observed experimentally in problems involving NAPLs (e.g., Lenhard and Parker, 1987; Lenhard, 1992; Kueper, et al., 1993; Steffy, et al., 1997; Johnston and Adamski, 2005) and hysteretic characteristic curves have been developed and used for modeling in the petroleum and environmental industries.
A disadvantage of the hysteretic functions shown in Figure 2.3 is that the derivatives of the relative permeability and the capillary pressure are discontinuous at the turning points. In multiphase flow codes such as TOUGH2-ECO2N, this leads to numerical difficulties, and the curves must be artificially smoothed as the turning points are approached. As an alternative, we have developed a simple technique for incorporating the key hysteretic effects in relative permeability and capillary pressure curves. The resulting curves are always smooth, and continuously differentiable, leading to good numerical performance in TOUGH2-ECO2N.
Numerical Methods
The new method was used with a nonwetting phase relative permeability function (Equation 2.3) recently proposed by Charbeneau (2007). Effects of hysteresis are included through the calculation of Snr as a linear function (Equation 2.5) of the maximum nonwetting phase saturation (Snmax) in each grid-block. Experiment data for a NAPL-water system from Johnston and Adamski (2005) show residual nonwetting phase saturation to have a linear relationship with maximum nonwetting phase saturations Snmax. The slope of the linear function (ƒr) is dependent on the properties of the porous medium, increasing with finer-grained texture (Johnston and Adamski, 2005). Surface tension between differing fluids also may have an effect on the value of ƒr, although the need for further investigation exists.
|
In contrast to other hysteretic formulations (Kaluarachchi, 1992; Doughty, 2007; Fagerlund, et al., 2008) our model (Figure 2.4) calculates a non-zero Snr for the main wetting-phase drainage process. The value of Snr is updated continuously using the maximum Snmax value of the respective grid-block, regardless of the process. This ensures a continuous function at turning-point saturations. Although this alters the shape of the main drainage curve, the relative permeability curve parameters can be adjusted to move this curve as desired. As the nonwetting phase is displaced from its maximum value (wetting phase imbibitions), the nonwetting phase relative permeability smoothly approaches the residual saturation. In the case of CO2 dissolution, if the value of Sn drops below the value of Snr, Snmax is reset according to Equation 2.6.
|
The hysteretic trapping model was applied to a capillary pressure function based on van Genuchten’s equation (Equation 2.7). The linear trapping model is the same approach used for the calculation of Snr in the relative permeability model.
|
Figure 2.4 shows examples of the new hysteretic characteristic curves. Notice the continuous slope where the imbibition branches meet the main drainage branch—an effect of continuously updating the value of Snr during the wetting-phase drainage process. The curves in Figure 2.3 (Kaluarachchi’s hysteretic curves) were calculated using an Snr value of 0 for the main drainage and a non-zero value for imbibition—the effect of this abrupt calculation change can be seen by the discontinuous slope where the imbibition branches meet the main drainage branch.
Figure 2.4. Examples of the new, smooth hysteretic characteristic curves.
Model Assessment
The first test made with the new hysteretic model was to simulate a NAPL/water experiment published by Johnston and Adamski (2005). The objective of their experiment was to explore the relationship between turning point saturation and residual NAPL saturation using a decane retention cell. Because this experiment involves nonwetting phase trapping in a two-phase system, it is physically similar to a CO2/water system, in which the CO2 is the nonwetting phase.
Figure 2.5 shows the set-up of their experiment. The retention cell held a core sample with a water reservoir connected to one end and a decane reservoir connected to the other end. A hydrophilic ceramic plate separated the core from the water reservoir while a hydrophobic ceramic plate separated the core from the decane reservoir. The plates were placed to ensure that decane would not leak into the water reservoir and that water did not leak into the decane reservoir. The pressure of the decane reservoir was held constant, while the pressure of the water reservoir was varied to apply and release suction to the core. The retention cell was subjected to cycles of subsequently increasing water suction to allow decane to flow into the core under increasing capillary pressures. The suction was released in between cycles to allow water to flow back into the core, pushing the decane out of the larger pores, leaving behind trapped nonwetting residual saturations.
Figure 2.5. Diagram of a retention cell experiment. (Steffy, 1997)
The maximum saturation of decane measured in the core before suction was released was recorded as the turning point saturation, while the minimum saturation of decane measured in the core before suction was re-applied was recorded as the residual NAPL saturation. Over the span of about 90 days, 6 turning point saturation-residual NAPL saturation measurements were taken for the Texas City 3.66 – 4.27 m core sample (Figure 2.6). The measurements of residual NAPL saturations showed an approximate linear relationship with the measurements of the turning-point saturations (Figure 2.7). An estimate of 0.394 was given as the slope of the linear function (ƒr).
Figure 2.6. Residual saturation-turning point data for the Texas City 3.66-4.27 m core sample. Estimates of Sn were made from the measured amount of decane that flowed into the cell and the amount of water that flowed out of the cell. (Johnston and Adamski, 2005)
Figure 2.7. Residual nonwetting phase saturation (Snr) plotted against turning-point nonwetting phase saturation (Sni) from experimental data. (Johnston and Adamski, 2005)
Figure 2.8 shows the 1-dimensional model that was run with TMVOC version of TOUGH2 to simulate this experiment. The Texas City 3.66 – 4.27 m core sample was represented by ten grid-blocks in the center of the model, the ceramic plates were each represented by 2 grid-blocks on both ends of the core cells, and the reservoirs were each represented by 1 grid-block on both ends of the model. The cross-sectional area of the model was kept consistent with that of the core sample from Johnston and Adamski (2005). Because the hydraulic properties of the ceramic plates and the intrinsic permeability of the core sample were unknown, estimated values were used for these parameters in the simulation. The measured main drainage capillary pressure curve for the core sample was fit with Equation 2.7 and used in the model.
Figure 2.8. Model set-up for the experiment conducted by Johnston and Adamski (2005). The cells are colored according to material type.
The grid-blocks representing the hydrophilic plate were given a high nonwetting phase entry pressure in order to repel the decane, and the grid-blocks representing the hydrophobic plate were given a large negative capillary pressure to repel water. The grid-blocks representing the water and decane reservoirs were initially fully saturated in water and decane, respectively, and were held at the desired conditions—the decane reservoir grid-block kept at a constant pressure while the pressure of the water reservoir was controlled to match the experimental conditons. The actual times and values of pressure change were obtained through personal communication with Colin Johnston and were used in the pressure control for the water reservoir grid-block.
The value of ƒr from Johnston and Adamski (2005) was used to calculate Snr in the hysteretic relative permeability function. The fitted van Genuchten parameters from Johnston and Adamski (2005) were used for the capillary pressure function (which was not hysteretic in this simulation). Results from the simulation plotted with the experimental data are shown in Figure 2.9.
Figure 2.9. Experimental data from Johnston and Adamski [2005] plotted with the results from the TMVOC simulation using the new hysteretic nonwetting phase trapping model.
The second numerical test made with the new hysteretic model was to attempt to recreate a 2-dimensional hysteretic simulation of CO2 injection published by Doughty (2007). Doughty’s objective was to compare the behavior of a simulation of CO2 injection and leakage from the storage formation to the surface using her hysteretic characteristic curves and using typical non-hysteretic characteristic curves.
Doughty’s model was a 2-dimensional injection into a formation 100 m thick with leakage through an overlying formation 1000 m thick. Her model extended radially to 40,000 m, and was divided into 61 layers each with 41 grid-blocks. The layers were 20 m thick except for a few thinner layers near the surface. The grid-blocks were 20 m wide extending to 600 m, and then increased steadily to produce an infinite acting model. Both her injection formation and overlying formation were homogeneous, with a vertical permeability of 100 md and horizontal permeability of 200 md. The model was initially fully brine-saturated with a salinity of 100,000 ppm. The model followed a geothermal gradient of 30°C/km with the surface and bottom of the model held constant at 15°C and 48°C, respectively. Her simulation began by injecting 900,000 tons of CO2 for a period of 30 days, and simulation time continued to 1,000 years.
To ensure an accurate comparison of the behavior of our hysteretic formulations, Doughty provided the TOUGH2 input files for her simulation. It was noticed that her published simulations were carried out using ECO2, an older version of the TOUGH2 module ECO2N. The major difference between these modules is that ECO2N includes the density effect that CO2 has in the aqueous phase once dissolved. This turns out to have a dramatic effect on the results at late simulation times. Doughty also provided an updated simulation using ECO2N and shared her new (unpublished) results. Our comparison study will continue using these new simulation results.
We adjusted the parameter values in our relative permeability and capillary pressure curves so that the primary drainage curves are similar to those used in Doughty’s model. We also used a slightly different trapping model (the Land model) to be consistent with Doughty’s approach. The simulation results compare favorably (Figure 2.10), and confirm that the new and simplified hysteretic CO2 trapping function is able to reproduce the essential features of more complicated hysteretic models.
Over the next year, we plan to extend our work on relative permeability and capillary pressure functions to include the unique behavior that we have observed in the recent laboratory exsolution experiment.
Figure 2.10. Results from Doughty (2007), top, and results from our new hysteretic model, bottom.
3. Pore-Scale and Core-Scale Modeling
Our efforts so far in this area have focused on continuum modeling of the laboratory experiments using TOUGH2-ECO2N (see Figure 2.2 above) for an example. As we develop new relative permeability functions, we will continue this type of core-scale modeling.
4. Regional-Scale Variable Density Groundwater Modeling
We have been performing studies to compare the numerical efficiency of the compositional multiphase flow simulator, TOUGH2, with the variable density groundwater flow and transport codes that are based on MODFLOW and MT3DMS. Although the TOUGH2 codes are designed for highly nonlinear multiphase flow problems, they can easily solve single phase variable density flow problems. Our initial comparisons show that it may be computationally feasible to run regional-scale groundwater modeling simulations using the TOUGH2 code. Use of these codes for the groundwater modeling part of this project would greatly simplify the overall modeling effort, because only a single modeling platform is needed to run both the multiphase flow and the groundwater flow simulations. Furthermore, we have parallel computing versions of the TOUGH2 code that can run on multiple processors. This allows us to make use of the Clemson Palmetto supercomputer (which has ~12,000 processors). Using this machine, we can run simulations on dozens to hundreds of processors.
5. Multiphase Flow Simulations of Field-Scale CO2 Injection
Our efforts on this task have been directed towards two main areas: 1) evaluation of dissolved CO2/brine transport from storage formations to drinking water formations through open abandoned wells, and 2) comparison of storage formation sweep efficiency and CO2 distribution during field-scale injection of supercritical versus dissolved CO2. These two areas are described below in separate sections.
A. Wellbore Leakage of Brines Containing Dissolved CO2
Introduction
In order to examine wellbore leakage, it is useful to explore the history of how wells have been completed and subsequently abandoned. Whenever a well is drilled, an open connected pathway is created between the surface and deep formations. Therefore, if wells are not properly sealed upon abandonment, there is potential for fluid migration up the wellbore, potentially contaminating overlying formations.
In the United States, oil and gas wells have been drilled for nearly 150 years. This has led to hundreds of thousands of wells that penetrate the subsurface. Because many of the suitable sites for carbon sequestration are depleted oil and gas formations or are in the vicinity of these formations, it is inevitable that abandoned wells will have to be dealt with when choosing sequestration sites.
Fortunately, many abandoned wells are not likely to pose a major leakage risk to sequestered CO2. Since 1952, well cements have had the appropriate additives to create a proper plug (Ide, et al., 2006). In addition, by the 1950s, most states had put in place sufficient well abandonment regulations. However, prior to this, well abandonment procedures were questionable. Before the 1950s, it is uncertain whether cement plugs were even effective. Early cements lacked the proper additives for the cement to harden at down-hole pressures and temperatures. In fact, such primitive techniques as pouring ice down the well to lower borehole temperature were used to try to achieve proper cementation. Furthermore, it is documented that many of the cement plugs from the Gulf Coast before the 1930s were contaminated with drilling mud (Smith, 1976).
In the early stages of oil and gas drilling, it is likely that many boreholes were abandoned without ever being effectively plugged. When early wells were plugged, they often were filled with whatever material was readily available. Well plugs discovered from the early 1900s have been found to contain materials such as logs, mud, and animal carcasses (Ide, et al., 2006).
Although many of these early wells are poorly sealed, many of them were drilled at relatively shallow depths. Therefore, they may not pose a serious risk to sequestration formations. However, poorly sealed wells do penetrate formations at a depth feasible for carbon sequestration (Gass, et al., 1977). Therefore, the possibility that gaseous CO2 may leak up them is a real concern. While leakage up a completely open wellbore may be unlikely, if CO2 does leak in this way, large volumes of CO2 could be released. In addition to leaking through these unsealed wells, it is possible that injected CO2 could degrade the integrity of properly sealed wellbores. Such leakage pathways include leakage through the cement, through corroded casing, and through the well annulus (Gasda, et al., 2004).
Over time, it is expected that injected CO2 will eventually dissolve into the storage formation brine with time-scales of hundreds to thousands of years (McPherson and Cole, 2000; Ennis-King and Paterson, 2003). In addition, some have proposed injecting CO2 as a dissolved phase (Burton and Bryant, 2007; Leonenko and Keith, 2008; Burton and Bryant, 2009). Fortunately, once CO2 is dissolved into storage formation brines, a buoyant gaseous CO2 phase no longer exists. In addition, upon dissolution, the brine becomes about 1% denser (Enick and Klara, 1990; Bachu and Adams, 2003); therefore, it will have a tendency to sink to the bottom of the formation.
While it seems less likely that dissolved CO2 will leak up abandoned wellbores, there is still a potential danger. If the storage formation is over pressured, or if an overlying aquifer is pumped, the pressure differential could force the dissolved CO2 upward through permeable pathways such as abandoned wells. As it rises, the CO2 has the potential to come out of solution as solubility decreases due to an increase in pressure.
Approach
In order to better understand the factors that control dissolved CO2 leakage such as well permeability, degree of overpressure, and well diameter, computer simulations were performed using the TOUGH2-ECO2N simulator (Pruess, 2007). In addition, a simple analytical model was developed to predict radial single-phase flow from a storage formation into a well, and then into a confined aquifer above. Although the analytical solution is only capable of calculating the volumetric flow rate of the aqueous phase, it is a useful predictive tool for estimating the order of magnitude of the effect of given parameters on the movement of brine containing dissolved CO2.
Analytical Model
The analytical model was developed by coupling steady state radial flow under Darcy’s law with linear Darcy or pipe flow in the wellbore. The well connects a storage formation to an overlying confined drinking water aquifer (DWA). In order to model a wide range of well conductances, from a completely open pipe to a permeability equal to the formations being modeled, two different wellbore flow models were used. For flow in a open or nearly open well, the Darcy-Weisbach equation for open pipe flow was used. For lower well permeabilities similar to the formation, Darcy’s law was used. The model assumes radial flow in both formations. Figure 5.1 shows a conceptual drawing of what the analytical model predicts as well as the variables used in the formula.
Figure 5.1. a) Arrows indicate the direction of flow in the analytical model. b) Drawing showing the variables used in the analytical model.
Solving for the radial volumetric flow rate (Q) into the well from the storage formation results in:
The radial flow out of the well into the DWA is similarly:
If the abandoned well is blocked enough to cause a significant flow resistance, Darcian flow up the well can be assumed:
These equations can be rearranged in terms of the head loss:
At steady state, Qs=Qa=Qw , and the sum of the three head drops (storage formation, well bore, DWA) is equal to the overall head difference between the storage formation and the DWA. Solving for the overall head drop gives:
Equation 5.8 gives the water flow rate up the well as a function of the overall head difference between the two aquifers (assuming constant water density). The term multiplied by the head drop is the overall conductance of the system.
For instances where well permeabilities are high and the flow can no longer be assumed to be Darcian, the wellbore flow term (Equation 5.6) is replaced with the Darcy-Weisbach equation:
Where f is the Darcy-Weisbach friction factor (typically between 0.035-0.1) and g is acceleration due to gravity. Solving for the head drop as in Equation 5.7 results in the following quadratic equation:
Both forms of the model are valid only for single phase aqueous flow of uniform density because hydraulic head is no longer applicable once the fluid density changes due to dissolution of CO2, the presence of brine, or variable temperature. Therefore, the model is not able to exactly predict the flow rate of CO2 laden brine into the DWA, especially if the CO2 comes out of solution. However, flow of CO2 dissolved in brine behaves as a single aqueous phase. Therefore, it is useful in estimating the sensitivity of different parameters on the overall leakage of CO2 up a leaking wellbore.
Numerical Model
The analytical model was helpful for designing a TOUGH2-ECO2N simulation to model dissolved CO2 movement up a leaking well to a drinking water aquifer. The sophistication of the ECO2N simulator allows for modeling the development of dissolved as well as a gaseous plumes in the DWA due to over pressurization of the storage formation and leakage through the well. A radially symmetric numerical grid design was selected for the TOUGH2-ECO2N simulations. Figure 5.2 shows the design of the numerical model.
Figure 5.2. a) Drawing showing the design of the numerical model. b) 2-d representation of the materials used in the simulation.
The model’s outer radius was set to 5,000 meters so that any developing leakage plumes could be contained inside the model and so that end effects could be minimalized. The model was given an overall vertical length of 1,200 meters. Both the storage formation and the drinking water aquifer are 100 meters thick. Between the two formations is a 1000 meter thick impermeable layer. The top of the model is set at a depth of 200 meters below the ground surface.
Hydrostatic pressure was used to generate the pressure gradient such that the pressure in the top formation is ≈2.0x106 Pa and the bottom formation is ≈1.4x107 Pa. In addition, a geothermal gradient of 30°C/km with a surface temperature of 15°C was used. This gradient is typical of the western United States and has been used in other numerical models (Pruess, 2008).
The impermeable layer represents a case where only flow between the storage formation and the DWA are considered. The layer does not transmit fluids, but it does allow for thermal conduction and dissipation of heat from the warmer storage formation water as it rises through the wellbore.
The wellbore consists of the center radial elements, which are given a different permeability from the surrounding elements.
The model was discretized by using 87 grid-blocks in the radial dimension. The thickness of each successive ring increases starting with the well radius (4"; 0.1 m) as the innermost grid-block and extending to a thickness of ≈95 m in the outermost ring. Beyond this outermost ring, a fixed state condition (constant pressure, temperature, concentration, and saturation) was maintained. This allows for flow into the storage formation as well as flow out of the DWA. In the z dimension, a total of 70 grid-blocks were used. The storage formation uses ten 10-meter thick model layers, while the DWA is more refined using twenty 5 meter thick model layers. The impermeable layer was given forty 25-meter thick model layers. This discretization results in a total of 6090 grid-blocks for the entire model.
Table 5.1 shows the material properties for the storage formation, the DWA, and the impermeable layer. All three layers have a rock density of 2600 kg/m3, wet heat conductivity of 2.51 W/m·°C, and a specific heat of 920 J/kg·°C.
Table 5.1. Material properties of the ECO2N model
|
The material properties do not represent any exact geologic location but instead represent typical properties for potential carbon sequestration sites. The relative permeability and capillary pressure parameters come from Doughty, 2007.
For the wellbore, it was desired to minimalize relative permeability and capillary pressure effects. Therefore, the relative permeability parameters for the well were changed so that Slr and Sgr were both 0.02. Capillary pressure was set to zero. In addition, the well porosity was set very high to 0.98. The well permeability varies between simulations.
Results
Analytical Model
In order to test the reliability of the analytical model, the numerical model was used to run simple flow tests with pure water. The volumetric flow rates observed in the numerical model then were compared to the analytical solution for accuracy.
The first comparison was the sensitivity of the flow rate to changes in well permeability. In the numerical model, the bottom formation was given an overpressure of 20 bars to induce flow up the well. Once steady state was obtained, the flow rates given by the model were compared to the calculated flow rate given by the analytical solution.
Table 5.2 shows good agreement between what is predicted by the analytical solution and the flow rate observed in the numerical model. The differences in the two models at lower well permeabilities are probably due to the fact the TOUGH2-ECO2N calculates the water density and viscosity as a function of pressure, while the analytical model assumes these parameters to be constant. The result from Equation 5.11, which uses the pipe flow equation also is shown. The numerical model uses Darcy’s law for all flows. For the case presented in Table 5.2, the overall conductance of the system becomes dominated by the formation conductances once the well permeability exceeds about 10-6 m2.
For another comparison, a numerical simulation was performed at three different storage formation overpressures. A well permeability of 1x10-8 m2 was used in both models. Table 5.3 shows good agreement between the numerical and analytical models at steady state.
Effects of Wellbore Permeability on Dissolved CO2 Leakage
The conductance term in Equation 5.8 is useful for showing the effect of well permeability on the overall flow rate of water through the model. Conductance, which is the inverse of resistance, describes the ease with which the fluid can flow. Thus, the higher the conductance, the higher the flow for a given head difference. Figure 5.3 plots the individual conductances of formations with a permeability of 1x10-12 m2 and a wellbore of various permeabilities against the overall conductance of the entire model. This shows the effect of the wellbore permeability on the overall conductance. At well permeabilities below 1x10-7 m2 the overall conductance is dominated by the wellbore’s conductance. This indicates that the well is the most resistant portion of the model and therefore limits the flow rate. At higher well permeabilities, it can be seen that the overall flow rate is no longer dependent on the well permeability. Instead, the formations provide the most resistance to flow.
|
Figure 5.3. Model conductances as a function of the well permeability.
|
If the numerical model flow rate from Table 5.2 above is plotted against well permeability, it can be seen that the flow rate is controlled in the same way that the conductance term suggests (Figure 5.4), even though this model switches over to the pipe flow equation when the well permeability exceeds 10-8 m2. Over most of the range where the well permeability controls the overall flow rate, the analytical model suggests a 1:1 relationship between flow rate and well permeability. As can be seen in Table 5.2, when the permeability is increased by two orders of magnitude, the flow rate also increases by roughly the same amount. As the well becomes even more permeable, the formation begins to limit the flow rate. As a result, the increase in flow rate is minimal between higher well permeabilities because the formations are limiting the overall flow rate. This indicates that in a leakage scenario, the maximum possible flow rate of the leaking fluid is limited by the permeability of the formations it is flowing out of or into.
Figure 5.4. Flow rate as a function of well permeability.
The analytical model can predict only the behavior of constant density water moving up the wellbore, and it cannot account for CO2 exsolution in the wellbore or formation. Therefore, in order to see more realistically how dissolved CO2 will behave as it moves up the wellbore, the numerical model was used. In the ECO2N model, the storage formation was given a salinity of 2% by mass, which is representative of a low salinity brine. The CO2 concentration in the storage formation was 4.4% by mass. As can be seen in Figure 5.5, this solubility is slightly lower than the maximum solubility at the initial P/T conditions of the storage formation (≈1.4x107 Pa; ≈55°C).
|
Figure 5.5. Dependence of aqueous solubility of CO2 on pressure and temperature.
Upward flow of dissolved brine was induced in all of the models by adding an additional 20 Bar (2x106 Pa) of pressure on the storage formation. This overpressure is consistent with predicted overpressures due to CO2 injection, and it is approximately equivalent to a 200 m head increase. Zhou, et al. (2010) showed overpressures as high as 35 Bar after 50 years of CO2 injection in their models. A range of different well permeabilities was simulated starting with a permeability equal to the DWA of 1x10-12 m2, increasing to 1x10-4 m2.
After 50 years of overpressure, the CO2 moves up the wellbore and accumulates in the DWA. In addition to the dissolved phase leaking into the DWA, the CO2 also forms a separate gas phase near the top of the formation. As can be seen in Figure 5.5, CO2 solubility decreases as the dissolved CO2 rises vertically through the well and pressure drops. Therefore, the mass fraction of CO2 above the new solubility conditions will form a separate gas phase. In the model, the gas phase first appears in the lower portion of the wellbore and continues to increase in volume as pressure decreases up the wellbore. Figure 5.6 shows the dissolved CO2 plume as well as the separate gas phase that forms in the DWA due to leakage up the wellbore at various wellbore permeabilities.
Figure 5.6. a) Dissolved CO2 plume in the drinking water aquifer at various well permeabilities at 50 years. b) Gas plume in the drinking water aquifer at 50 years. Vertical exaggeration is 1.5.
For all well permeabilities, the dissolved plume has two components to its shape. Horizontal accumulation at the top of the DWA is due to the gaseous CO2 dissolving into the resident water. The remainder of the plume represents the dissolved CO2 that has entered the formation through the wellbore. This plume is larger at the bottom of the formation because the dissolved brine is denser than the resident water of the DWA. The gaseous CO2 plume is concentrated largely along the top of the DWA due to its buoyancy and the fact that it comes out of the well at the top of the formation.
Figure 5.6 again illustrates that the amount of CO2 that escapes to the DWA is highly dependent on the well permeability up to a certain point. Table 5.4 and Figure 5.7 show the mass accumulation of CO2 in the DWA at 50 years. The mass accumulations plotted in figure 5.7 show the same dependence on permeability as Figure 5.4. At lower well permeabilities, the mass of CO2 that escapes to the DWA is dependent on the permeability of the well. As well permeability increases, the storage formation begins to limit the mass of CO2 that moves up the wellbore. In addition, as the analytical solution predicts, as the well permeability drops by an order of magnitude, so does the total mass accumulation of CO2.
Table 5.4. Well permeability (m2)
|
1.00x10-4 |
1.00x10-6 |
1.00x10-8 |
1.00x10-10 |
1.00x10-12 |
Gaseous CO2 (kg) |
1.14x108 |
9.29x107 |
2.30x106 |
1.21x103 |
0 |
Dissolved CO2 (kg) |
4.78x108 |
3.76x108 |
1.78x107 |
2.09x105 |
1.59x103 |
Total CO2 (kg) |
5.92x108 |
4.69x108 |
2.01x107 |
2.10x105 |
1.59x103 |
Figure 5.7. Mass accumulation of dissolved and gaseous CO2 in the drinking water aquifer plotted as a function of well permeabilities at 50 years.
Effects of Wellbore Diameter on Dissolved CO2 Leakage
Many different types of wells have been drilled into potential sequestration sites and have the potential to be improperly sealed. In addition to oil and gas production wells, there are exploration wells, test borings, and many other well types (Gass, et al., 1977). Therefore, understanding CO2 leakage effects due to wellbore diameter is important. The ECO2N model again was used with the same properties of 2% salinity and 4.4% dissolved CO2. However, the model discretization was changed such that three different wellbore diameters could be modeled. The model was given a storage formation overpressure of 20 bar and a well permeability of 1x10-8 m2. Figure 5.8 shows the model and the dissolved phase plume in the DWA after 50 years of overpressure with different wellbore diameters. As the wellbore diameter increases, the amount of CO2 that leaks into the DWA increases as well. Table 5.5 compares the increases in total mass accumulation of CO2 after 50 years due to different well diameters. The table shows that the increase is closely related to the increase in the cross-sectional area of the well due to increasing diameter. Thus, doubling the wellbore cross-sectional area would in turn double the mass of CO2 leaked into the DWA.
Figure 5.8. Mass accumulation of dissolved CO2 in the drinking aquifer at various well diameters.
|
Gaseous CO2 (kg) |
Dissolved CO2 (kg) |
Total CO2 (kg) |
Increase in Total Mass |
Well Area (m2) |
Increase in Area |
8" (0.2 m) Well |
1.76x106 |
1.73x107 |
1.91x107 |
|
0.03 |
|
18" (0.457 m) Well |
1.56x107 |
9.42x107 |
1.10x108 |
5.8x |
0.16 |
5.2x |
24" (0.61 m) Well |
3.00x107 |
1.50x108 |
1.80x108 |
9.5x |
0.29 |
9.3x |
Similar to Figure 5.6, the horizontal accumulation of dissolved CO2 along the top of the model in Figure 5.8 occurs due to the presence of gaseous CO2. Therefore, looking at Figure 5.8, it can be seen that the gas phase extends radially much further than the dissolved phase as the well diameter increases. Table 5.5 shows mass accumulations in the DWA after 50 years. It is apparent that the gas phase becomes a larger fraction of the total mass of leaked CO2 as well diameter increases. For the 8" well, the gas phase is 9.3% of the total mass but increases to 14.2% and 16.7% for the 18" and 24" wells.
Effects of Overpressure on Dissolved CO2 Leakage
In any formation where CO2 is injected, there will be some degree of overpressurization. Storage formations likely will be overpressurized the most during injection with the pressure dissipating over time. Simulations by Zhou, et al. (2010) showed overpressures as high as 35 bar with overpressures above 10 bar persisting after injection has ceased. Therefore, to evaluate the effect of differing overpressures on dissolved CO2 leakage through a wellbore, simulations were carried out at three different overpressures. The well was kept at a permeability of 1x10-8 m2 and a diameter of 8" (0.2 m). Figure 5.9 shows the dissolved plume in the DWA at 50 years.
Figure 5.9. Mass accumulation of dissolved CO2 in the drinking water aquifer at various storage formation overpressures.
Due to increasing overpressure, larger amounts of dissolved CO2 are forced into the wellbore and eventually into the DWA. Table 5.6 shows the mass accumulation of CO2 in the DWA due to the three overpressures modeled. As overpressure increases, the gaseous fraction of the total leaked CO2 increases as well constituting 2.3%, 9.3%, and 10% for 10, 20, and 30 bar of overpressure, respectively.
Table 5.6 Overpressure
|
10 bar |
20 bar |
30 bar |
Gaseous CO2 (kg) |
1.28x105 |
1.76x106 |
3.42x106 |
Dissolved CO2 (kg) |
5.45x106 |
1.73x107 |
3.08x107 |
Total CO2 (kg) |
5.58x106 |
1.91x107 |
3.43x107 |
Comparison of Sweep Efficiency Between Supercritical and Dissolved CO2 Injection
The distribution of CO2 in the storage formation following injection is highly relevant to our overall study of risk to drinking water aquifers. Understanding where the injected CO2 is likely to migrate in the storage formation helps us better understand possible failure modes for the injected fluids. One way to characterize the distribution of CO2 in a storage formation is through the sweep efficiency, which is defined below.
There are several proposed methods for geologic CO2 injection. It may be injected, for example, as a supercritical phase (Ennis-King and Paterson, 2002), as CO2 saturated water (Eke, 2009), as brine alternating CO2 cycles (Eke, et al., 2009), brine and CO2 co-injection followed by a chase brine (Qi, 2007; Qi, 2008), or as a CO2 saturated brine (Lake, 1989; Burton and Bryant, 2007; Eke, et al., 2009; Fang, et al., 2010). Of these methods, supercritical CO2 injection and CO2 saturated brine injection have been chosen for sweep efficiency analyses in the current work.
The volumetric sweep efficiency is a measure of the ratio of stored CO2 in an aquifer to available aquifer volume (Wang, et al., 1998). The sweep efficiency is dependent on the injection pattern, reservoir thickness, permeability, areal and vertical heterogeneity, mobility ratio, density differences between fluids, and flow rate. In the following cases the available aquifer volume is calculated disregarding pore compressibility and assuming maximum specific yield equal to that of the effective porosity. This calculation is given in Equation 5.12. The ratio representing volumetric sweep efficiency is given by Equation 5.13.
Here, h is the storage formation thickness, r is the radial extent of CO2 migration from the injection well, and Vinj is the volume of injected fluid (either supercritical CO2, or brine containing dissolved CO2). In addition to the sweep efficiency, we would like to compare the total storage volume required for supercritical versus dissolved CO2 injection.
The prevailing technique for CO2 injection is in the supercritical phase. This method optimizes the mass of CO2 injected by volume. In this phase there is a density difference of about 300 kg/m3 between the CO2 and resident brine. This density difference drives a buoyant migration of CO2 upward toward the sealing layer (Bachu, 2002; Eke, et al., 2009; Benson, et al., 2005). Ideally, the seal would have low permeability and high capillary entry pressure. These properties then would act as a barrier and cause the CO2 to migrate laterally.
In addition to density differences, the viscosity of supercritical carbon dioxide is less than the viscosity of the aqueous phase by a factor of 15. Therefore, supercritical CO2 injection method is classified as immiscible displacement of an aqueous phase by a less dense and less viscous phase (Garcia, 2003; Benson, et al., 2005). Because of the lower density and viscosity of CO2 compared to water, the injection of CO2 will have a tendency to produce viscous and gravitational instabilities, leading to viscous fingering and gravity override (Garcia and Pruess, 2003).
Viscous fingering is caused when a less viscous fluid displaces a more viscous fluid, producing an unstable saturation front taking on the shape of waves or rounded lobes of saturation (Homsy, 1987; Hagoort, 1988; Garcia, 2009). In the case of carbon sequestration, this fingering may allow only a small fraction of the pore volume of a brine-saturated formation to be available for sequestration, while providing a large interfacial region between the carbon dioxide and the brine (Ferer, 2002).
Gravity override refers to the supercritical CO2 overrunning the resident brine because of the lower density of the CO2 (Hagoort, 1988). Because the viscosity ratio for supercritical CO2 and brine is in excess of unity (Garcia and Pruess, 2003), the resultant gravity override is unstable. The combined effect is the formation of a single viscous finger commonly referred to as a gravity tongue, traveling along the top of the aquifer and bypassing the brine (Hagoort, 1988).
These hydrodynamic instabilities may lead to bypassing and poor sweep efficiency (Garcia, 2003). Previous simulation studies have implied that only about 2% of the pore space will contain CO2 if it is injected alone (Fang, et al., 2010). To enhance storage performance, either the probability of leakage must be decreased or the rate at which CO2 is immobilized within the aquifer must be increased (Leonenko, 2007). As suggested by Burton and Bryant (2007), Leonenko and Keith (2007), Eke, et al. (2009), and Fang, et al. (2010) injecting CO2-saturated brine may be a solution to the risks associated with supercritical CO2 injection and a way to enhance storage performance of saline aquifers.
Ex situ dissolution is achieved through surface mixing of CO2 and brine within a pipeline operating at the storage formation’s pressure (Leonenko, 2007). This brine is retrieved from the formation via production wells and re-injected once saturated with CO2. For ex situ dissolution, the volumetric flux of brine must be much larger than that of the CO2. For example, for pressure and temperature conditions that lead to a CO2 solubility of 5%, the brine flow rate must be 20 times that of the CO2 (Leonenko, 2007). The brine production wells used for mixing enable management of reservoir pressure, which can increase storage capacity (Leonenko, 2007). In addition, the production wells can regulate the direction of plume flow and eliminate the large-scale displacement of native brine (Burton and Bryant, 2007).
Dissolving the CO2 prior to injecton changes the nature of the CO2 movement in the subsurface. The density of the solution is greater than separate phase CO2, and it is also denser than the resident brine that does not contain dissolved CO2. This produces a downward buoyancy drive (Burton and Bryant, 2007). This is an advantage over supercritical injection because it removes the need for a perfect seal and can be injected safely at shallower depths than separate phase supercritical CO2 (Fang, et al., 2010). Only single phase flow occurs, eliminating complications associated with saturation fronts, mobility contrasts, fingering and reduced phase permeabilities (Burton and Bryant, 2007). This causes the CO2 to travel at the same rate as the injected brine and it leads to a more uniform sweep efficiency of the reservoir (Eke, et al., 2009; Fang, et al., 2010). Once dissolved, the CO2 may remain underground for millions of years before discharging at the surface (Eke, et al., 2009).
Methods
Using Lawrence Berkeley National Lab’s TOUGH2-ECO2N, local and regional scale multiphase flow models are being created to simulate CO2 injection into deep saline aquifers using both supercritical and dissolved injection. The first of these flow models compares CO2 sweep efficiency following injection into a homogeneous formation. The second set of models compares CO2 sweep efficiency following injection into a stratified heterogeneous formation. The third set of models incorporates heterogeneous, random, spatially correlated permeability fields generated from Lawrence Berkeley National Lab’s iTOUGH2-GSLIB. Model parameters and hydrogeologic characteristics are held constant between comparisons. The same mass of CO2 is injected in each case; therefore, a much greater total mass (brine plus CO2) is injected for the dissolved CO2 case. To accommodate for the resulting pressure increase, injection schemes are being designed so as not to exceed the fracture pressure of the formation. Injection is modeled for 20 years and is followed by a 980 year monitoring period.
Injection into a homogenous storage formation
In initial models, a 200 meter thick homogenous storage formation is represented using a radial grid extending out from the well 60 kilometers. Pressure and temperature gradients are established between the bottom of the storage formation at 20 MPa and 50ºC and the top at 18 MPa and 46.1ºC. All models are run isothermally. The formation is uniformly given a salinity of 50,000 mg/L. The van Genuchten-Mualem Model was chosen to represent relative permeabilities and the van Genuchten function was used for capillary pressure curves. Relevant hydrogeologic properties of the storage formation and saturation curve parameters are given in Table 5.7. Parameter values were determined from Birkholzer et. al. (2008) and Barnes et. al. (2009) which represent studies of similar formations. These values are characteristic of the deep saline aquifers under consideration for carbon capture and storage activities.
Table 5.7. Hydrogeologic properties for the storage formation
Approximately 200,000 tons of CO2 per year are injected along a 200 meter well screen for 20 years. The rate of CO2 injection was chosen as the maximum possible amount injected from any one well in the dissolved phase without overcoming the fracture pressure of the storage formation, set at 3.17 MPa above hydrostatic pressure. In the case of CO2 dissolved in brine, the brine is assumed to have the same composition as the resident fluid. In this case, a volume of brine 20 times that of the CO2 is needed for full dissolution.
Figures 5.10a and 5.10b show the sweep of CO2 at the end of a 20-year injection period. Considering a CO2 mass fraction of 0.005 the cutoff, the farthest the supercritical CO2 travels from the well is 775 meters. This results in a volumetric sweep efficiency of 0.094. The dissolved CO2 travels 1,280 meters at farthest and achieves a volumetric sweep efficiency of 0.632.
Figures 5.10c and Figure 5.10d show the sweep of CO2 80 years after the end of the injection period. Buoyancy and viscous forces drive the supercritical CO2 to the top of the formation and 1,275 meters from the well. This behavior is an example of a gravity tongue, which reduces volumetric sweep efficiency to 0.035. This is a decrease of 63% from the end of injection. The dissolved CO2 travels 1,290 meters from the well and shows only a small buoyancy effect due to its slightly higher density compared to the resident brine. The sweep efficiency is reduced to 0.622, a decrease of only 1.5%.
Overall, the dissolved CO2 sweep efficiency is about 18 times higher than the supercritical CO2 sweep efficiency, but the radial extent of CO2 in the two schemes is similar, due to the dilute nature of the dissolved CO2.
Notably, the dissolved CO2 sinks to the bottom of the aquifer much more slowly than the supercritical CO2 rises to the top. In the first 80 years after end of injection, supercritical CO2 increases its areal footprint on the cap rock by 70%. Dissolved CO2 increases its areal footprint on the bottom of the aquifer by only 1.6%.
Injection into a stratified formation
In order to more realistically simulate regional pressure effects, a system of stacked aquifers and aquitards 2000 meters thick was modeled radially 60 kilometers. Each aquifer is 200 meters thick and each aquitard is 100 meters thick. Within the storage formation (the bottom 200 m zone), three 2-meter thick confining layers were added to increase vertical heterogeneity. It is anticipated that this layering in the storage formation will improve the sweep efficiency during supercritical CO2 injection. The base of the storage formation and below is considered to be impermeable. A cross section of the area surrounding the injection well is shown in Figure 5.11a. Pressure and temperature gradients are established between the bottom of formation at 20 MPa and 50ºC and the top at 3.81 Pa and 12.18ºC. Models are run isothermally. The formation is uniformly given a salinity of 50,000 mg/L. Aquifer properties were held the same as the case of injection into a homogeneous storage formation. Relevant hydrogeologic properties of the aquitards are given in Table 5.8. Again, parameter values were determined from Birkholzer, et al., and Barnes, et al., and are characteristic of the hydrogeology around deep saline aquifers.
Figure 5.11a. Cross section of the area surrounding the injection well.
Table 5.8. Hydrogeologic properties for aquitards
Again, 200,000 tons of CO2 per year are injected along a 200 meter well screen for 20 years. Figures 5.11b and 5.11c show the sweep of CO2 at the end of a 20-year injection period. Measurements are based on a dissolved CO2 mass fraction of 0.005 as a cutoff. The farthest the supercritical CO2 travels from the well is 680 meters. This results in a volumetric sweep efficiency of 0.124. The dissolved CO2 travels 1,290 meters at farthest and achieves a volumetric sweep efficiency of 0.625.
A pressure wave cutoff is set at 0.01 MPa greater than original hydrostatic pressure. In the case of supercritical CO2, a lateral pressure wave reaches 39.25 km from the source and a vertical pressure wave is felt 350 meters above the top of the storage aquifer. In the case of dissolved CO2 injection, a lateral pressure wave is felt 57.5 km from the source and a vertical pressure wave is felt 735 meters above the storage aquifer.
Figures 5.11d and 5.11e show the sweep of CO2 at 80 years after the end of the injection period. The supercritical CO2 travels along the borders of confining layers a maximum of 960 meters from the well. This results in a volumetric sweep efficiency of 0.062, a decrease of 50% from end of injection. The confining layers within the storage aquifer act as barriers to the sinking dissolved CO2 and slow the creep along the bottom of the formation. In 80 years, the dissolved CO2 reaches 1,295 meters from the well. The volumetric sweep efficiency is reduced to 0.62, a decrease of 0.8%. At this point, the pressure waves in all directions of both cases have returned to negligible amounts. Again, the dissolved case has a much higher sweep efficiency (a factor of 10), but the overall volume of the storage area (defined by the outer radius) is 1.8 times larger than for the supercritical injection case.
The homogeneous and stratified injection cases show the impact vertical heterogeneity has on the transport of injected CO2. Particularly in the case of supercritical CO2, vertical heterogeneity increased the storage efficiency by 32% initially. While the dissolved CO2 had slightly decreased storage efficiency in the stratified case, it also had a decreased rate of buoyancy driven flow. This decrease in buoyancy driven flow acts to immobilize the CO2 and enhances the overall storage performance.
Injection into a storage formation with heterogeneous, random, spatially correlated permeability fields
In reality, aquifers are highly heterogeneous in all directions. To increase the hydrogeologic complexity, Lawrence Berkeley National Lab’s iTOUGH2-GSLIB is being used to generate heterogeneous, random, spatially correlated permeability fields for these storage sites. Highly heterogeneous permeability fields promote viscous fingering and may lead to poorer sweep efficiency.
A 200 meter thick heterogeneous storage formation is represented using a three dimensional Voronoi grid covering an area of 3 km by 3 km (shown in Figure 5.12a). Pressure and temperature gradients are established between the bottom of the storage formation at 20 MPa and 50ºC and the top at 18 MPa and 46.1ºC. All models are run isothermally. The formation is uniformly given a salinity of 50,000 mg/L.
|
Figure 5.12a. 3km by 3km by 200m storage formation represented using a Voronoi grid. This image is subject to a vertical exaggeration factor of 2.
Using parameters from Barnes, et al. (2009), the intrinsic permeability is assumed to be log-normally distributed with a variance of 2.06 about an average permeability of 1.14e-13 m2. Sequential Gaussian Simulation is used to create a spherical variogram model with a vertical range of 1 meter for log permeability. A very long horizontal to vertical ratio of 10,000:1 is assumed for the variogram length. Figure 5.12b shows logarithmic isosurfaces of the permeability modifiers applied to the storage formation. All other hydrogeologic properties and saturation curve parameters match those from the injection into a homogeneous aquifer.
|
Figure 5.12b. Permeability modifiers applied to the storage formation. This image is subject to a vertical exaggeration factor of 2.
Figures 5.12c and 5.12d show the three-dimensional sweep of CO2 after 20 years of supercritical injection. The farthest the supercritical CO2 travels from the well is 1,220 meters. This results in a volumetric sweep efficiency of 0.038. The dissolved CO2 travels 1,450 meters at farthest and achieves a volumetric sweep efficiency of 0.492.
Figure 5.12c. Supercritical CO2 saturation after 20 years of injection.
Figure 5.12d. Dissolved CO2 concentration after 20 years of injection. Images subject to a vertical exaggeration factor of 5.
Figure 5.12e shows supercritical CO2 40 years after the end of injection. At this point, the CO2 has reached the model boundaries, increasing its areal extent by 51%. The sweep efficiency 40 years after injection is about 0.03. Figure 5.12f shows the dissolved CO2 100 years after the start of injection. It has not traveled a significant distance since the end of injection and maintains its volumetric sweep efficiency of 0.492. Future simulations will require a larger model domain so that the supercritical CO2 migration can be fully assessed.
Figure 5.12e. Supercritical CO2 40 years after start of injection.
Figure 5.12f. Dissolved CO2 concentration 100 years after start of injection. Images are subject to a vertical exaggeration factor of 5.
Within decades after injection, the supercritical CO2 overcomes the areal footprint of the dissolved CO2, but it is mostly concentrated in a small layer at the top of the injection formation, reflecting its poor sweep efficiency. In contrast, the dissolved CO2 is stagnant, and is fairly uniformly distributed throughout the storage formation. It appears the addition of heterogeneous, random, spatially correlated permeability fields has enhanced the effects of viscous fingering and has lead to poorer sweep efficiency for the supercritical case.
6. Remediation Designs
We have not yet begun work on this task.
Future Activities:
To increase the accuracy of these sweep efficiency models, future simulations will include highly resolved grids and will be run using a parallel version of TOUGH2-ECO2N. The case of injection into a storage formation with heterogeneous, random, spatially correlated permeability fields will be run on a stratified regional scale. Hysteretic relative permeability and capillary pressure curves from Task 3 will be implemented. To reduce numerical error, the grid orientation effect of Voronoi grids versus rectangular grids will be compared in the case of supercritical CO2 injection.
Research will be continued by comparing the analytical model to the results of the variable well diameter and overpressure numerical models.
It is plausible that dissolved CO2 could be forced up an abandoned well without the existence of storage formation overpressure. If an overlying aquifer that is connected to a formation containing dissolved CO2 is pumped, CO2 may be drawn up the abandoned well due to the drop in pressure of the overlying formation. The ECO2N model will be used to model this scenario and examine the effects.
In addition, the complexity of the numerical model will be increased by replacing the impermeable layer with a stratigraphic series of aquitards and aquicludes. This will allow for evaluation of leakage into formations other than just the topmost drinking water aquifer. Also, the effect of different capillary entry pressure and relative permeability properties in the DWA on CO2 leakage plumes will be explored.
References:
Journal Articles:
No journal articles submitted with this report: View all 16 publications for this projectSupplemental Keywords:
Relevant Websites:
Progress and Final Reports:
Original AbstractThe perspectives, information and conclusions conveyed in research project abstracts, progress reports, final reports, journal abstracts and journal publications convey the viewpoints of the principal investigator and may not represent the views and policies of ORD and EPA. Conclusions drawn by the principal investigators have not been reviewed by the Agency.