Estimates of CO2 leakage along abandoned wells constrained by new data

International Journal of Greenhouse Gas Control 84 (2019) 164–179

Contents lists available at ScienceDirect

International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc

Estimates of CO2 leakage along abandoned wells constrained by new data Tom J.W. Postma , Karl W. Bandilla, Michael A. Celia ⁎

T

Department of Civil and Environmental Engineering, Princeton University, 59 Olden St., Princeton, NJ 08540, USA

ARTICLE INFO

ABSTRACT

Keywords: Geological storage Leakage Semi-analytical solutions Abandoned wells Risk analysis Field data

The viability of carbon capture and geological storage (CCS) projects depends in part on the risk that injected CO2 or displaced pore fluid will leak out of the target formation into surrounding formations or to the surface. Abandoned oil and gas wells, of which millions exist both throughout the United States and globally, form a potential conduit for this leakage. Recently, specific field measurements have been made to quantify the range of effective permeabilities that can be expected in abandoned wells, enabling us to, for the first time, combine fieldscale numerical simulations of CO2 sequestration in deep saline aquifers with real data on effective permeabilities of leaky wells. Using a previously developed semi-analytical reservoir simulator that can accommodate an arbitrary sequence of alternating aquifers and aquicludes, as well as an arbitrary number of leaky wells, we investigated how the amount of CO2 that leaks out of the target formation depends on the spatial density of nearby abandoned wells and their effective permeability. Furthermore, we assess the influence that variations in pressure and temperature found between geological targets have on this dependency. We find that the observed differences in leakage between geological targets are controlled almost exclusively by differences in density of CO2 at the local subsurface conditions, causing the CO2 plume to contact a different number of wells when injecting at the same constant mass rate. We quantitatively assess the results obtained from our numerical experiments by combining them with the permeability data that have recently become available, typical spatial densities of abandoned wells, and performance requirements put forward in the literature. Our results indicate that leakage of CO2 through abandoned wells is unlikely to be a major limitation in storage security of CCS projects.

1. Introduction According to the IPCC, the “global carbon budget” that can be emitted while retaining a two-thirds probability of keeping global temperature rise from pre-industrial levels below two degrees Celsius is around 1000 gigatons (Gt) of CO2 (IPCC, 2014). With yearly energyrelated emissions reaching a record high of 32.5 Gt in 2017 (IEA, 2017), transformative changes to the global energy system will have to be made on a time scale of years to decades. Carbon capture and storage (CCS), in which large volumes of CO2 are sequestered in porous formations in the deep subsurface, has been identified as a key technology in almost every viable plan for handling the global carbon problem, providing a way to keep atmospheric concentration of CO2 at acceptable levels as the world transitions into a carbon-neutral energy economy (Pacala and Socolow, 2004). In terms of the “stabilization wedges” put forward in the paper by Pacala and Socolow, the amount of CO2 that needs to be permanently sequestered to contribute roughly one tenth of the emissions reduction required is 25 Gt, occupying a volume of around 2 × 1011 m3 at reservoir ⁎

conditions. In their analysis, this volume will have to be stored over a period of 50 years, and it is worth noting that this volume is almost as large as the projected total global oil production in that same period (Celia et al., 2015). The viability of any CCS project depends in part on storage security, i.e. the risk that injected CO2 (or displaced pore fluid) will leak out of the target formation into surrounding formations or to the surface. To successfully isolate the injected CO2 from the atmosphere, four trapping mechanisms have been identified, each with its own associated time scale (Fig. 1). During the time of active injection and until the CO2 has been immobilized (either by dissolution or capillary trapping), the injected CO2 exists as a mobile, CO2-rich phase distinct from and of significantly lower density than brine. This means that the CO2 plume is being retained only by the integrity of the seal above the injection target. It is during this phase of the CCS operation that the risk of unintended migration is at its highest, especially in the time of elevated pore fluid pressure associated with the active injection period (Nordbotten et al., 2009). Various potential pathways for unintended migration across the seal

Corresponding author. E-mail address: [email protected] (T.J.W. Postma).

https://doi.org/10.1016/j.ijggc.2019.03.022 Received 14 October 2018; Received in revised form 14 March 2019; Accepted 22 March 2019 1750-5836/ © 2019 Elsevier Ltd. All rights reserved.

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

plugged and abandoned wells with the injected CO2 plume will be difficult to avoid – particularly considering the scale at which CCS needs to be deployed to achieve climate change mitigation targets, as discussed in the beginning of this section. Various studies have been conducted exploring the multiple leakage pathways through plugged and abandoned wells, as well as the possibility of CO2–brine–cement interactions that might increase or decrease the risk of leakage (for reviews, see e.g. Carey (2013), Choi et al. (2013), Zhang and Bachu (2011)). In general, these reviews state that the risk of geochemical reactions creating defects in successfully plugged and abandoned wells is small, but a flow of CO2-rich brine along pre-existing defects might change the condition of the plug depending on the circumstances. Various simulation studies have been put forward that assess leakage through abandoned wellbores at specific CCS sites, either explicitly modelling well abandonment and completion details (e.g. Pawar et al. (2009)), or by representing the leaky wells as one-dimensional flow conduits with a certain effective permeability, thereby not having to specify details about abandonment, completion or leakage mechanisms (e.g. Celia et al. (2011), Nicot et al. (2013)) In another study, Cihan et al. (2012) estimated far-field pressure buildup and possible migration of resident brine when injecting CO2 in a multilayer system, using a generalized analytical framework that operates under the assumption that far-field pressure influence of CO2-injection (i.e. a two-phase flow process) can be described in the same way as one would in the case of single-phase flow. Among other leakage pathways, the study also considered leakage through abandoned wellbores. The authors found that accumulation of leaked brine is largest in aquifers directly above or below the target aquifer, with leakage rates through abandoned wells decreasing significantly with larger vertical distance from the injection formation. Their results also confirmed the expectation that the rate of brine leakage from the injection formation increases as the spatial density of leaky wells increases. For simulations where leaky wells were modelled as open wellbores (with a permeability of 1010 mD), this dependency showed asymptotic behaviour: at high spatial densities a further increase in the amount of leaky wells did not cause a corresponding increase in leakage. This is because the conductivity of the aquifer becomes the limiting factor in transport of fluids to the leaky well. For simulations considering leakage through the casing-hole annulus, no asymptotic behaviour was observed when the effective permeability of leaky wells was on the order of 105 mD or lower. Celia et al. (2011) used a highly computationally efficient semianalytical reservoir simulator to assess the risk of leakage through existing wells in the Alberta Basin in Canada, making use of the welldocumented locations and penetration depths of the wells in the area of interest. Little to no quantitative data were available on effective permeabilities of abandoned wells (either in Alberta or anywhere else) at that time, causing the authors to have to rely on an estimated bimodal distribution to populate their simulator with effective well permeabilities in a stochastic framework. Since then, specific field measurements have been made to quantify the range of effective permeabilities that can be expected in abandoned wells (see e.g. Crow et al. (2010), Duguid et al. (2013), Gasda et al. (2013), Kang et al. (2015), Tao and Bryant (2014)), with general agreement between studies that the effective permeabilities span a range from nanodarcies to hundreds of millidarcies. With these new data in hand, the present work employs the same semi-analytical simulator used in Celia et al. (2011), with the aim to provide more general statements on the risk of leakage through abandoned wellbores in mature sedimentary basins. To that end, we consider an idealized layered geology, and assess the amount of CO2 leaking out of the target formation during the injection period. Specifically, we investigate how this amount depends on (1) the spatial density of nearby abandoned wells, (2) the effective permeability of these wells, and (3) variations in geothermal characteristics. Combining results from numerical experiments with the spatial

Fig. 1. Schematic overview of the contributions of different trapping mechanisms, and their evolution over time. Note that formation properties and various other factors can significantly alter the percentages shown. (Reprinted from IPCC (2005).)

