Environ. Sci. Technol. 2009, 43, 743–749
Model for CO2 Leakage Including Multiple Geological Layers and Multiple Leaky Wells JAN M. NORDBOTTEN,* DMITRI KAVETSKI, MICHAEL A. CELIA, AND STEFAN BACHU Department of Mathematics, University of Bergen, Bergen, Norway
Received May 13, 2008. Revised manuscript received October 3, 2008. Accepted October 15, 2008.
Geological storage of carbon dioxide (CO2) is likely to be an integral component of any realistic plan to reduce anthropogenic greenhouse gas emissions. In conjunction with large-scale deployment of carbon storage as a technology, there is an urgent need for tools which provide reliable and quick assessments of aquifer storage performance. Previously, abandoned wells from over a century of oil and gas exploration and production have been identified as critical potential leakage paths. The practical importance of abandoned wells is emphasized by the correlation of heavy CO2 emitters (typically associated with industrialized areas) to oil and gas producing regions in North America. Herein, we describe a novel framework for predicting the leakage from large numbers of abandoned wells, forming leakage paths connecting multiple subsurface permeable formations. The framework is designed to exploit analytical solutions to various components of the problem and, ultimately, leads to a grid-free approximation to CO2 and brine leakage rates, as well as fluid distributions. We apply our model in a comparison to an established numerical solver for the underlying governing equations. Thereafter, we demonstrate the capabilities of the model on typical field data taken from the vicinity of Edmonton, Alberta. This data set consists of over 500 wells and 7 permeable formations. Results show the flexibility and utility of the solution methods, and highlight the role that analytical and semianalytical solutions can play in this important problem.
1. Introduction Anthropogenic emissions have increased atmospheric concentrations of CO2 by about 35% during the past 200 years, to a current value of about 385 ppm, well above the preindustrial level of 280 ppm (12). If the relentless increase is to be reduced or reversed, technological solutions must be implemented on a massive scale. While many options are being considered, one attractive approach is carbon capture and storage, or CCS (12, 13, 21). CCS involves the capture of CO2 before it is emitted into the atmosphere and subsequent storage away from the atmosphere. Within the storage options for CCS, geological storage appears to be the most promising. Geological storage involves injection of the captured CO2 into deep geological formations, where it would reside for centennial to millennial time scales or longer. There is general * Corresponding author phone:
[email protected]. 10.1021/es801135v CCC: $40.75
Published on Web 12/17/2008
+47
55584869;
2009 American Chemical Society
e-mail:
agreement that the scale of the carbon problem will require deep saline aquifers to be major subsurface repositories for the CO2, in addition to depleted oil and gas reservoirs (with or without enhanced oil recovery) and other subsurface options (12). Success in a geological CCS program will depend on the ability to choose formations from which only minimal amounts of CO2 will leak. This requires formations with sufficient permeability to accept the CO2 injection rates, and which are overlaid by one or more low-permeability units referred to as caprock formations. One leakage pathway of concern involves the large number of wells that may penetrate otherwise excellent caprock formations (12). The centurylong legacy of oil and gas exploration, especially in North America, has left many millions of exploration and production wells, most of which perforate otherwise intact caprock formations. In some locations, a typical plume of injected CO2 might encounter tens to hundreds of existing wells (6). If the large-scale injection of CO2 is viewed in this perspective, then we may have a system similar to that shown schematically in Figure 1. Simulation of flow and leakage in this system constitutes a severe computational challenge that cannot be solved easily and efficiently by standard simulation methods (1). Addressing CO2 storage from a risk assessment (3) or system level (24) perspective will further increase the computational effort required. Herein we describe a semianalytical, grid-free approximation to the CO2 migration and leakage problem. Within the context of a strong set of assumptions, the solution captures the large-scale migration of both CO2 and brine while still resolving the impacts of all of the potentially leaky wells in the system. The model expands on earlier work of the authors in several important ways. It builds from a simpler version of this approach for the system of one injection well and one leaky well, over two permeable layers (20) to include arbitrary numbers of wells, and an arbitrary number of layers of permeable formations. We have incorporated a new formulation for interface upconing along a leaky well (18), and we have developed robust solution methods for both the pressure fields and saturation. Several example calculations are presented to demonstrate the flexibility and utility of the model.
2. Fluid Migration and Leakage Model Injection of CO2 into a deep saline aquifer involves many complicated processes, including two-phase flow dynamics, interphase mass transfer (dissolution), geochemical reactions, and possible nonisothermal effects. While these processes are all interesting and potentially important over different length and time scales, we will focus on the most
FIGURE 1. Diagrammatic representation of pathways for CO2 migration and possible leakage through and along degraded wells. Reprinted with permission from the work of Bachu and Celia (1) (modified from the work of Gasda et al. (6)) 2008 AGU. VOL. 43, NO. 3, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
743
critical processes associated with possible leakage of injected CO2 with an emphasis on leakage along existing wells. Since the largest driving force for leakage occurs during injection, when both pressure and strong buoyancy govern the flow, we focus on the injection phase of a storage operation. In addition, since dissolution and geochemical reactions usually occur over longer time scales (12), we will ignore these processes. Thermal effects are only considered weakly, in that the fluid properties are allowed to vary between aquifers. 2.1. Governing Equations and Model System. We consider injection of CO2 into a deep saline aquifer and therefore will model two fluid phases, namely, brine and a CO2-rich phase that for simplicity we will refer to as CO2. We assume the flow field can be described by the multiphase extension of Darcy’s law, vR ) -
KkR (∇pR + FRg ∇z) µR
(1)
In eq 1, vR [L T-1] is the volumetric flux of phase R, pR [M L-1 T-2] is the fluid pressure, K and kR are the intrinsic [L2] and relative [L0] permeabilities, FR and µR are the density [M L-3] and dynamic viscosity [M L-1 T-1], respectively, g is the gravitational constant [L T-2], and z is the vertical coordinate [L], taken positive upward. The equations of mass conservation can be stated for the carbon (R ) c) and brine (R ) b) phases, ∂φFRSR + ∇ · (FRvR) ) QR, ∂t
R ) c, b
(2)
The additional notation introduced represents the fluid saturation SR [L0], time t [T], porosity φ [L0], and the source (or sink) term QR [MT-1], which in our application will be zero everywhere except at wells. We close the system by requiring that the fluid saturations fill the media, and by ignoring capillary pressure, the effect of which is considered small at large spatial scales. This allows us to use the relations Sc + Sb ) 1,
p c ) pb
(3)
As initial conditions, we take the system to be fully saturated with brine (Sb ) 1) at hydrostatic pressure, pb(x, 0) ) p0 - Fbgz
(4)
for some datum pressure p0, and coordinate vector x ) [x, y, z]T. We model the subsurface as a stack of alternating permeable (aquifer) and low-permeable (aquitard) formations, which is typical of many sedimentary basins. In the current work, we consider all aquitards to be impermeable; however, we are working to relax this assumption in the future. The system is considered to be (partially or fully) penetrated by a number of wells. The aquitards will be impermeable, thus we can omit them from the domain, and take the domain Ω as the Nf aquifers and Nw wells. To enable analytical treatment, we will limit the discussion to homogeneous and horizontal aquifers, and we will consider each well segment (that is, each segment of a well passing through an aquifer or aquitard) as having uniform properties. A discussion of the assumption of homogeneous horizontal aquifers, for this problem can be found in ref 8. Note that each aquifer may have different properties; similarly, each well segment (whether along the same well or different wells) may also have different properties. A CO2 plume arises due to injection into a given permeable layer. The plume spreads into the formation and eventually encounters some number of existing wells. Each well is characterized by a set of effective permeability values, with 744
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 3, 2009
one value assigned to each segment of the well that crosses a caprock formation. If there are Nf - 1 caprock (or aquitard) formations, then each well in the system will have Nf - 1 values of effective permeability assigned to it. Leakage across the caprock formations is restricted to flow along these wells, and is therefore governed by these assigned values of effective permeability. If CO2 leaks along one or more of the wells, the CO2 will migrate upward, where it will encounter one or more of the overlying permeable formations. Depending on the flow conditions, some of the leaking CO2 may flow into these intervening permeable formations, creating new CO2 plumes. These secondary plumes may in turn intersect some number of existing wells, leading to additional leakage, and so on. Assuming this basic model of flow and leakage, the important parameters for the system include the permeability and thickness of each of the permeable layers, the thickness of each caprock formation, the CO2 injection rate, the fluid properties (density and viscosity) of the CO2 and brine, the size of the overall domain and the associated boundary conditions along the boundaries, and the assignment of effective well permeabilities along each segment of each well. Note that in the solution for flow and leakage, the dynamics of both the CO2 and the brine will be modeled. This allows leakage dynamics of both fluids to be simulated and analyzed. 2.2. Approximation of Solution. Most standard multiphase flow simulators can solve eqs 1 and 2 subject to the given boundary conditions (see, e.g, 4, 5, 10). However, for the simulations of interest in the context of CCS, the overall problem requires very large computational domains, while simultaneously requiring local grid refinement around leaky wells. Because we may need to deal with hundreds of wells within the domain of influence of the injection (6) and on the order of 10 aquifers between the injection formation and the surface, high-quality grids conforming to this geometry easily require many millions of grid cells, implying prohibitive computational cost even for large-scale parallel computing. In this section, we will outline a different approach, which uses a much sparser spatial representation, where each node corresponds to a leaky well, rather than to an arbitrary spatial location. This dramatically reduces the computational cost of the simulation, essentially because the number of unknowns is the total number of well segments, e.g., 10 000 for a system with 10 layers and 1000 wells, rather than the millions of unknowns if the flow equations are solved for closely spaced spatial locations (as required by standard multiphase simulators). In deep terrestrial systems, the injected CO2 will be less dense than the resident brine, with CO2 density between 25% and 75% of the brine density (19). The resulting strong buoyant drive leads to gravity segregation, which makes the assumption of a macroscopic sharp interface reasonable. Under the associated assumption of essentially horizontal flow, we can write the vertically integrated form of eq 2 for an aquifer of thickness H as ∂ ∂t
∫
H
0
φFRSR dz + ∇ ·
∫
H
0
FRvR dz )
∫
H
0
QR dz,
R ) c, b (5)
We assume a sharp interface separating CO2 (with residual brine saturation) above the interface and brine below the interface. Define the thickness of the CO2 layer as h(x, t). With the sharp interface assumption, h(x, t) fully describes the saturation distribution in the aquifer. Assuming negligible vertical variation in porosity and densities, eq 5 simplifies to
( ∫
∂ (φFc(1 - Sres)h) + ∇ · Fc ∂t
H
0
) ∫
vc dz )
H
0
Qc dz
(6)
( ∫
∂ (φFb(H - h) + φFbSresh) + ∇ · Fb ∂t
H
0
) ∫
vb dz )
H
0
G(x, xi, h, t) )
Qb dz (7)
Applying the vertical equilibrium assumption, which is consistent with the assumption of essentially horizontal flow, yields the vertically integrated flux from eq 1 as,
∫
H
0
vc dz ) -
∫
H
0
Kkch ∇ (pB + h(pb - pc)g) µc
(8)
K(H - h) ∇ pB µb
(9)
vb dz ) -
∂p ˜i ∂Q
(17)
These quantities represent the sensitivity of the pressure field to the sources and sinks. As such, they are effectively Green’s functions for a linearization of the pressure solution. Such a linearization is appropriate for a single well operating at a constant rate of injection (see eq 20 later and ref 16). Assuming that this linearization is appropriate for the full system, we can construct the solution by superposition, Nw
p(x, t) ) p0 +
∑ Q˜ G(x, x , h, t) + F(x, h, t) i
i
(18)
i)1
Here kc is the relative permeability of CO2 at the residual brine saturation, pB(x, y, t) ) p(x, y, zl ) 0, t), where zl is the local vertical coordinate in the aquifer. The pressures at the top and bottom of the aquifer, pT and pB, respectively, are related through pT ) pB - (H - h)Fbg - hFcg
(10)
Summation over eqs 6 and 7 weighted by FR-1 leads to h
(
)
1 - Sres ∂ Sres ∂ 1 ∂ 1 φF + φF + (H - h) φF + ∇ · Fc ∂t c Fb ∂t b Fb ∂t b Fc H H H 1 ˜ dz (11) Fc vc dz + ∇ · Fb vb dz ) Q 0 0 0 Fb
( ∫
)
( ∫
) ∫
We denote the total volumetric source (or sink) term as ˜) Q
Qc Q b + Fc Fb
(12)
To proceed to an analytically tractable pressure equation, we make two assumptions. First, we approximate the compressibility of CO2 and brine as equal and define the effective compressibility as ceff ≡
(
)
Sres ∂ 1 - Sres ∂ 1 ∂ φFc + φF ) φF Fc ∂p Fb ∂p b Fb ∂p b
(13)
Second, let the spatial variability of density be negligible within a given fluid phase, within a given aquifer. Then, eq 11, simplifies, ceff
(
)
Kkch(Fb - Fc)g 1 ∂p - ∇ · Kλeff ∇ p + ∇h ) ∂t µcH H
∫
H
0
˜ dz Q (14)
Here we have introduced the effective mobility given by (note that the relative permeability of brine kb ) 1) λeff )
kch H - h + µcH µbH
(15)
Finally, within an aquifer, the source terms Qc and Qb appear as point sources located at the wells, Nw
QR(x) )
∑ δ(x - x )Q i
R,i
(16)
The offset F appears as a consequence of considering the pressure at the bottom of the aquifer (recall that p ) pB), which is thus dependent on the vertical fluid distribution at that point. In terms of standard pressure-saturation formulations, pB - F is analogous to the so-called global pressure for the system. We discuss the offset function F further in section 2.3.1. 2.3. Analytical Solutions. To complete our model description, we apply two additional analytical components. The choices of these components are not intrinsic to the modeling concept and could be substituted for other models depending on the accuracy needed (see, e.g., ref 9 for preliminary work on integrating traditional numerical methods). As examples, the injection response described in section 2.3.1 might be generalized to include caprock flow following the ideas of Hunt (11). Conversely, the upconing model given in section 2.3.2 can be simplified by using e.g. Muskat’s approximation (15). And we could replace the assumed Darcy flow in leaky wells by more complex pipeflow models if that were appropriate. Our first component is a similarity solution for two-phase injection into porous media. This solution provides an expression from which estimates of the two primary variables, pressure pB(x, t) and interface location h(x, t), can be made for any individual plume. These estimates allow the functions G and F to be calculated, thereby forming the basis for an overall solution algorithm. The second component is a generalization of traditional upconing models for two-phase flow into a well. This solution provides us with the near-well flow field and saturation distribution needed to more accurately compute flows along any given well. 2.3.1. Fluid Injection into Confined Aquifers. In a previous investigation (16), we considered the case of the constantrate injection into confined aquifers initially filled with brine. The solution used the sharp interface assumption. It was shown that the saturation distribution, given by h(x, t) has a self-similar solution in time, such that the solution could be expressed in terms of a single space-time variable h(x, t) ) h(ri2/t), where ri ) ||x - xi||2 is the distance to the well. In self-similar coordinates, the equations governing the fluid distribution and pressure are (see eqs 11 and 17 of ref 16) dh′ 1 - h′ dh′ 2 d 2Γλhχ ′ ) +1 dχ χ dχ h′(λ - 1) + 1 dχ
[
∆p′ ) -
i)1
where δ(x - xi) is the Dirac delta function centered at the well location x ) xi. Consider the formal solution p of eq 14 at some time t. If the inverse of compressibility is large relative to the perturbation from initial pressure in the system, which means that pressure equilibrates fast relative to saturation changes, we can approximate the pressure equations by superposition and introduce the quantities
1 2Γ
(
∫
ψ
χ
)]
dχ + F(h′) (h′(λ - 1) + 1)χ
(19) (20)
The dimensionless variables and parameters in these equations are defined as χ) Γ)
2πφH(1 - Sres)ri2 Mc(t) ⁄ Fc
2π(Fb - Fcw)gKλbH2 , ˜ Q
(21) λ)
λc λb
VOL. 43, NO. 3, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
(22)
9
745
∆p′ )
p - p0 , (Fb - Fc)gH
h′ )
h H
(23)
where Mc(t) is the cumulative injected mass of CO2. By comparison of eq 20 to eq 17 of the work of Nordbotten and Celia (16), we see that the pressure offset F in eq 19 due to the vertical fluid distribution given by F(h′) ) -
ln[(λ - 1)h′ + 1] λ h′(χ) λ-1 λ-1
[
]
(24)
The compressibility of the system, ceff, enters through the outer limit of integration in the pressure equation, and is analogous to successive steady-state approximations (see refs 18 and 20), Ψ)
4.5πHφλbK(1 - Sres) ˜ ceffQ
(25)
The advantage of using the self-similar approximation given in eq 19 is that it replaces a system of 3D nonlinear partial differential equations (eqs 6 and 7) by a single nonlinear ordinary differential equation. Even if the latter is solved numerically, this strategy produces tremendous computational savings. For the limit of Γ f 0, analytical solutions to eqs 19 and 20 exist, which allows for simple closed-form analytical solutions. 2.3.2. Fluid Upconing near Wells. Understanding the nearwell behavior of multiphase flow has a long history in both groundwater and petroleum engineering. This is also an important issue in the CO2 leakage problem, because nearwell flows may strongly influence the saturation in the wells and, hence, larger-scale flow patterns. Analytical solutions to near-well flow are typically obtained by either (i) assuming vertical equilibrium or (ii) applying perturbation equations around a known solution. Since potentially large vertical flows near leaky wells may cause significant local upconing effects, neither of these traditional approaches gives satisfactory results in the current application. Therefore we use a new solution approach, which partially relaxes the vertical equilibrium assumption by explicitly approximating vertical flow (17). Single-phase flow in the well (Sc ) 1) occurs when the upconing, near the well, is not sufficient to overcome the fluid layering. This will be the case when the flow rates in the well are small, or when the thickness of a CO2 plume at a well is large. Using eqs 30 and 31 of ref 17, we see that this situation occurs when the flow of the dense fluid (brine) in the well is zero, Qd′ ) 0. Conditions for this are given by
(
ηouter2 < ηwell2 +
k′ k′ exp((h0)2Γλ) 4Γλ 4Γλ
)
(26)
Here η ) ri/Kh is a dimensionless radial distance scaled by the horizontal permeability, and the dimensionless (inverse) vertical permeability is given by k′) H2/Kz. The outer radius of influence of the upconing is denoted by ηouter. For singlephase flow, it is well-known (see, e.g., ref 2) that at steady state
( )
ηouter ) cη
˜ Q vˆ H2
Rη
(27)
where vˆ is the background flow rate (neglecting the effect of the well), and the exponent Rη ) 1. For two-phase flow, we obtain good results using Rη ) 2. When inequality 26 is violated, two-phase flow occurs in the leaky well, as brine is pulled up into the well and flows with the CO2. We describe this with the approach of ref 17, referring to the online appendix for the relevant equations. 2.3.3. Well Flow Model. Flow in abandoned wells is in general a complex process, which may include open well 746
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 3, 2009
flow and significant thermal effects (22). These complex effects may be important and may be incorporated into our modeling framework. However, in the applications presented herein, we will not focus on the complexities of the well flow, and ignore nonisothermal effects in the flow. We assume that flow along wells occurs mainly through fractured or degraded porous media along the outside of the well casing, including well cements and drilling fluids, and that the smallscale flow paths can be represented at the larger scale by an effective Darcy permeability (see right-hand part of Figure 1). 2.4. Overall Solution Algorithm. The major components of our overall solution algorithm have been presented in the previous sections, including solutions for the space-time evolution of both the CO2 plume and the pressure field, and the determination of leakage conditions along a well through the use of an upconing model. These form the basis of our solution algorithm, which we outline as follows. Let the state of the system be well defined at time tn (e.g., at the initial condition t0 all CO2 masses are zero and the initial pressure distribution is hydrostatic in all aquifers). To calculate the solution at time tn+1 ) tn + ∆t, we use the known CO2 plume masses from time tn to calculate a new pressure field at time tn+1. This pressure field is defined for all space, using the Green’s functions and the flow rates at the individual wells. This should be contrasted to traditional multiphase simulators, which discretize the pressure field on a much denser computational grid. Mass conservation along segment l of well i is imposed using d l(t) M ) QR,i dt R,i
∑
l+ QR,j (t)
(28)
l (t) j∈SR,i
This equation is written for cumulative masses MR,i corresponding to plumes of the CO2 (R ) c) and brine (R ) b) phases associated with well segment i. It states that the mass added to a plume associated with segment i is equal to the sum of the fluxes into the plume, minus the sum of fluxes out of the plume. Vertical flows along any well segment are described by the one-dimensional multiphase version of Darcy’s law, written in discrete form as l QRl ) -πrw2K w
(
)
kR pBl+ - pTl+ FRg FR µR Hl
(29)
Equation 29 is used to describe the flux through a well perforating caprock (or aquitard) layer l. We denote the well by a subscript w, and the permeable formations immediately above and below this layer are denoted by l+ and l-, respectively. Note that we distinguish between pressures at the top and bottom of permeable formation through the use of subscripts T and B, respectively. As the wells are treated as 1D objects intersecting the aquifers, their internal volume is considered negligible. Therefore, saturation of each phase within the well segment becomes a key parameter in the Darcy equation (eq 29) which governs the value of the relative permeability kR and hence controls the flow rate. Whenever a neighbor plume intersects the well in the formation immediately below caprock section l (that is, formation l- in eq 29), the upconing solution of section 2.3.2 is used to determine flow rates for both fluids in the well segment. Given the flow rates, the fluid saturation can be determined using the relative permeability function (which is assumed to be knownshere we used linear functions, but nonlinear forms can also be handled). For wells that are not intersected by neighboring plumes (so there is no upconing), we use Scl ) (Scl-1)ν where Scl-1 is the saturation in the segment below, and ν ≈ 0.5. This specifica-
tion was determined by comparison with saturation profiles observed in detailed, fully numerical simulations of the system. The overall solution algorithm involves the following steps: (1) solve the linear system of pressure equations, which yields pressures at each of well segments in each permeable layer in the system, (2) use these pressures, in conjunction with the upconing equations, to determine fluxes and saturations in each well segment, (3) update all plumes and the associated cumulative masses, and (4) move on to the next time step. This numerical algorithm is analogous to the implicitpressure-explicit-saturation (IMPES) algorithm sometimes used in petroleum reservoir engineering. Several additional assumptions are needed when carrying out the actual calculations. First, when calculating saturations within the domain, we deal with overlapping plumes by simply taking the maximum thickness, h ) max hi i
(30)
Second, when calculating the pressure solutions, the assumption of radial symmetry inherent in the solutions of section 2.3.1 requires that all parameters in the equation have radial symmetry. However, since CO2 and brine have different viscosities, and there are multiple (locally radial) plumes throughout the domain, the pressure field will fail to have global radial symmetry. We construct an approximate radially symmetric field using the procedures outlined in ref 20 for the two-well problem. Third, note that the solution of ref 16 was derived for a constant injection rate. This is clearly not the case for leaky wells (which act as injection wells in upper layers), where their flow rates vary in time, increasing from zero to a nearconstant value. We account for this in two ways: In eq 21, the denominator of the key dimensionless group χ involves the cumulative mass in the plume, and not simply the (constant) flow rate multiplied by time. This corresponds to using the current rate dMc/dt at time t to determine an effective time t* that would be required, at that flow rate, to reach a total mass of Mc. This effective time t* is then used in the similarity solution. Additionally, when calculating the outer limit for the pressure perturbation, denoted by Ψ in eq 25, the correction to account for the variable total flow rates is given by the modified equation Ψi )
4.5πHφλbK(1 - Sres) |Vtot,i| ˜ i| ceffMc,i/Fc |Q
(31)
Here Vtot,i represents the total fluid volume which has entered (or left) the formation through the given well segment i. The absolute value signs appear because we use the same approach for both positive and negative flows.
3. Representative Results The approximation and solution strategies detailed in the previous sections provide the basis for an efficient numerical simulator. This section illustrates some of the capabilities of our implementation. We present two examples in this section. The first example serves as a validation of the model by comparing it to the results of a traditional numerical model. The second example demonstrates the applicability to large field-scale calculations, which is far beyond the capabilities of traditional methods. 3.1. Eight Well, Five Aquifer Example. Our first example is a representative case from a suite of validation experiments we have conducted. This case has five aquifers, penetrated by nine wells which form a three by three square lattice when seen from above. We consider an isotropic permeability of 10 mD in all formations, with a porosity of 10%. All formations
FIGURE 2. Comparison of leakage into successive formations between the new model (blue) and an established numerical simulator (green). Aquifer numbering is from the bottom up and is indicated by the patterns shown with black color in the legend. are modeled as 20 m thick, with an areal extent of 250 km2. Boundary conditions are constant hydrostatic pressure. The densities of brine and CO2 are set to constant throughout the domain, at typical values of 1000 and 600 kg/m2, respectively, with viscosities of 0.5 and 0.05 mPa · s. The distance between the central well, which we take as the injection well, and the four nearest (abandoned) wells is 1000 m. All abanonded wells have identical properties, with permeabilities of 100 mD and radii of 20 cm. This problem therefore has an 8-fold symmetry. Such a symmetrical setup allows for sufficient grid resolution in the numerical simulator that will be used to validate our new approximation method. We simulated 2500 days of injection at 50 kg/s and compare results obtained by two methods: (i) the semianalytical model described herein and (ii) the industry-standard simulator Eclipse (23). For the numerical Eclipse simulation, we discretized a quarter of the domain (exploiting symmetry), with 37 6320 grid cells. On comparable computational resources, our model reduces the run-time by more than 3 orders of magnitude, with similar savings in memory usage. In practical terms, this replaces days of computing with less than a minute on a standard desktop machine. In Figure 2, we show accumulated mass from all wells in each layer. Note the logarithmic scale on the y-axis. The aquifer numbering is from bottom up, starting with the injection aquifer. For early times, there is a substantial deviation between the models, due to CO2 arriving at the leaky wells significantly earlier in the numerical simulation. We attribute this, at least in part, to numerical diffusion in Eclipse. At late times, there is excellent agreement between the numerical simulation and our model. 3.2. Full-Scale Field Example. To demonstrate some of the capabilities of our new model for CO2 injection and leakage, we have applied the model to a data set associated with a field location close to Edmonton, Alberta, Canada. At this location, which is shown in Figure 3, four large power plants collectively emit about 35 million tons of CO2 per year. As such, this is a potential location for carbon capture and storage. Within a domain that is 30 km × 30 km, 508 existing wells were identified, and the local layered structure of the subsurface was characterized (see http://www.ags. gov.ab.ca/activities/wabamun/Wabamun_base.html). These data form the basis for a model system with seven permeable layers separated by six impermeable layers. These layers were assigned nominal permeability values of 10 milli-Darcy (mD), which is characteristic of the Alberta Basin. The existing wells were assumed to penetrate all of the layers in the model. While the statistical properties of effective permeability values along well segments are currently unknown, we use VOL. 43, NO. 3, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
747
FIGURE 3. Characteristics of the Wabamun Lake area, Alberta, Canada: (a) location and major stationary CO2 sources and (b) distribution of deep wells (oil and gas and injection wells). Reprinted with permission from the work of Bachu and Celia (1) 2008 AGU.
FIGURE 4. Leakage into successive formations from a hypothetical injection at a location in Alberta. Aquifer numbering is from the bottom up. values herein as illustrations of the variability we expect to see. Experimental methodologies for measuring the permeability of abandoned wells and ongoing field campains are discussed elsewhere (7). Thus, each well segment was assigned an effective permeability value sampled randomly from a log-normal probability distribution. Here, we assumed that log10 k ∼ N(µ ) log10(10 mD), σ2 ) 1), which corresponds to a most likely permeability value of 10 mD and a variance of about an order of magnitude. In this example, we show a single realization of parameters from this distribution; a detailed investigation of the model behavior in response to the input parameter will be presented in a forthcoming publication. Injection was modeled over a characteristic time period of 30 years at 500 kg/s. Boundary conditions are constant hydrostatic pressure. Note that while we only model a single injection well in this run, this may also be considered as a proxy for a cluster of injection wells. The accumulation of CO2 in the six permeable layers above the injection layer, formed from a series of secondary plumes, is shown in Figure 4. Similarly to the previous results, we see a clear reduction in the amount of leakage that occurs as additional layers are encountered moving upward in the stratigraphic sequence. This shows the importance of intervening permeable layers above the first impermeable layer and their ability to mitigate leakage rates. The only layer where this pattern is not followed is the upper-most permeable layer, where a no-flow boundary above it leads to a somewhat artificial accumulation of mass. In this sense, the upper layer in our model represents the 748
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 43, NO. 3, 2009
FIGURE 5. Leakage statistics in terms of plume masses for successive formations from a hypothetical injection at a location in Alberta. Aquifer numbering is from the bottom up, and the units on the x-axis refer to the logarithm of mass in kilograms. The left y-axis is for the probability density curves, and the right axis refers to number of well segments with plumes of mass equal to 1 kg or less. total mass that may leave the geological system and should be thought of as a proxy for the shallow groundwater and atmosphere. These results show the importance of including all layers in an analysis of an injection scenario, with the concomitant implications for site characterization and regulatory requirements. We note that the rate of brine flow is also captured in the model and shows a similar trend to CO2 leakage. However, due to lack of buoyancy, in some cases we observe negative brine fluxes, associated with countercurrent flow of brine. A discussion of brine leakage patterns is included in the online Supporting Information for this publication. Additional insights given by our simulator are illustrated in Figure 5, which shows the distribution of secondary plume sizes, sorted by layer. The largest plumes occur in the layer immediately above the injection layer, with a distribution clearly skewed toward higher values. Plumes become progressively smaller as we move upward in the stratigraphic succession, and the distributions become more symmetric. Again, the highest permeable layer is somewhat anomalous in its number of plumes, due to the boundary condition at the top of the domain. The atom of probability at 1 kg (labeled 100 in the figure) represents all smaller plumes, almost all of which are insignificant. This information about plume sizes and their overall spatial distribution can provide practical guidance for activities such as the design of subsurface monitoring and leakage detection strategies.
Supporting Information Available Upconing equations used and the flow of brine predicted for the final example. This material is available free of charge via the Internet at http://pubs.acs.org.
Literature Cited (1) Bachu, S. and Celia., M. A. Assessing the Potential for CO2 Leakage, Particularly through Wells, from CO2 Storage Sites. The Science and Technology of Carbon Sequestration, AGU Monograph; American Geophysical Union: Washington, DC, 2008, in press. (2) Bear, J. Dynamics of Fluids in Porous Media; Elsevier: New York, 1972. (3) Celia, M. A.; Bachu, S.; Nordbotten, J. M.; Kavetski, D.; Gasda., S. Modeling critical leakage pathways in a risk assessment framework: Representation of abandoned wells. In Proceedings of the 4th Annual Conference on Carbon Capture and Sequestration, May 2-5, 2005, Alexandria, VA.
(4) Class, H.; Bielinski, A.; Helmig, R.; Kopp, A.; Ebigbo, A. Numerical simulation of CO2 storage in geological formations. Chem. Eng. Tech. 2006, 78 (4), 445–452. (5) Doughty, C.; Pruess, K. Modeling supercritical CO2 injection in heterogeneous porous media. Vadose Zone 2004, 3 (3), 837– 847. (6) Gasda, S. E.; Bachu, S.; Celia, M. A. Spatial characterization of the location of potentially leaky wells penetrating a deep saline aquifer in a mature sedimentary basin. Environ. Geol. 2004, 46 (67), 707–720. (7) Gasda, S. E.; Nordbotten, J. M.; Celia, M. A. Determining the effective wellbore permeability from a field pressure test: Numerical analysis of detection limits. Environ. Geol. 2008, 54 (6), 1207–1215. (8) Gasda, S. E.; Celia, M. A.; Nordbotten, J. M. Upslope Plume Migration and Implications for Geological CO2 Sequestration in Deep, Saline Aquifers. IES J. A: Civil Struct. Eng. 2008, 1 (1), 2–16. (9) Gasda, S.; Nordbotten, J. M.; Celia., M. A. Vertical Equilibrium with Sub-scale Analytical Methods for Geological CO2 Sequestration. Comput. Geosci., submitted for publication. (10) Juanes, R.; Spiteri, E. J.; Orr, F. M.; Blunt, M. J. Impact of relative permeability hysteresis on geological CO2 storage. Water Resour. Res. 2006, 42 (12), W12418. (11) Hunt, B. Flow to a well in a multiaquifer system. Water Resour. Res. 1985, 21 (11), 1637–1641. (12) Metz, B., Davidson, O., de Coninck, H. C., Loos, M., Meyer, L. A., Eds. IPCC Special Report on Carbon Capture and Storage; Cambridge University Press: New York, 2005. (13) The Future of CoalsOptions for a Carbon-constrained World; Massachusetts Institute of Technology: Cambridge, MA, 2007. (14) U.S. Greenhouse Gas Abatement Mapping Initiative. Reducing U.S. Greenhouse Emissions: How Much and at What Cost?; McKinsey and Company, 2007.
(15) Muskat, M. Physical Principles of Oil Production; McGraw-Hill: New York, 1949. (16) Nordbotten, J. M.; Celia, M. A. Similarity solutions for fluid injection into confined aquifers. J. Fluid Mech. 2006, 561, 307– 327. (17) Nordbotten, J. M.; Celia, M. A. An improved analytical solution for interface upconing around a well. Water Resour. Res. 2006, 42, W08433. (18) Nordbotten, J. M.; Celia, M. A.; Bachu, S. Analytical solutions for leakage rates through abandoned wells. Water Resour. Res. 2004, 40 (4), W04204. (19) Nordbotten, J. M.; Celia, M. A.; Bachu, S. Injection and storage of CO2 in deep saline aquifers: Analytical solution for CO2 plume evolution during injection. Transport Porous Media 2005, 58 (3), 339–360. (20) Nordbotten, J. M.; Celia, M. A.; Bachu, S.; Dahle, H. K. Semianalytical solution for CO2 leakage through an abandoned well. Environ. Sci. Technol. 2005, 39 (2), 602–611. (21) Pacala, S.; Socolow, R. Stabilization Wedges: Solving the Climate Problem for the next 50 Years with Current Technologies. Science 2004, 305, 968. (22) Pruess, K. On CO2 fluid flow and heat transfer behavior in the subsurface, following leakage from a geologic storage reservoir. Environ. Geol. 2007, 54 (8), 1677. (23) Schlumberger Information Systems. Eclipse technical description; Report, Houston, TX, 2004. (24) Stauffer, P. H.; Viswanathan, H. S.; Pawar, R. J.; Guthrie, G. D. A system model for geologic sequestration of carbon dioxide. 2009, 43, 565-570.
ES801135V
VOL. 43, NO. 3, 2009 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
749