have been identified, with many valuable contributions towards effective risk assessment made by previous studies. In the rest of this section, we will highlight some of these recent advances. Firstly, given that many sealing formations have a low but nonzero permeability (Li et al., 2005); the combination of the large areal extent of a typical CO2 plume, high injection pressures, and low CO2 phase viscosity has the potential to cause significant leakage from the target formation. However, given the fact that caprocks are usually water-wet and capillary entry pressure for CO2 is very large, advective leakage of CO2 via the intact caprock is very unlikely. Previous studies on the long-term fate of CO2 in deep saline aquifers found diffusive transport to be slow enough to be of lesser concern (Lindeberg et al., 2003; Lu et al., 2009). A review on leakage through intact caprock can be obtained in the relevant sections of Song and Zhang (2013). A second potential leakage mechanism is flow along permeable regions associated with non-sealing fault zones or fracture networks in the cap rock. Although an ideal injection target will not include any transmissive faults, the large areal extent of the CO2 system and the ubiquitousness of faults and fractures in sedimentary basins (Gluyas and Swarbrick, 2004) makes a thorough understanding of their leakage potential important. Earlier simulation studies of CO2 leakage along permeable fault zones showed complex, non-monotonous leakage behaviour due to phase transitions and non-isothermal effects (Pruess, 2005, 2008), indicating that fully capturing the relevant physics may be important to accurately assess the risks associated with this leakage mechanism. A more detailed review of recent advances regarding leakage through fault zones can be found in the relevant sections of Song and Zhang (2013). In general, the authors find that many questions remain in this particular area of leakage risk assessment. Song and Zhang highlight the uncertainties in both the spatial distribution of fault zones and fracture networks as well as their characteristics, uncertainties that may be difficult to reduce. They also note complexities that may arise from self-limiting or self-enhancing processes in these leakage flows, e.g. through clay shrinkage or precipitation of minerals. This study is concerned with another type of concentrated leakage pathway, namely leakage that may occur along pre-existing wellbores that were drilled through the caprock. Because of the sheer scale at which the drilling campaigns of the upstream petroleum industry have penetrated sedimentary basins across the globe, contacting existing 165

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

density of abandoned wells found in reality and the data on the ranges of effective well permeabilities that is now available, we aim to quantify the likelihood that leakage through plugged and abandoned wells will form a significant risk to the viability of a CCS operation. The structure of the paper is as follows. First, we present a brief overview of the numerical simulator used for numerical experiments. Next, we describe the recent specific field measurements of leaky well permeability. This is followed by a description of the geological model used as a simulation domain, including the way we vary the parameters of interest. We then present some detailed results of a typical simulation, including pressure responses in the different aquifers, evolution of the CO2 plume over time, and an analysis of where leaked brine and CO2 end up. Next, we present the full dataset that shows how leakage depends on the parameters under investigation, and we end by discussing the results and drawing general conclusions about the possible risk of leakage through abandoned wells.

Nordbotten and Celia (2006b) derived a similarity solution for the twophase system, eliminating the need for a full numerical multiphase flow simulation and greatly improving computational efficiency. The solution involves the location of the CO2–brine interface h(r, t) behind which brine is at residual saturation, with r denoting radial distance from the injector and t the time since injection commenced. If the injection stream consists of dry CO2, the residual brine trapped behind the CO2–brine interface will eventually evaporate, leaving only CO2 in the pore space. This region is bounded by the drying front i(r, t). The solutions developed by Nordbotten and Celia (2006b) solve for the two fronts h(r, t) and i(r, t) as well as the pressure along the bottom of the target formation pbot(r, t), noting that the assumption of vertical equilibrium allows for the pressure to be determined everywhere in the three-dimensional system from pbot(r, t), h(r, t) and the known phase densities. The set of equations they obtain can be rewritten in terms of a single similarity variable, in which the three unknowns are only dependent upon this variable χ, defined by:

2. Computational model The numerical experiments presented in this work are conducted using ELSA (Estimated Leakage using Semi-analytical Analysis), a previously developed reservoir simulator that uses a set of analytical and semi-analytical solutions for CO2 injection into confined aquifers (Nordbotten and Celia, 2006b), a 1D-description of brine and CO2 flow along leaky well segments (Nordbotten et al., 2009), and a description of interface upconing near leaky wells (Nordbotten and Celia, 2006a). The simulator is capable of handling an arbitrary number of alternating aquifers and aquicludes, as well as an arbitrary number of leaky wells. The combined simulation algorithm is described in Nordbotten et al. (2009). A brief recapitulation of the governing equations, key assumptions and solution strategy, paraphrased from Celia et al. (2011), will be included here; for more in-depth information on the component parts the reader is referred to these earlier publications.

=

=

We consider a two-phase system, initially at hydrostatic pressure and filled with brine, into which CO2 is injected at a constant mass rate. The mass balance reads:

S

+

·( q ) = Q

kkr, ( p + µ

(1)

g z)

2 (

b

c )g

kH 2

µb Q well

(4)

with the phase densities for brine and CO2 denoted by ρb and ρc [M L−3], respectively, the gravitational acceleration g [L T−2], the absolute permeability of the formation k [L2] and the brine viscosity μb [M L−1 T−1]. When the grouping Γ is small (Γ ⪅ 0.1), the solutions for plume shape and pressure reduce to closed-form expressions, further improving computational efficiency. (For details, see Nordbotten and Celia (2006b).) However, these conditions are not met for any of the numerical experiments conducted in this study, so the resulting ordinary differential equation needs to be solved numerically.

where ϕ is the porosity of the formation with dimension [L0], the fluid density is ρ [M L−3], the fluid saturation is S [L0], time is t [T], and the source term is Q [M L−3 T−1], with subscript α denoting either CO2 (c) or brine (b) phases. The Darcy flux is given by:

q =

(3)

where H [L] denotes the thickness of the aquifer, the porosity of the target formation is ϕ [L0], the vertically integrated residual brine saturation is Sbres [L] and the volumetric injection rate of CO2 is Qwell [L3 T−1]. As a result of this manipulation, the set of equations describing the pressure in the system and the plume shape need only be solved once, providing solutions h(χ), i(χ) and pbot(χ). These variables can then be evaluated for any combination of r and t by determining the corresponding value of the similarity variable χ. In this specific study, residual saturation of brine was set to zero. Consequently, the interface i(r, t) is not considered. The solutions simplify further for conditions where the viscous forces on the injected CO2 dominate the forces associated with buoyancy, quantified by evaluating the dimensionless group Γ defined as:

2.1. Governing equations

t

2 H (1 Sbres) r 2 Q well t

(2)

2.3. Leakage along well segments

2

with permeability tensor k [L ], relative permeability kr [–], dynamic viscosity μ [M L−1 T−1], pressure p [M L−1 T−2], gravity constant g [L T−2] and z [L] positive upward.

Leakage through wells is considered as one-dimensional Darcy flow through aquiclude well segments, i.e. that part of the existing well that connects two neighbouring aquifers. When leakage occurs, the CO2 accumulating in an adjacent aquifer forms a secondary CO2 plume in that aquifer, which in turn can reach more leaky wells, etcetera (Fig. 4). The true mechanisms and flow paths that underlie leakage along wellbores can be complex and varied (as visualized in Figs. 2 and 3 ), but as ELSA is concerned with leakage along relatively large distances, it opts to lump the ability of fluid to flow along a certain aquiclude well segment into an effective permeability coefficient that applies to the bulk of that specific section of leaky well. We note that in this formulation, the leaky well does not form a continuous conduit along its length, but rather a set of one-dimensional connections between neighbouring aquifers (Fig. 4, marked in green). Physically, this is consistent with leakage on the outside of the well casing, either through the cement and drilling mud residue found in the annulus or through

2.2. CO2 plume evolution The core of ELSA is a solver for injection of a constant mass flow of CO2 into a homogeneous, horizontal, confined aquifer. The solver assumes strong buoyant segregation due to the density difference between the two phases, as well as vertical equilibrium with respect to the pressure field in a given formation. Additionally, the capillary transition zone is assumed to be small relative to the height of the formation, allowing the boundary between the CO2 and brine phases to be represented by a sharp interface. Although the viscosity and density of both brine and CO2 can have different values in each formation, fluid properties within any single formation are assumed uniform and constant for the length of the simulation. Under these assumptions, 166

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

along the well with dimension [L T−1], k the effective permeability of the well segment [L2], z positive upwards and kr,α the relative permeability with respect to phase α, taken to be a linear function of saturation. As flow is one-dimensional with a single value for effective permeability characterizing the whole length of the aquiclude well segment, (5) can be rewritten in discrete form as:

q

,i

=

ki k r, µ

,i

p

p ,i

,i +

li

+

g

(6)

with the subscript i referring to the ith aquiclude in the vertical sequence and i+ and i− referring to the values found at the top and bottom of that aquiclude, respectively. The effective permeability is assumed to describe the flow as if it occurs over the entire wellbore area Aw [L2], which means that the volumetric flow rate of phase α along an aquiclude well segment is given by qα,iAw [L3 T−1]. From the point of view of the mass balance, the amount of mass in the aquiclude well segments is assumed to be negligible.

Fig. 2. Schematic representation of CO2 plume migration and possible leakage mechanisms through an abandoned wellbore. (Reprinted from Bachu and Celia (2013), modified from Gasda et al. (2004).)

2.4. Interface upconing When one of the abandoned wells in the domain is leaking, the outflow of mass from the formation in question causes a local decrease in pressure around that well. If the CO2 plume has not yet reached the well, only brine will leak. If the CO2 plume has reached the leaky well, this pressure drop will make the CO2-brine interface move upwards, locally reducing the thickness of the CO2 plume h(r, t). If the thickness of the plume is still larger than zero, only CO2 will leak. If the thickness of the plume does go to zero, both CO2 and brine will leak and an extra step will calculate the fractional flow of both phases along the well segment in question. For details, see Nordbotten and Celia (2006a).

Fig. 3. Schematic representation of injected CO2 plume into an aquifer of height H, with plume thickness h(r, t) and drying front i(r, t). (Reprinted from Nordbotten and Celia (2006b).)

2.5. Simulator algorithm The components of ELSA as described in the previous sections are put together in a highly computationally efficient numerical simulator, the details of which can be found in earlier publications. See Nordbotten et al. (2009) for a description emphasizing the mathematical formulation, and Dobossy et al. (2011) for a description of the algorithms used. No spatial grid is required due to the analytical solutions employed to describe the pressure field, plume evolution, leakage and the local perturbations around leaky wells. As the equations involved are nonlinear, temporal discretisation is still required. At any time step, the pressure is computed at each well location in each aquifer based on superposition of solutions for individual wells, with coupling between aquifers due to the one-dimensional Darcy flow along aquiclude well segments. This results in a set of algebraic equations with the number of unknowns equal to the number of wells times the number of aquifers, which can be solved to obtain the pressure at each well in each aquifer. Combined with the vertical equilibrium assumption and the analytical solutions for the pressure field around any well, the pressure field can be determined in the entire spatial domain without the need for a spatial discretisation. Next, this pressure field is used to update the flows in the system, including the local effect of leakage on the CO2–brine interface using the upconing solution around leaky wells. Lastly, all the CO2 plume shapes are updated to redistribute mass according to the previously determined pressure and flow fields. With the pressures, flow rates and plume shapes updated, the algorithm can move to the next time step.

Fig. 4. Schematic representation of plume evolution and leakage in a multilayer system, with the primary plume in the lowest aquifer and secondary plumes in the aquifers above. (Adapted from Kavetski et al. (2006).)

zones of the aquiclude that were damaged during drilling. If leakage were to happen on the inside of the casing and the rest of the casing is intact, communication with surrounding aquifers would not occur. As abandoned wells usually contain one or multiple plugs, the flow rate on the inside of the casing will be determined by the effective permeability of the sequence of well plugs. Mathematically, the flow along a well segment across one aquiclude is modelled by the one-dimensional multiphase extension of Darcy's law:

q =

kk r, µ

p z

+

g

3. Field measurements of leaky well permeability Multiple studies have been published recently that aim to quantify the permeability of leaky wells. In this section, we will provide an overview of the results of these studies, and highlight the fact that different measurement techniques provide different information about

(5)

where qα is the volumetric flux of phase α (either brine (b) or CO2 (c)) 167

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Table 1 Summary of available field data on leaky well permeability. Measurement type

Length scale [m] −2

Permeability [mD]

−1

Core plugs CHDT VIT

10 to 10 10−2 to 10−1 0 10 to 101

SCP/SCVF Gas emissions

102 to 103 (between casing shoes) 102 to 103 (estimated well depth)

−2

10 to 10−4 to 0.5–1.0 10−3 to 25–170 10−2 to 10−6 to

the leaky well. A summary of the available measurements and their associated length scale can be found in Table 1. In one of the first efforts to quantify leaky well permeability with real field data, Crow et al. (2010) investigated the long-term integrity of wells in a CO2 rich environment by assessing a 30 year-old well from a natural CO2 production reservoir. Data were obtained from sidewall core analysis, as well as vertical interference testing (VIT), a downhole test which measures vertical hydraulic communication outside of the casing along a certain section of the well, typically a few metres in length. The authors found evidence for CO2 migration along the wellbore in the occurrence of carbonated cement, and estimated an effective permeability of 0.5–1.0 mD for an 11 ft. well section using VIT. The complete set of observations and measurements led them to conclude that the amount of fluid migration was probably small, and that the wellbore barrier system has performed well over the lifetime of the well. Gasda et al. (2013) conducted field surveys of potentially leaky wells using VIT. Taking measurements on sections of 3 m, Gasda et al. (2013) found permeabilities spanning many orders of magnitude, ranging from microdarcies to darcies. In another study, Duguid et al. (2013) used Schlumberger's Cased Hole Dynamics Tester (CHDT), core plug measurements, and VIT to estimate the permeability of the cement outside the well casing. VIT measurements indicated permeabilities up to 170 mD. CHDT measurements (which measure cement permeability on the scale of a few centimetres) gave values of cement permeability ranging from 10−4 to 104 mD. Cement core plug measurements in the laboratory showed lower permeabilities in the range of 10−2 to 1 mD, indicating that the interfaces between casing and cement or cement and formation are more significant to leakage than the quality of the cement itself. Tao and Bryant (2014) estimated the permeability of about 300 wells in six different fields that exhibit sustained casing pressure (SCP) or surface casing vent flow (SCVF). Both SCP and SCVF are associated with flow of gas from the underground formations, across cement bonds, into an intermediate annular space and to the surface, and are indicators of an imperfect cement job. Using a Monte Carlo-style framework to account for uncertainties in parameters such as the depth of the leaking formation, Tao and Bryant (2014) reconstructed expected effective permeabilities for the leaky wells. While the results spanned six orders of magnitude, the large majority of leaky wells had a permeability in the range of 10−2 to 10 mD. High outliers were found with an estimated permeability of 102 to 103 mD. More recently, Kang et al. (2015) combined data on methane emissions from 42 abandoned oil and gas wells in western Pennsylvania with estimated well depths from historical documents to estimate the wells’ effective permeabilities, capturing the combined leakage pathways in and around the wellbore. They found effective permeability values ranging from 10−6 to 102 mD, with the highest estimated permeability being 130 mD. It should be noted that because of the different length scales associated with the various measurement techniques (a few centimetres for core plugs, a few metres for VIT, and the entire length of the well for emissions data), these measurement techniques provide different information about the well. The smaller the length scale, the more

1.0 104 103 10 102

Reported in Duguid et al. (2013) Duguid et al. (2013) Crow et al. (2010) Gasda et al. (2013) Duguid et al. (2013) Tao and Bryant (2014) Kang et al. (2015)

localised the measurement is. Core plugs and VIT provide an effective permeability over relatively small distances which can exhibit a large spread based on local factors such as heterogeneities in the cement or an imperfect bond between cement and casing or cement and rock. Over larger distances, the effective permeability can be viewed as the harmonic average of the set of local permeabilities along the length of the well, with lower values dominating the effective permeability of the combined system. For a model that uses a single effective permeability value along any given abandoned well, the effective permeabilities presented by Kang et al. (2015) provide the most appropriate values for well permeability, and we will use those to guide our modelling study. We also note that the upper limit of values are surprisingly consistent among all of the measurements reported using the different techniques. 4. Numerical experiments Numerical experiments were conducted in an idealized, horizontally layered geometry consisting of a stack of alternating aquifers and aquicludes. Two geothermal settings were considered, denoted as a “warm” setting with high surface temperature and high geothermal gradient, and a “cold” setting with low surface temperature and low geothermal gradient. We note that because no thermal effects are considered and fluid properties do not vary in time or within formations, these only affect the properties of formation fluids as set at the beginning of a numerical experiment. In both geothermal settings, situations involving injection into a deep formation as well as into a shallower formation were considered, resulting in four different CCS scenarios. In each scenario, we then varied the spatial density of leaky wells in the domain, as well as their effective permeability. 4.1. Geological model The simulation domain consists of a laterally extensive alternating sequence of 13 aquifers and 12 aquicludes, spanning a depth from 500 m at the top of the shallowest aquifer to 3000 m at the bottom of the deepest aquifer Fig. 5. The lateral extent of the domain is chosen such that the pressure perturbation resulting from fluid injection does not reach the boundaries, resulting in an effectively infinite lateral extent. No-flow boundaries were imposed at the top of the uppermost aquifer and bottom of the lowermost aquifer. Each layer in the domain has a thickness of 100 m, and with the exception of the topmost aquifer, every aquifer was assigned an absolute permeability k = 100 mD and a porosity φ = 0.12. The topmost aquifer was assigned exceptionally high permeability (k = 50 D) and porosity (φ = 0.50), so as to represent an infinite sink where influx of CO2 or brine from lower formations has no effect on pressure. Any CO2 that has leaked into this formation is assumed to continue to move to the shallow subsurface and eventually to the atmosphere. 4.2. Fluid properties The cold scenario was chosen as having a surface temperature Tsurf = 10 °C and a geothermal gradient ΔT = 25 °C/km. The warm scenario was chosen as having a surface temperature Tsurf = 20 °C and a 168

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 5. Schematic representation of the layered geological model considered in this study.

geothermal gradient ΔT = 40 °C/km. Each formation in the stack of aquifers and aquicludes was assigned a single uniform temperature taken as the temperature computed in the middle of that layer using these values for Tsurf and ΔT. An analogous approach was used to estimate the salinity of the resident brine in each formation, assuming a salinity of 0 ppm at the surface and a salinity gradient of 5 × 104 ppm/ km. Lastly, the pressure used to estimate fluid properties was taken as the initial pressure in the middle of each formation, using a hydrostatic pressure gradient of 107 Pa/km. Based on these values of temperature, pressure, and brine salinity, the viscosity and density of both CO2 and brine were computed using empirical correlations found in literature. For CO2 density we used the correlation put forward in Duan et al. (1992), valid at temperatures 273 < T < 1273 K and pressures p < 8 ×108 Pa. For brine density, we used the correlation from Phillips et al. (1981), (283 < T < 623 K, p < 5 ×107 Pa, salinities < 2.6 × 105 ppm). For CO2 viscosity, we used the correlation from Fenghour et al. (1998), (200 < T < 1000 K, p < 3 ×108 Pa). Lastly, the brine viscosity was computed using the relation put forward in Batzle and Wang (1992), valid for T < 523 K and salinities < 4.6 × 105 ppm. For figures showing the variation of fluid properties with depth in both scenarios, please refer to Appendix A; Figs. 15–18 .

Fig. 6. Plan view of the simulation domain for one particular configuration of leaky wells. Injection well in blue, leaky abandoned wells in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

well to decrease dramatically. These three factors cause the leakage estimates computed for each well configuration in this study to be conservative. 5. Results 5.1. Characteristics of leakage behaviour To investigate the basic characteristics of leakage behaviour in our multilayered geology, we consider a simple scenario with a single leaky well, located at a distance r = 5000 m from the injector. The leaky well is assigned an effective permeability kw = 100 mD along its entire length. At t = 0, we begin injection of pure CO2 at a rate of Qinj = 5.0 Mt/year (159 kg/s) into one of the shallower aquifers, L11. In this section, we will first explore when and to what extent leakage of brine and/or CO2 occurs, followed by a closer inspection of where the leaked reservoir fluid migrates. Next, we explain this behaviour by identifying the character and magnitude of the driving forces that promote leakage across the aquiclude well segments. Lastly, we discuss to what extent the observed results can be generalized towards more complex configurations of leaky wells in the system.

4.3. Configuration and properties of wells In addition to the injection well located in the centre, the simulation domain was populated with a number of abandoned wells, the spatial density and effective permeability of which are the main variables under investigation. All leaky wells have a radius rw = 0.1 m. The leaky wells are distributed on an evenly spaced square grid measuring 40 × 40 km in extent, using an even number of rows and columns so as not to have a leaky well at the same location as the injector (Fig. 6). The total number of leaky wells in the domain was increased roughly linearly from nw = 196 (a 14 × 14 grid, 0.12 wells/km2) to nw = 4096 (a 64 × 64 grid, 2.56 wells/km2). Every leaky well in the domain reaches all the way to the bottom of the domain, i.e. the bottom of L25. The effective permeability of wells was taken as constant along the length of each well, as well as having the same value for all wells in a particular simulation. Its value was varied in steps of an order of magnitude from 0.1 mD to 104 mD. We note that in real-life scenarios, not all existing wells will have been drilled to the same vertical depth. Instead, the number of wells penetrating a particular aquifer will decrease as the depth of that aquifer increases. Also, only a portion of the plugged and abandoned wells in a real-life scenario are expected to be leaky. Furthermore, these leaky wells are likely to exhibit variation in effective permeability along the length of the well, with the presence of a lowpermeability section causing the effective permeability of the entire

5.1.1. Onset and magnitude of leakage Tracking the evolution of the CO2 plume over time (Fig. 7), we observe that the outer radius of the plume increases with t , as described in Nordbotten et al. (2009). Indicative of the fact that buoyancy dominates the flow regime, the inner radius of the CO2 plume (i.e. the radial coordinate where h(r, t) = H) grows at a much slower rate. Flow of brine is almost immediately induced in the leaky well because of the pressure increase associated with injection, which changes to a flow of CO2 at the time when the CO2 plume reaches the radial distance of the leaky well (Fig. 8). We note that the flow rate of CO2 through the leaky well is on the order of 10−5 · Qinj. The sudden increase in flow rate observed around t = 14 years is due to the way the secondary CO2 plume around the leaky well forms in the aquifer above the target formation (L9), and the effect that has on the pressure increase computed in L9 associated with the influx of CO2 in that aquifer. The sudden increase occurs at the time when the plume thickness reaches 169

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 7. Inner and outer radius of the CO2 plume versus time for an injection rate Qinj = 5.0 Mt/year.

Fig. 9. Mass of CO2 leaked into surrounding aquifers versus time for a scenario with one leaky well with kw = 100 mD, at radial distance r = 5000 m. Inset axes display the same quantities as the main axes, but note the change of scale.

Fig. 8. Flow rates of brine and CO2 phases across L24, the aquiclude overlying the target formation, plotted versus time.

Fig. 10. Mass of brine leaked into surrounding aquifers versus time for a scenario with one leaky well with kw = 100 mD, at radial distance r = 5000 m. Lines representing L7 and L15 have a constant value of zero.

the thickness of the formation (h = H) just outside the radius of the leaky well. As the flow rate is very low, the secondary plume in L9 resembles a step function, with the almost vertical CO2 front causing the plume thickness h at the location of the leaky well radius to move from close to zero to h = H within a single time step. As the vertical pressure profile is reconstructed based on the assumption of vertical equilibrium which in turn uses the plume height h, this causes an apparently discontinuous decrease in pressure at the bottom of L9 at the location of the leaky well, increasing the flow rate of CO2.

migrate downwards to aquifers below the target formation (inset, Fig. 9). Again, the onset of CO2 leakage corresponds to the time after which the radius of the plume has reached the leaky well. Tracking the mass of brine leaked into the surrounding aquifers (Fig. 10), we observe that brine begins to leak almost instantaneously into the aquifer directly above (L9) and below (L13) the target aquifer. Very soon after the outer radius of the CO2 plume reaches the leaky well, only the CO2 phase leaks upward and the cumulative amount of brine leaked to L9 plateaus. As for the downward leakage, since the inner radius of the plume does not reach the leaky well in the period considered, only brine contacts the lower aquiclude well segment and brine continues to leak into the aquifer directly below the target formation (L13). No leakage of brine is observed into aquifers further down (L15 and deeper) or up (L7 and shallower), with the change in mass of brine in those layers remaining at a constant value of zero.

5.1.2. Fate of leaked formation fluids If we monitor the amount of CO2 that accumulates in the various surrounding aquifers, we observe that most of the leaked CO2 accumulates in the aquifer directly overlying the target formation (L9), with a smaller amount migrating all the way to the shallowest aquifer (L1), as shown in Fig. 9. Quantities of CO2 two orders of magnitude smaller accumulate in the intermediate aquifers, while no CO2 is observed to 170

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

CO2 storage is in large part due to the assumption of uniform wellbore permeability. The details of the fate of leaked formation fluid will depend on the variation of wellbore permeability along the length of a leaky well. 5.1.4. Lessons learned from these observations The system considered in the previous sections is simplified in the sense that it contains only one leaky well. However, the qualitative characteristics of leakage behaviour are expected to be applicable to systems with more leaky wells, and to systems where the effective permeability of the leaky wells is higher. Firstly, the highest pressure increase will be observed in the target aquifer, causing the largest unintended fluid accumulations to be in the aquifers directly above and below the target aquifer. Secondly, due to buoyancy and the fact that the outer radius of the CO2 plume is much larger than the inner radius, CO2 will be more prone to migrate upwards, with a certain amount migrating all the way to the topmost aquifer. We note that in principle, the pressure increase in the target aquifer can be strong enough to overcome the resistance to downward flow associated with buoyancy of the CO2 phase, especially at locations close to the injector where this might be relevant due to the extent of the inner plume radius. However, this was not observed in our simulations. In cases with a larger number of leaky wells or a higher leaky well permeability, some differences may occur. Firstly, an increase in these two variables will cause more brine and CO2 to leak, possibly inducing more substantial pressure increases in the aquifers directly above and below the target aquifer. This pressure increase can in turn promote leakage to the aquifers over- and underlying those, respectively. Although this pressure increase can in principle cascade from one aquifer to the next by the mechanism described, significant attenuation is expected with each intermediate aquifer as leakage rates are still very small relative to Qinj. A further possible difference is that if leaky wells are present close to the injector, the inner radius of the CO2 plume may also contact leaky wells. This may cause CO2 to migrate downwards to the aquifer below the target formation, although it should be noted that the effect of buoyancy, in this case, will work to reduce the leakage rate.

Fig. 11. Pressure and buoyancy contributions to driving force behind leakage across L10, for a scenario with one leaky well with kw = 100 mD, at radial distance r = 5000 m.

5.1.3. Driving forces promoting leakage As becomes clear from the mathematical description of flow along aquiclude well segments (see Eq. (6)), leakage is driven by a combination of a pressure difference across the aquiclude at the location of the leaky well and a term associated with buoyancy. As an example, let us consider the forces driving flow upward along the aquiclude well segment crossing L10, the aquiclude above the target zone, as visualized in Fig. 11. The pressure and buoyancy contributions visualised in this figure are the values of the first and second of the bracketed terms in Eq. (6), respectively. When the system is at hydrostatic pressure and filled with brine, pressure and buoyancy terms are exactly equal and opposite and no flow occurs. When the pressure at the leaky well increases due to injection, a net driving force is observed, at consistent timing with the start of brine leakage described in the previous section. If the brine that was in place at the location of the leaky well is replaced by the lighter CO2 phase, we observe an increase in net driving force as the buoyancy term changes value. This occurs around t = 6 years. The pressure effect leading to the apparent discontinuity in flow rate around t = 14 years, as described in Section 5.1.1, is also clearly observed. We note that, for this case of upward leakage out of the injection formation, the pressure contribution to the driving force is about an order of magnitude greater than the buoyancy contribution (Fig. 11). Because the rate of leakage is five orders of magnitude smaller than the injection rate, the pressure increase observed due to fluid migration into any of the surrounding aquifers is virtually nonexistent. This means that apart from leakage directly out of the target aquifer, any further fluid migration is driven chiefly by buoyancy. This also provides insight into why, as described in Section 5.1.2, very little accumulation of CO2 is observed into any of the intermediate aquifers between L1 and L9, as the flux in and out of these formations along the leaky well pathways will be practically equal. Due to the no-flow condition imposed at the top of the uppermost aquifer (L1), the CO2 that migrates to this formation will accumulate there. By the same logic, brine is only found to migrate towards the aquifers directly above (L9) and below (L13) the injection formation but not any further, as there is no driving force associated with buoyancy. This result implies that in this model, intermediate aquifers between L1 and L9 will not provide significant secondary storage for CO2, due to the fact that buoyancy will cause the CO2 to continue to move upward. We note that the lack of secondary

5.2. Full results From the previous paragraphs, we know that most of the leaked CO2 will accumulate in the formation directly overlying the injection target, with a smaller amount rising to the shallowest zone. Furthermore, we have observed that the resident brine is unlikely to migrate along large vertical distances from the injection target. As the goal of any CCS operation is to retain the CO2 in the targeted formation, we will consider the leakage variable of interest to be the percentage of total injected CO2 that has leaked out of the target formation, regardless of whether it is subsequently retained in an over- or underlying aquifer. In every simulation, we consider injection of pure CO2 at a rate of Qinj = 5.0 Mt/year (159 kg/s) for a period of 50 years. The “shallow” scenarios entail injection into an intermediate aquifer (L11), whereas the “deep” scenarios consider injection into the lowermost aquifer (L25). The range of spatial densities and leaky well permeabilities simulated was chosen to be especially broad and on the high end of expected values, with the objective of finding some subset of parameter combinations for which significant leakage occurs. In other words, we aim to find out how many wells a CO2 plume has to contact, and in what deteriorated state these wells need to be, before leakage becomes significant. 5.2.1. CO2 leaked as percentage of total CO2 injected Fig. 12 visualize the total amount of CO2 that has leaked out of the target formation after 50 years as a percentage of the total amount of CO2 injected, plotted as a function of the spatial density of leaky wells and their effective permeability kw. For readability, contour lines have been added that represent the same data as the colour scale. From these 171

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 12. Logarithmic surface plots visualizing leakage of CO2 as percentage of total CO2 injected, as a function of the spatial density of leaky wells and their effective permeability kw. (a) Cold, shallow scenario, (b) warm, shallow scenario, (c) cold, deep scenario, (d) warm, deep scenario.

leakage figures, it may be observed that generally, the amount of CO2 that leaks out of the target formation is lower for the cold scenarios than the warm scenarios, as well as lower for the deep scenarios than the shallow scenarios. It should be noted that these differences are relatively small. We also find that the amount of leakage is more strongly dependent on the spatial density of wells on the low end of the range of spatial densities considered. From these leakage results, it becomes evident that before leakage enters the 1% range at the end of injection, the domain of investigation will have to be populated by on the order of a thousand leaky wells (corresponding to a spatial density of around 0.6 wells/km2) with an effective permeability of 103 mD.

the size of the CO2 plume, and the number of leaky wells contacted as a result of that. 6. Discussion To put the results in perspective, we need to anchor our view in several ways. First, we use the recently available data of effective permeabilities in abandoned wells to determine a range of permeabilities that can be expected in real-life scenarios. Second, we use data on the spatial distribution of existing oil and gas wells to determine a similar expected range for the spatial density of leaky wells. Using these two constraints, we can determine an area on the leakage surfaces presented in the previous section to determine a range of expected leakage in real-life scenarios. We will then comment on the question of to what degree these leakage rates are acceptable, before discussing the possible limitations imposed by our solution strategy and suggestions for further inquiry.

5.2.2. Percentage of CO2 leaked scaled by number of leaky wells contacted One of the main differences imposed by both the depth and temperature of injection is the density of the CO2 phase in the target formation, and thereby the volume required by an equal mass of injected CO2. In other words, injection at a mass rate of Qinj = 5.0 Mt/year can lead to varying plume sizes among each of the four scenarios, causing a different number of leaky wells to be contacted by the plume. However, if we scale the percentages plotted in the previous paragraph by the number of leaky wells contacted by the CO2 plume after 50 years (Fig. 13), we find that (1) this scaled leakage remains essentially constant as the spatial density of wells varies, and (2) this scaled leakage varies linearly with the effective permeability of the leaky wells. Importantly, we find that all four scenarios show practically identical results. This indicates that the key difference across the four scenarios is

6.1. Data on spatial well density and permeability As mentioned in Section 1, several studies have been published in the past years that describe effective permeability of leaky wells based on measurements on different length scales: gas emissions data for the entire length of the well, vertical interference testing (VIT) on a scale of a few metres, and cement core plug measurements on the scale of a few centimetres. Given the way leaky wells are represented in our model, 172

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 13. Logarithmic surface plots visualizing leakage of CO2 as percentage of total CO2 injected divided by the number of leaky wells contacted by the CO2 plume after fifty years, as a function of the spatial density of leaky wells and their effective permeability kw. (a) Cold, shallow scenario, (b) warm, shallow scenario, (c) cold, deep scenario, (d) warm, deep scenario.

i.e. one-dimensional flow conduits with a single constant effective permeability, the method used by Kang et al. (2015) to determine effective well permeability based on gas emissions data was found to form the closest match. As an upper bound, the authors found an effective permeability of 130 mD. With regards to the spatial density of leaky wells, we may refer to the map of worldwide drilling density as published in the IPCC Special Report on CCS (IPCC (2005), Fig. 5.27) As the most extreme category is all but absent from the map, we take an intermediate value of the second highest category of spatial well density to be a representative very high value (roughly 1.5 wells/km2). For sake of comparison at a higher resolution, we can expect the well spacing as found in the southwestern corner of the Canadian site described in Celia et al. (2011) to have a very high spatial density, as the area has been subject to especially dense drilling patterns due to the tight oil and gas deposits found in the region (Peachey, 2014; Michael et al., 2009). The well density in this area amounts to around 1.6 wells/km2, consistent with the value obtained from the map of worldwide drilling density. It should be noted that these values all consider the density of wells at the surface, but that with increasing formation depth, the number of wells penetrating that formation decreases.

(DOE) has set a climate mitigation goal to achieve 99% storage permanence (Bielicki et al., 2015), i.e. a total leakage of 1% of the amount of CO2 injected. In an earlier study, Hepple and Benson (2005) have attempted to quantify the largest acceptable annual rate of CO2 seepage back into the atmosphere to ensure CCS projects are both safe and effective. They found that for scenarios with small to moderate amounts of CO2 sequestered globally (tens to hundreds of GtC, in their models), an annual seepage rate on the order of 0.1% of the total injected CO2 per year would ensure the effectiveness of CO2 sequestration as a climate change mitigation tool, while still providing for some safety margin. Another study, (Pacala, 2003), which attempted to set performance requirements with a similar strategy, found that in some scenarios even a seepage rate of 1.0%/year could be effective, although the author was careful to note that safety constraints other than effectiveness in terms of climate change mitigation could then become important. To compare these performance requirements to our results, let us compute the total percentage of CO2 that would have leaked in 50 years, if after every year of injection 0.1% of the total mass of CO2 sequestered since the start of injection leaks from the target formation. After 50 years, this amounts to a leakage of 2.55% of the total amount injected. Using this percentage, we can draw a new contour on our leakage results, indicating the region outside of which the effectiveness of the CCS project will be ensured (Fig. 14). Combining this with a rectangle denoting the range of spatial density and permeability of leaky wells we just determined, we can become quite confident that even in case our model CCS site is very densely populated with very

6.2. Performance requirements Although clear regulations regarding acceptable levels of CO2 leakage have yet to be put in place, the U.S. Department of Energy 173

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 14. Logarithmic surface plots visualizing leakage of CO2 as percentage of total CO2 injected after 50 years, as a function of the spatial density of leaky wells and their effective permeability kw. Red contour marks the boundary of effective storage of CO2 as put forward in Hepple and Benson (2005), yellow hatched region marks the range of expected maximum leakage. (a) Cold, shallow scenario, (b) warm, shallow scenario, (c) cold, deep scenario, (d) warm, deep scenario.

leaky wells, they are unlikely to cause leakage to a degree that will endanger the effectiveness of the CO2 storage operation.

qc =

kkr,c µc

g Aw Nw,t

(7)

where the CO2 phase viscosity μc and the density difference between the two phases Δρ are taken as the values for the aquiclude overlying the target formation, Aw is the area of a single leaky well and Nw,t is the number of leaky wells touched by the CO2 plume. As an example, let us estimate for each scenario the leakage rate for a well configuration of intermediate spatial well density (1.44 wells/km2) and leaky well permeability (k = 100 mD) (Table 2). The percentage presented in the last column is calculated as a fraction of the total mass of CO2 in the target formation at the end of the period of active injection. To estimate leakage rates for different spatial well densities or different leaky well

6.3. Leakage on longer time scales The numerical experiments described herein only consider leakage in the period of active injection. To make qualitative statements on leakage behaviour on longer time scales, we must consider (1) the magnitude of the driving forces promoting leakage, and (2) the degree to which the injected CO2 remains mobile. When injection ceases, the pressure in the target aquifer will re-equilibrate, thereby negating most of the pressure rise that was due to active injection. In Section 5.1.3, it became clear that at a distance r = 5000 m, the pressure rise due to injection contributed around 90% of the driving force promoting leakage. This percentage will decrease with radial distance from the injector, but it can still safely be stated that the driving force promoting leakage of CO2 will decrease considerably when injection ceases. To estimate the leakage rate when the pressure drive has subsided but there still exists a mobile CO2 phase, let us assume that the pressure profile in the domain has reverted back to hydrostatic some time after the period of active injection. The total leakage rate of CO2 across the aquifer overlying the target formation at that time may then be approximated by:

Table 2 Relevant parameters and estimated buoyancy-driven leakage rates for a spatial well density of 1.44 wells/km2 and a leaky well permeability k = 100 mD.

174

Scenario

μc [Pa s]

Δρ [kg/ m3]

Nw,t [–]

qc [m3/s]

qc [mass %/year]

Cold, shallow Cold, deep Warm, shallow Warm, deep

6.23 × 10−5 6.24 × 10−5 3.30 × 10−5 4.21 × 10−5

301 327 583 506

468 348 912 448

6.87 × 10−5 5.54 × 10−5 4.90 × 10−4 1.63 × 10−4

6.4 × 10−4 5.1 × 10−4 2.7 × 10−3 1.1 × 10−3

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

10−4%, consistent with the range found in Celia et al. (2011). The computed average spatial density of leaky wells will obviously change as one changes the area under review, but this is to show that to a first approximation, we can estimate leakage from our figures based on information about the injection target. We note that the formation in question would still be rejected as a CCS target because of its poor injectivity. A more suitable injection target identified by Celia et al. (2011) is the Nisku formation (top depth 1881.8 m, H = 72 m, k = 170 mD, φ = 0.12). Based on the amount of wells penetrating that formation (39 wells), we estimate the spatial density of wells to be around 0.02 wells/km2. This spatial density falls well below the smallest spatial density simulated in our study (0.12 wells/km2), indicating very low expected leakage. The fact that this formation is penetrated by such a small number of leaky wells, even though the spatial density of wells at the surface is very high, is a good example of how the spatial density of leaky wells decreases with depth.

permeabilities, these estimates can be scaled linearly with the appropriate factors. By this estimation, we arrive at mass leakage rates on the order of 10−3 to 10−4%/year, so long as the CO2 plume remains mobile. Assuming for illustrative purposes that the CO2 plume stays mobile indefinitely, these estimated leakage rates mean that it will take on the order of 1000–10,000 years before a further 1% leaks out of the target formation. Regarding the mobility of the CO2 phase, there are two further trapping mechanisms that will act to immobilize CO2 as time progresses. Firstly, if further lateral migration of the CO2 plume occurs it will leave behind CO2 at residual saturation, immobilized by capillary effects. Secondly, CO2 will begin to enter the brine phase by dissolution. As CO2-rich brine is more dense then pure brine, this also provides further storage security. 6.4. Applicability of results in different scenarios

6.5. Model limitations and opportunities for future research

The presented results were obtained at a specific injection rate for a specific geometry with certain hydromechanical characteristics, with wells that penetrate all the way down and have a single uniform effective permeability. To determine to what extent these results will be applicable in other scenarios, we note that key variables in the leakage analysis are (1) the radius of the CO2 plume, (2) the pressure rise due to injection and (3) the number of leaky wells touched by the plume. Regarding the extent of the CO2 plume, we predict that so long as buoyancy dominates over viscous effects (Eq. (4)), the extent of the plume will be controlled by the volumetric injection rate of CO2 Qwell, with the layer height H of the injection formation being of lesser influence. Specifically, we expect the area of the plume (and thereby the number of wells contacted for a given spatial density of leaky wells) to increase roughly linearly with the volumetric injection rate. In our results, we noted that scaling the percentage of CO2 leaked by the number of wells contacted by the plume revealed a linear dependence between these two variables. From this it follows that for a domain containing a relatively large number of evenly distributed leaky wells, both the amount of CO2 leaked and the amount of CO2 injected depend linearly on volumetric injection rate. Consequently, the fraction of CO2 leaked will only depend on Qwell through the pressure increase associated with injection, not by the varying extent of the plume. The pressure increase due to injection will rise with increasing volumetric injection rate, as well as with decreasing permeability of the target formation. The number of leaky wells touched by the plume evidently depends, apart from the extent of the plume, on the spatial density of leaky wells in the target formation. We note once more that this density usually decreases with increasing vertical depth, and that the spatial densities simulated in our numerical experiments are very much in the range of the worst possible scenario that is likely to occur in real CCS projects. From our results, we find that if there is a large number of wells in the domain that are evenly spread out, the amount of CO2 leaked out of the target formation shows a linear dependence on the number of wells contacted by the plume. As a comparison, let us consider the geology and leaky well configuration as presented in Celia et al. (2011). The Nordegg/Banff formation described therein is found at a top depth of 1537.8 m, has a thickness H = 80 m, intrinsic permeability k = 4 mD, porosity φ = 0.10. Based on the amount of wells penetrating the layer (719 wells), we can estimate the spatial density of wells to be around 0.3 wells/km2. Its depth and location make this particular injection target a cold, shallow formation in our four-way description. If we look at the corresponding leakage figure (Fig. 12a), at a spatial density of 0.3 wells/km2 and a spread of leaky well permeabilities ranging from 0.1 to 100 mD, we would expect to see the percentage of injected mass leaked out of the target formation to be somewhere in a range of 0.1 to

There are a few simplifications in our simulations that deserve a closer look. Firstly, a plugged and abandoned well with a length of multiple kilometres is unlikely to have a constant effective permeability all the way up. Although our results show that the effective permeability of the overlying aquiclude well segment will be the controlling factor for leakage out of the injection target, the presence of lowerpermeability sections somewhere along the leaky well may prevent leaked CO2 from moving all the way up the well to the shallower zones. Instead, the CO2 will then migrate into intermediate zones, which then provide secondary storage capacity and decrease the risk of unintended migration of CO2 to shallow aquifers. This is in agreement with the results presented in Celia et al. (2011), where variation in effective permeability along leaky wells returned much lower estimates for leakage to the shallower zones in the system. Another simplification in our model is the fact that we did not update any fluid properties in the course of the simulation. However, we expect a full numerical simulation with varying fluid properties will not provide us with different conclusions than the ones drawn here. Scaling the leakage data by the number of wells contacted reduced the four scenarios to practically identical results. This suggests that the size of the CO2 plume is the variable that dominates variation in leakage among different injection depths and geothermal settings, which in turn is determined by the in-situ density of CO2. Incorporating the effect of pressure increase due to injection will, if anything, increase CO2 density and thereby decrease plume radius, but these effects are expected to be small. Even so, models that take varying fluid properties into account more fully can give more insight into that claim, if only for very simple geometries that are not computationally prohibitive. 7. Conclusions Unlike previous studies investigating CO2 leakage risk through abandoned wells in CCS applications, this study combines real data on effective well permeability with field-scale numerical simulations. This was achieved using a previously developed highly efficient simulator, which offers great flexibility in the vertical sequence of aquifers and aquicludes, as well as the number, locations and properties of leaky wells in the domain. In a vertical stack of 13 aquifers and 12 aquicludes, we considered 50 years of CO2 injection at a mass rate of Qinj = 5.0 Mt/ year, at different injection depths and for different geothermal settings. The domain was penetrated by an evenly spaced grid of leaky wells, the spatial density and permeability of which were chosen as independent variables. All leaky wells reach all the way to the lowermost formation and have a constant effective permeability along the length of the well. We note that this will result in conservative estimates of

175

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

leakage, as leaky wells are likely to exhibit variation in effective permeability along their length, and the presence of a single low-permeability segment will dramatically decrease the effective permeability across the entire well. Furthermore, not all wells are drilled to the same vertical depth, and the number of leaky wells contacted by the CO2 plume will decrease as the injection target formation is chosen deeper in the subsurface. Examining the details of leakage in a simple scenario with a single leaky well, we find that most of the CO2 accumulates in the aquifer directly overlying the target aquifer, with a smaller amount migrating all the way to the shallowest zone under the influence of buoyancy. Quantities of CO2 two orders of magnitude smaller accumulate in the intermediate aquifers. Downward migration of CO2 is not observed because the inner radius of the CO2 plume has not reached a leaky well, but may in principle occur for leaky wells close to the injector as the pressure drive may be strong enough to overcome the resistance of buoyancy to downward migration of CO2. Brine was observed to leak only to the formations directly over- and underlying the target formation, as there is no buoyant drive for the brine phase. From this we conclude that unintended migration of brine across large vertical distances is unlikely. We note that this leakage behaviour is due in part to the representation of leaky wells used in our model, i.e. a set of leakage pathways across each aquiclude, not a single continuous pathway across all formations. Physically, this is consistent with leakage outside the well casing. The estimated amounts of CO2 leaked, evaluated as the mass percentage of total CO2 injected that leaks from the target formation, are in agreement with earlier estimates for a particular CCS site in Alberta, Canada (Celia et al., 2011). More importantly, we arrive at the conclusion that, while not all formations investigated by Celia et al. (2011) are ideal targets for CCS, none of them would be deemed unsuitable because of unacceptable leakage risks associated with leaky wells. Based on the fact that scaling the leakage percentages by the

amount of leaky wells contacted by the CO2 plume produces four practically identical graphs, we conclude that differences in amount of CO2 leaked among injection at different depths or geothermal characteristics of the formation are due chiefly to the varying size of the CO2 plume, driven by the differences in CO2 phase density in the injection target and the constraint of constant mass injection of CO2. With regards to the magnitude of the estimated leakage, we have sought to analyse the results of our numerical experiments by combining them with typical spatial densities of abandoned wells, measurements of effective leaky well permeabilities in the field, and performance requirements put forward in the literature. In our estimation, leakage of CO2 through abandoned wells is unlikely to be a major limitation in storage security of CCS projects. This study only considers leakage during the time of active injection. Because the driving forces promoting leakage through abandoned wells decrease upon cessation of injection, and solubility and capillary trapping become more predominant over time, risk of CO2 leakage through abandoned wells will likely decrease after the injection period. Estimates of buoyancy-driven leakage rates find low rates, on the order of 10−3 to 10−4 mass%/year, so long as the CO2 plume remains mobile. Further research might include full numerical simulations that take into account changing fluid properties, phase changes, thermal effects, or any other added complexity, assuming computational challenges therein can be overcome. However, these will not, in our estimation, provide counterarguments to the conclusions drawn here. The only possible exception to this would be strong geochemical reactions in the cement that increase well permeability with time. Acknowledgements This work was supported by the Carbon Mitigation Initiative at Princeton University and by the Princeton Environmental Institute.

Appendix A. Fluid properties Figs. 15–18

Fig. 15. Variation of CO2 phase density ρc with depth, for cold and warm geothermal settings.

176

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 16. Variation of brine phase density ρb with depth, for cold and warm geothermal settings.

Fig. 17. Variation of CO2 phase viscosity μc with depth, for cold and warm geothermal settings.

177

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al.

Fig. 18. Variation of brine phase viscosity μb with depth, for cold and warm geothermal settings.

References

seepage. Environ. Geol. 47 (4), 576–585. IEA, 2017. Global Energy & CO2 Status Report. Paris, France. IPCC, 2005. In: Metz, B., Davidson, O., de Coninck, H.C., Loos, M., Meyer, L.A. (Eds.), IPCC Special Report on Carbon Dioxide Capture and Storage. Prepared by Working Group III of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK 442 pp. IPCC, 2014. In: Core Writing Team, Pachauri, R.K., Meyer, L.A. (Eds.), Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovenmental Panel on Climate Change. IPCC, Geneva, Switzerland 151 pp. Kang, M., Baik, E., Miller, A.R., Bandilla, K.W., Celia, M.A., 2015. Effective permeabilities of abandoned oil and gas wells: analysis of data from Pennsylvania. Environ. Sci. Technol. 49 (7), 4757–4764 PMID:25768798. Kavetski, D., Nordbotten, J.M., Celia, M.A., 2006. Analysis of potential CO2 leakage through abandoned wells using a semi-analytical model. In: Binning, P.J., Engesgaard, P., Dahle, H., Pinder, G., Gray, W. (Eds.), Proceedings of the XVI International Conference on Computational Methods in Water Resources. Copenhagen, Denmark. Li, S., Dong, M., Li, Z., Huang, S., Qing, H., Nickel, E., 2005. Gas breakthrough pressure for hydrocarbon reservoir seal rocks: implications for the security of long-term CO2 storage in the Weyburn field. Geofluids 5 (4), 326–334. Lindeberg, Erik, Bergmo, Per, 2003. The long-term fate of CO2 injected into an aquifer. In: Gale, J., Kaya, Y. (Eds.), Greenhouse Gas Control Technologies – 6th International Conference. Pergamon, Oxford, pp. 489–494. Lu, J., Wilkinson, M., Haszeldine, R.S., Fallick, A.E., 2009. Long-term performance of a mudrock seal in natural CO2 storage. Geology 37 (1), 35. Michael, K., Bachu, S., Buschkuehle, B.E., Haug, K., Talman, S., 2009. Comprehensive characterization of a potential site for CO2 geological storage in central Alberta, Canada. In: In: Grobe, M., Pashin, J.C., Dodge, R.L. (Eds.), Carbon Dioxide Sequestration in Geological Media – State of the Science: AAPG Studies in Geology Vol. 59. pp. 227–240. Nicot, J.P., Oldenburg, C.M., Houseworth, J.E., Choi, J.W., 2013. Analysis of potential leakage pathways at the Cranfield, MS, U.S.A., CO2 sequestration site. Int. J. Greenhouse Gas Control 18, 388–400. Nordbotten, J.M., Celia, M.A., 2006a. An improved analytical solution for interface upconing around a well. Water Resour. Res. 42 (8). Nordbotten, J.M., Celia, M.A., 2006b. Similarity solutions for fluid injection into confined aquifers. J. Fluid Mech. 561, 307–327. Nordbotten, J.M., Kavetski, D., Celia, M.A., Bachu, S., 2009. Model for CO2 leakage including multiple geological layers and multiple leaky wells. Environ. Sci. Technol. 43 (3), 743–749 PMID:19245011. Pacala, S., 2003. Global constraints on reservoir leakage. In: Kaya, Y., Ohyama, K., Gale, J., Suzuki, Y. (Eds.), Greenhouse Gas Control Technologies. 6th International Conference. Proceedings of the 6th International Conference on Greenhouse Gas Control Technologies, vol. 1. 1–4 October 2002, Kyoto, Japan. pp. 267–272. Pacala, S., Socolow, R., 2004. Stabilization wedges: solving the climate problem for the next 50 years with current technologies. Science 305 (5686), 968–972. Pawar, R.J., Watson, T.L., Gable, C.W., 2009. Numerical simulation of CO2 leakage through abandoned wells: model for an abandoned site with observed gas migration in Alberta, Canada. Energy Procedia 1.1. Greenhouse Gas Control Technologies 9, 3625–3632. Peachey, B., 2014. Mapping Unconventional Resource Industry in the Cardium Play Region: Cardium Tight Oil Play Backgrounder Report. Petroleum Technology

Bachu, S., Celia, M.A., 2013. Assessing the potential for CO2 leakage, particularly through wells, from geological storage sites. Carbon Sequestration and its Role in the Global Carbon Cycle. American Geophysical Union (AGU), pp. 203–216. Batzle, M.L., Wang, Z., 1992. Seismic properties of pore fluids. Geophysics 57 (11), 1396. Bielicki, J.M., Peters, C.A., Fitts, J.P., Wilson, E.J., 2015. An examination of geologic carbon sequestration policies in the context of leakage potential. Int. J. Greenhouse Gas Control 37, 61–75. Carey, J.W., 2013. Geochemistry of wellbore integrity in CO2 sequestration: portland cement–steel–brine–CO2 interactions. Rev. Mineral. Geochem. 77 (1), 505. Celia, M.A., Bachu, S., Nordbotten, J.M., Bandilla, K.W., 2015. Status of CO2 storage in deep saline aquifers with emphasis on modeling approaches and practical simulations. Water Resour. Res. 51 (9), 6846–6892. Celia, M.A., Nordbotten, J.M., Court, B., Dobossy, M., Bachu, S., 2011. Field-scale application of a semi-analytical model for estimation of CO2 and brine leakage along old wells. Int. J. Greenhouse Gas Control 5 (2), 257–269. Choi, Y.S., Young, D., Nešić, S., Gray, L.G.S., 2013. Wellbore integrity and corrosion of carbon steel in CO2 geologic storage environments: a literature review. International Journal of Greenhouse Gas Control 16. The IEAGHG Weyburn-Midale CO2 Monitoring and Storage Project S70–S77. Cihan, A., Birkholzer, J.T., Zhou, Q., 2012. Pressure buildup and brine migration during CO2 storage in multilayered aquifers. Groundwater 51 (2), 252–267. Crow, W., Carey, J.W., Gasda, S., Williams, D.B., Celia, M.A., 2010. Wellbore integrity analysis of a natural CO2 producer. International Journal of Greenhouse Gas Control 4.2. The Ninth International Conference on Greenhouse Gas Control Technologies 186–197. Dobossy, M.A., Celia, M.A., Nordbotten, J.M., 2011. An efficient software framework for performing industrial risk assessment of leakage for geological storage of CO2. Energy Procedia 4. 10th International Conference on Greenhouse Gas Control Technologies 4207–4214. Duan, Z., Møller, N., Weare, J.H., 1992. An equation of state for the CH4-CO2-H20 system: I. Pure systems from 0 to 1000 °C and 0 to 8000 bar. Geochim. Cosmochim. Acta 56, 2605–2617. Duguid, A., Butsch, R., Carey, J.W., Celia, M.A., Chugunov, N., Gasda, S.E., Ramakrishnan, T.S., Stamp, V., Wang, J.Z., 2013. Pre-injection baseline data collection to establish existing wellbore leakage properties. In: Energy Procedia 37. GHGT-11 Proceedings of the 11th International Conference on Greenhouse Gas Control Technologies. 18–22 November 2012, Kyoto, Japan. pp. 5661–5672. Fenghour, A., Wakeham, W.A., Vesovic, V., 1998. The viscosity of carbon dioxide. J. Phys. Chem. Ref. Data 27 (1), 31–44. Gasda, S.E., Bachu, S., Celia, M.A., 2004. Spatial characterization of the location of potentially leaky wells penetrating a deep saline aquifer in a mature sedimentary basin. Environ. Geol. 46 (6), 707–720. Gasda, S.E., Celia, M.A., Wang, J.Z., Duguid, A., 2013. Wellbore permeability estimates from vertical interference testing of existing wells. In: Energy Procedia 37. GHGT-11 Proceedings of the 11th International Conference on Greenhouse Gas Control Technologies. 18–22 November 2012, Kyoto, Japan. pp. 5673–5680. Gluyas, J.G., Swarbrick, R., 2004. Petroleum Geoscience. Blackwell Science Ltd., Malden, MA, USA. Hepple, R.P., Benson, S.M., 2005. Geologic storage of carbon dioxide as a climate change mitigation strategy: performance requirements and the implications of surface

178

International Journal of Greenhouse Gas Control 84 (2019) 164–179

T.J.W. Postma, et al. Alliance Canada, Calgary, AB, Canada. Phillips, S.L., Igbene, A., Fair, J.A., Ozbek, H., 1981. A Technical Databook for Geothermal Energy Utilization. Lawrence Berkeley Laboratory, Berkeley, CA, USA. Pruess, K., 2005. Numerical studies of fluid leakage from a geologic disposal reservoir for CO2 show self-limiting feedback between fluid flow and heat transfer. Geophys. Res. Lett. 32 (14). Pruess, K., 2008. On CO2 fluid flow and heat transfer behavior in the subsurface, following leakage from a geologic storage reservoir. Environ. Geol. 54 (8), 1677–1686.

Song, J., Zhang, D., 2013. Comprehensive review of caprock-sealing mechanisms for geologic carbon sequestration. Environ. Sci. Technol. 47 (1), 9–22 PMID:23020638. Tao, Q., Bryant, S.L., 2014. Well permeability estimation and CO2 leakage rates. Int. J. Greenhouse Gas Control 22, 77–87. Zhang, M., Bachu, S., 2011. Review of integrity of existing wells in relation to CO2 geological storage: what do we know? Int. J. Greenhouse Gas Control 5 (4), 826–840.

179