A contribution to risk analysis for leakage through abandoned wells in geological CO2 storage

Advances in Water Resources 33 (2010) 867–879

Contents lists available at ScienceDirect

Advances in Water Resources j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / a d v wa t r e s

A contribution to risk analysis for leakage through abandoned wells in geological CO2 storage Andreas Kopp c, P.J. Binning a, K. Johannsen b, R. Helmig c, H. Class c,⁎ a b c

Institute of Environment and Resources, Technical University of Denmark, Bygningstorvet, 2800 Kongens Lyngby, Denmark Parallab, Bergen Center for Computational Science, University of Bergen, Hoeyteknologisenteret, Thormoehlensgate 55, 5008 Bergen, Norway Department of Hydromechanics and Modelling of Hydrosystems, Universität Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany

a r t i c l e

i n f o

Article history: Received 16 October 2009 Received in revised form 4 May 2010 Accepted 4 May 2010 Available online 15 May 2010 Keywords: CO2 storage Multiphase flow Risk assessment Leakage

a b s t r a c t The selection and the subsequent design of a subsurface CO2 storage system are subject to considerable uncertainty. It is therefore important to assess the potential risks for health, safety and environment. This study contributes to the development of methods for quantitative risk assessment of CO2 leakage from subsurface reservoirs. The amounts of leaking CO2 are estimated by evaluating the extent of CO2 plumes after numerically simulating a large number of reservoir realizations with a radially symmetric, homogeneous model. To conduct the computationally very expensive simulations, the ‘CO2 Community Grid’ was used, which allows the execution of many parallel simulations simultaneously. The individual realizations are set up by randomly choosing reservoir properties from statistical distributions. The statistical characteristics of these distributions have been calculated from a large reservoir database, holding data from over 1200 reservoirs. An analytical risk equation is given, allowing the calculation of average risk due to multiple leaky wells with varying distance in the surrounding of the injection well. The reservoir parameters most affecting risk are identified. Using these results, the placement of an injection well can be optimized with respect to risk and uncertainty of leakage. The risk and uncertainty assessment can be used to determine whether a site, compared to others, should be considered for further investigations or rejected for CO2 storage. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Carbon dioxide storage in subsurface reservoirs is a promising technology for managing anthropogenic carbon dioxide fluxes into the atmosphere. The design of a carbon dioxide storage system is subject to considerable uncertainty. For example, there is some chance that geological conditions will lead to leakage from the CO2 reservoir. But uncertainty is also involved in the engineering work, in possible leakage through manmade pathways back into the atmosphere and in the abandonment of the site. In order to select suitable geological formations for CO2 storage and design these systems it is therefore important to develop a concept of risk. The management of risk is an inevitable part of every CO2 storage attempt to mitigate potential harm to the health of humans and animals, and to the environment. Risk management must be site specific to be successful. Nevertheless, a general framework should be adopted to serve as a basis for site specific risk management. This work addresses the development of a novel risk assessment framework for CO2 storage. As

⁎ Corresponding author. E-mail addresses: [email protected] (A. Kopp), [email protected] (P.J. Binning), [email protected] (K. Johannsen), [email protected] (R. Helmig), [email protected] (H. Class). 0309-1708/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.05.001

a first step, CO2 plume evolution is investigated in a reservoir using numerical simulations based on a large statistical database of potential reservoir properties. The proposed approach to risk assessment is not specific to any site, region or basin. Rather it reflects the entire range of possible reservoir properties that can be encountered. In this study, the entire range of reservoir properties provided by the database is considered. Such an approach can be relevant in the very early project phase, when very little information about the actual reservoir properties is available. The assessment might lead to the conclusion, that a specific reservoir is more suitable for CO2 storage than another. Further investigations of a pre-selected reservoir (seismic investigations, well core analysis, detailed numerical simulations etc), must then be used to determine the suitability of the site or reject it. Hence, this risk assessment serves as an initial analysis method, and is intended to determine the need for additional reservoir investigations. A characteristic of CO2 storage in geologic formations is the influence of several physical and geochemical trapping mechanisms on very different time scales. On a short time scale, structural and stratigraphic trapping are the dominant trapping mechanisms. On the medium time scale, residual and solubility trapping become increasingly important [16,30,33,41,62], while at the long time scale mineral trapping may become the dominant mechanism provided that

868

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

geochemical conditions are favorable [25]. Storage security increases over time, since residual, solubility and mineral trapping represent a long-term fixation of CO2. The contribution of various trapping mechanisms on multiple time scales adds complication to risk assessments. Similarities can be found in radioactive waste disposal where performance assessment calculations must consider the return of radionuclides to the accessible environment over periods of longer than 10,000 years [46]. For CO2 storage, two time scales may be of interest, and these are directly related to the environmental impacts that might occur from leakage: local environmental effects or global effects. The scope of the risk assessment and the associated time scales are defined dependent on these impacts [31]. The medium to long timescale (20 years or longer) considers global effects, that is the return of the CO2 back into the atmosphere. The short time scale considers local environmental impacts, leading to risk directly associated with exposure to leaked CO2; that is local health, safety and environmental hazards. For example, Ref. [59] suggests a potential for self-enhancement of leakage rates, leading to a so-called ‘pneumatic eruption’; although the author says it is unlikely. Ref. [56] states that local risks (together with economic considerations) are likely to constrain allowable leakage rates more tightly than global impacts. Hence, the focus of the risk assessment method developed here is on the short time scale, i.e. up to 20 years. Only the injection period is simulated. Risk is usually defined in terms of three questions [34]: How can a system fail? What is the likelihood of failure? What are the consequences of failure? The evaluation of these questions is sometimes difficult because of the lack of data, poor knowledge of processes, etc. Refs. [5] and [4] investigated natural and industrial analogues (natural gas storage, liquid waste disposal, and nuclear waste disposal) that could be used to guide risk mitigation in subsurface CO2 storage. At present, geological CO2 storage risk assessments rely on expert opinions to rate individual scenarios, which then get combined to form an overall perception of risk. For example, Ref. [55] developed a screening and ranking framework to evaluate potential storage sites on the basis of risk arising from possible CO2 leakage to health, safety and the environment. In their spreadsheet-type analysis, the investigator has to rate, for example, the lithology of the storage formation by giving it a ‘weight’ (relative importance), an ‘assessment attribute relative to risk’ and a ‘certainty factor’ (judging how well the information is known). Estimates of the reliability and uncertainty are dependent on the level of knowledge of the investigator. The risk assessment outlined in the following, decreases subjectivity and therefore may increase confidence of investigators/experts using such comprehensive frameworks. The rather simple method presented here tries to evaluate risk, associated with CO2 leakage based on the uncertainty of reservoir properties. The first question addressed is ‘how can a storage system fail’. A CO2 reservoir can fail by leakage through fractures, wells or other geological weaknesses located at some distance away from the CO2 injection well. We focus here on CO2 leakage through abandoned wells, although the method can be extended naturally, for example, to fractures. Leakage through an abandoned well can occur by multiple pathways [8]. Since information is limited on the state of such a well (material, condition and set-up of cement plugs, casings, etc.), uncertainty of potential leakage is large. The leakage through one well is then multiplied by the number of wells encountered by the CO2 plume in the subsurface, yielding cumulative leakage out of the storage formation. Leakage rates have been investigated by Refs. [51], [53], and [52] by a semi-analytical approach. Ref. [23] analyzed the spatial characteristics of well locations in the Alberta basin (Canada) and states that a typical CO2 plume can encounter up to several hundreds of wells in high-density areas. Here, leakage is judged to be significant if it occurs within a given time period. Such a failure can be evaluated by numerical simulation. The numerical simulator (see

Section 4) is used to determine the distribution of CO2 in the reservoir, and failure is defined to occur if leakage is produced. This is equivalent to stating that in reservoirs where the CO2 is spread over a greater lateral area it is more likely that leakage occurs than in reservoirs with a compact CO2 volume. The second part of a risk assessment is to determine the likelihood that such a failure will occur. This is assessed here by considering the potential range of reservoir parameters, like porosity, geothermal gradient, depth, permeability anisotropy etc. by examining information recorded in a reservoir property database. Here the parameters are assumed to be uncorrelated. As is explained below in Section 2, a computationally costly methodology is proposed that allows random sampling of the parameter space and requires many simulations, for each of which failure is assessed by examining whether CO2 has leaked or not. For completing a total number of 50 large-scale simulations, the ‘CO2 Community Grid’ environment [32] was used. The final part of a risk assessment is to determine the consequence of failure, or the damage. In this paper, damage is defined to be equivalent to the potential mass of CO2 leaking out of a reservoir. Here the potential leakage mass is defined to be the total amount of CO2 that has passed by a leaky well at a given radial distance away from the injection well at a given time. This is equivalent to saying, that all CO2 that reaches a leaky well is damaging. Within the given framework, this can be considered a conservative approach to risk estimation (see Sections 6 and 7). For this risk analysis, only homogeneous and radially symmetric cases are studied, the proposed approach can easily be extended to more complex scenarios. However, their interpretation would be more difficult. The approach to risk assessment presented in this paper is quite general and can be applied using either numerical or analytical models of the distribution of CO2 in repositories. Here a numerical approach is employed, in part so that the method can be easily extended to more complex cases in subsequent work, and in part to demonstrate how advances in computing technologies can be applied to engineering problems. 2. Risk analysis concept Risk can be defined in many ways. According to Ref. [34] it can be defined in terms of three questions: • How can a system fail? • What is the likelihood of failure? • What are the consequences of failure? For the first question, there are many ways in which a CO2 storage system can fail. For example the system can fail to provide the necessary injectivity with respect to the CO2 mass delivered by the surface installations (power plant, pipe network, process engineering devices etc.). There are a number of engineering failures associated with the design of the installations and the injection process, for example well corrosion, untight well plugs, and formation clogging due to halite precipitation in the close vicinity of the injection well. There might be failure associated with the design of the CO2 monitoring devices to measure CO2 plume evolution and possible CO2 leakage. The system can also fail to provide the necessary capacity to store the intended CO2 production. Finally, the system can fail after injection due to CO2 leakage, e.g. through the injection well. The risk study conducted here focuses on CO2 leakage out of the storage reservoir. Ref. [31] summarized number of possible leakage pathways: (i) through the pore system in low-permeability caprocks such as shales, if the capillary entry pressure at which CO2 may enter the caprock is exceeded; (ii) through openings in the caprock or fractures and faults; and (iii) through anthropomorphic pathways, such as poorly completed and/or abandoned pre-existing wells. We focus on leakage through poorly completed and/or abandoned pre-

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

existing wells although the approach presented can also be applied to the other leakage pathways. Abandoned wells have been identified as one of the most probable leakage pathways for CO2 ([4,23]). A preexisting (leaky) well is assumed to exist at some distance from the injection well. It has a constant diameter of 1 m and is screened over the full reservoir thickness. The leakage pathway in the abandoned well is not explicitly modelled, as will be explained below in the paragraph on the consequences of failure. The likelihood of failure can be assessed by considering the potential range of reservoir parameters like porosity, geothermal gradient, depth, and anisotropy. The statistical characteristics of reservoir parameters were determined from a database including information from over 1200 reservoirs or by a model (permeability anisotropy). Based on this data, probability density functions (pdf's) have been derived. The parameter space of potential reservoir properties is then randomly sampled and simulations are conducted to assess the distribution of CO2 in the subsurface. Correlations between the parameters given in the database were investigated and are very low [39]. For each simulation, failure is assessed by examining whether CO2 has spread to a given radius within a given time. If CO2 has spread beyond the leaky well radius, this case has failed. By simulation of many such cases, a likelihood of failure can be determined by relating the number of cases that failed nf (dependent on distance r from the injection well to the leaky well and time t) to the total number of cases simulated (N). To determine the consequences of failure, we define the damage to be the amount of CO2 mass that has spread within the leaky well sector beyond the leaky well distance (compare to Fig. 1). The assumption is that damage (D) measures the mass that can potentially leak through the pre-existing well out of the system. This means that in this study the leakage process itself is not modelled. Ref. [31] outline a framework for assessing environmental risks. Two categories of potential environmental impacts from geological CO2 storage are defined: local environmental effects and global effects. Global effects arise from CO2 leakage to the atmosphere and relate to the uncertainty in the effectivity of CO2 storage. Local effects arise, for example, from the direct exposure of plant and/or life species to CO2 close to the site. This study does not distinguish between these impact categories. The impact of leaking CO2 is not explicitly considered. Leaking CO2 could possibly (i) reach the next shallower geological unit where it is safely stored, (ii) pollute a fresh-water aquifer, or (iii) leak back to the atmosphere. One should be aware that for an individual leakage scenario, the real damage may not be proportional to the amount of CO2 that has leaked.

869

Since a simple reservoir geometry is used in this study together with a random sampling of statistical parameter distributions, this risk assessment does not refer to any specific site. Rather risk, as it is defined here, refers to a statistical leakage probability that might occur at a preexisting well. Risk is defined by multiplying the likelihood of failure with the consequence of failure, which is expressed in Eq. (1).

Risk ½kg =

nf N |{z}

likelihood of failure

⋅

N D ∑ i; i = 1 nf |{z}

ð1Þ

consequence of failure

where nf [–] is the number of cases that failed (dependent on distance r [m] from the injection well to the leaky well and time t [d]), N [–] is the total number of cases simulated, and Di is the damage [kg] of case i at time t. The calculation of the damage Di is done in two steps. First, all mass of CO2 that passes the distance r from the injection point is d determined. This mass of CO2 is then multiplied with in order to 2πr count only that fraction of CO2 mass that swept over the leaky well of diameter d (m). Risk has the unit of mass (kg). For this study, the diameter of the leaky well is 1 m (as indicated in Fig. 1). Since leakage is not explicitly modeled (the leaky well is not included in the simulations), a single model run can be used to calculate damage and risk for many leaky wells radii r at different distances r. This allows a flexible and efficient evaluation of risk. The concept outlined above, calculates risk as an average over the entire parameter space. This becomes apparent because nf cancels out D of Eq. (1) and what is left is average risk ∑Ni = 1 i . The average risk, N using the statistical distribution within the parameter space, can be used to motivate important conclusions (as given in Sections 5 and 7), but does not determine the risk of an individual case. A procedure to calculate risk for an individual case, could be the following: Assume that each simulation has four independent input parameters, call them p1…p4, coming from the parameter distributions. When fixing p2, p3, and p4, the probability of failure P1 due to variation of parameter p1 can be determined by varying p1 and finding p1⁎ that partitions between failure and no failure. The probability P1 is then determined from the probability distribution to be the probability that p1 N p1⁎, where p1 N p1⁎ defines the range of parameters p1 leading to failure. The total failure probability P of this case can then be calculated by P = P1·P2·P3·P4. The risk associated with this case would be the damage produced by this case times the failure probability, i.e. Riski = Di·P1·P2·P3·P4. This individual risk is of course dependent on the distance r from the injection well to the leaky well and the time t. This study does not consider the risk of an individual case any further, but instead focuses on the average risk and the conclusions that can be drawn from it. 3. Parameters

Fig. 1. Sketch of the radially symmetric model domain and definition of the leakage that occurred after time t at a leaky well in distance rleakywell from the injection well. CO2 is injected into the center boundary and spreads laterally in the reservoir. Buoyancy forces drive the CO2 upwards where it accumulates below the caprock. Consequently the CO2 plume spreads faster at the top than at the bottom.

To define parameter input sets for simulations, four independent parameter distributions are randomly sampled. These four parameters are porosity φ, depth of the reservoir below surface D, average geothermal gradient at the site denoted dTdz, and the anisotropy between vertical and horizontal absolute permeability AnIso. These are called primary parameters. All other parameters, called secondary parameters, are functions of primary parameters or constants. The primary parameters have been selected based on previous sensitivity investigations [36] where they are shown to be the parameters which are most influential on CO2 plume evolution behavior in a reservoir. Other parameters are not selected as primary parameters for three reasons, (i) lack of data, (ii) correlation with the primary parameters, and (iii) computational effort would be too high. Detailed information on statistical characteristics of primary parameters, dependencies of secondary on primary parameters, mutual

870

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

interrelation of parameters, parameter sources, etc. is given in the following and summarized in Table 2. 3.1. Primary parameters Statistical characteristics of the primary parameters φ, D, and dTdz can be calculated from the U.S. National Petroleum Council public database [54]. The NPC database is part of the TORIS database, which contains data for over 2500 crude oil reservoirs in the US. TORIS is unfortunately not available to the public. Nevertheless, the public NPC database includes data from over 1200 reservoirs. Detailed information on the statistical analysis of the database is given in [37]. Table 1 summarizes the resulting statistical characteristics. In Fig. 2 the histograms of the NPC parameter distributions employed in the following are shown. The hypothesis that datasets follow either a normal or log–normal distribution (shown in Fig. 2 as lines), is rejected by a Kolmogorov–Smirnov test [43]. Consequently, standard probability distributions cannot be used. To calculate statistical characteristics of the 4th primary parameter, the absolute permeability anisotropy (AnIso), a layered reservoir is considered. Layers have varying thickness and an alternating constant absolute permeability of 3.24ot10−15 m2 and 1451.85·10− 15 m2. These permeabilities reflect the 5th and 95th percentiles of the absolute permeability distribution in the NPC database. When estimating upscaled or effective permeability for flow in a layered reservoir, the direction of the flow is of importance. The harmonic mean (see Eq. (2)) of the permeabilities is used to estimate effective permeability for flow perpendicular to the layers.

HM =

∑i mi ∑i

ð2Þ

mi Ki

where i is the layer index, mi [m] is the thickness of the layer having permeability Ki [m2]. The arithmetic mean (see Eq. (3)) of the permeabilities is used to estimate effective permeability for flow parallel to the layers,

AM =

∑i mi ⋅Ki : ∑i mi

ð3Þ

The anisotropy of effective flow permeability is calculated as the harmonic mean divided by the arithmetic mean. When varying the number and the thickness of the layers, the distribution shown in Fig. 2 (bottom-right) is the result. Statistical characteristics of the distribution are: minimum = 0.0089, maximum = 1.0, arithmetic mean = 0.0283, median = 0.0118, 5th percentile = 0.00891 and 95th percentile = 0.0842. Section 3.3 gives the procedure of sampling the primary parameter distributions in order to define a simulation case.

Table 1 Statistical characteristics for parameters in the NPC-database: porosity (φ), average geothermal gradient (dTdz) and depth below surface (D). Stated is the number of datasets used (n), minimum (Min), maximum (Max), arithmetic mean (A. Mean), median, 5th and 95th percentile of the data. Note that random sampling of the primary parameters occurs only in the range between the 5th and 95th percentile of the data, so that unrealistically high values of porosity and geothermal gradient are not relevant here. Parameter

n

Min

Max

A. Mean

Median

5th Perc. 95th Perc.

φ [%] 1222 7 58 21 20 9 dTdz [K/m] 1250 0.009 0.298 0.036 0.03 0.018 D [m] 1273 17 5502 1680 1524 386

34 0.062 3495

3.2. Secondary parameters Secondary parameters, for example absolute horizontal and vertical permeability, capillary pressure, residual phase saturations and relative permeabilities, are dependent on primary parameters. Functional dependencies are given in the following. In Fig. 3, correlation functions from the literature (Carman-Kozeny type [28,57,64]) along with NPC database values of porosity φ and absolute horizontal permeability Kh are given. Correlation coefficients of functions and database values are rather poor. However, the Carman-Kozeny function given in Eq. (4) with an average grain diameter of 35 μm fits the interpolated NPC database values reasonably well. The absolute horizontal permeability is therefore calculated from porosity by using: 3

Kh =

1 φ d 150 50 ð1−φÞ2

ð4Þ

where Kh is absolute horizontal permeability [m2] and d50 is average grain diameter [μm]. Absolute vertical permeability Kv [m2] is calculated by multiplying the absolute horizontal permeability Kh by the anisotropy factor AnIso (a primary parameter). Calculation of residual water saturation Sw,r as a function of absolute permeability and porosity is given by Ref. [28] for a specific site analyzed in this study. However, the permeability and porosity ranges and distributions leading to the functionality given in Ref. [28] are not identical to the data considered here. This results in unrealistic values for residual water saturation Sw,r. Since no other correlation function is known to the authors, residual water saturation is defined as constant. Likewise, residual gas saturation SCO2,r is constant. Relative water and gas permeability kr [–] is defined by a Brooks&Corey model ([7]). Since Brooks&Corey input parameters Sw,r, SCO2,r together with the sorting factor λ are constant, the relation is identical in all simulations. This is a simplification which is necessary due to lack of better knowledge (data). Only very few measured relative permeability relations of CO2-brine systems are known, e.g. [2]. But this sparse data does not allow the definition of dependencies on e.g. primary parameters. Relative permeability has a high influence on plume evolution behavior, as will be discussed in Section 7. Capillary pressure pc [Pa] is dependent on porosity, interfacial surface tension between the water and the CO2 phase, water-rock contact angle, and water saturation. A measured capillary pressure relation and the Leverett J-function are used to determine a distribution of capillary pressure relations as a function of the primary parameters. The measured data is from [58] Experiment no. 12 where the sand column has a permeability of approx. 2·10− 10 m2 and a porosity of 0.37. The experiment was conducted at a pressure of 8.5 MPa and at 300.15 K. The author also gives values for interfacial tension and contact angle. These data allow to use a Leverett J-function [44] to normalize the capillary pressure relation. This relation is representative for the investigated rock-type. The dimensionless Leverett J-function is given in Eq. (5): p ðS Þ JðSw Þ = c w σ⋅cosΘ

rffiffiffiffiffiffi Kh : Φ

ð5Þ

N m

where pc is capillary pressure [Pa], σ is interfacial surface tension [ ], and Θ is the contact angle [°]. The assumption, underlying Eq. (5) is, that the porous medium can be modelled as a bundle of non-connecting qffiffiffiffi ffi capillaries. In this model,

Kh Φ

[m] is proportional to the typical pore

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

871

K

Fig. 2. Histograms data show relative frequency of porosity [%] (top-left), geothermal gradient [ ] (top-right) and reservoir depth below surface [m] (bottom-left) derived from the m NPC database and relative frequency of the natural logarithm of anisotropy between vertical and horizontal intrinsic permeability [–] (bottom-right) derived from a anisotropy model. Lines indicate normal distributions having the same statistical characteristics as the respective histogram datasets.

throat radius. The Leverett J-function is used to scale the given capillary pressure relation to other values of permeability, porosity, temperature and pressure (since the interfacial tension σ changes with pressure and

Fig. 3. Correlation functions between porosity and absolute permeability and NPC database reservoir values. Given are correlations by a Carman-Kozeny model ([64]) for different average grain diameters, correlations by Ref. [57] for average sandstones and Rotliegend sandstone and a correlation given by Ref. [28].

temperature). This is done by equating the J-function for two cases and solving for the capillary pressure of interest (c.f. Eq. (6)). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Kh;1 Φ2 σ2 cosΘ2 pc;2 ðSw Þ = pc;1 ðSw Þ : Kh;2 Φ1 σ1 cosΘ1

ð6Þ

Subscript ‘1’ indicates the reference capillary pressure relation and rock properties (here experiment 12 in [58]), and subscript ‘2’ the capillary pressure relation for the reservoir of interest. An equation for the interfacial tension σ for water and CO2 systems, as a function of pressure and temperature, is given by Ref. [42] (see Fig. 4). Examples of resulting capillary pressure relations are shown in Fig. 5. Salinity S [kg/kg] is defined as a constant. This is because a previous sensitivity study has shown that this parameter had little influence on plume evolution [36]. However, statistical characteristics could have been calculated from the NPC database. Reservoir dip [°] is defined to be a constant at zero. Although, statistical characteristics could have been calculated from the NPC database, the computational effort required to include reservoir dip as a primary parameter would have been too high. In order to include reservoir dip in simulations, full 3D simulations would have become necessary and the computational advantage of employing radial symmetry would have been lost. Nevertheless, reservoir dip strongly influences plume evolution behavior and results will be discussed with respect to this simplification in Section 7. CO2 fluid properties (density [kg/m3] and viscosity [Pa·s]) are calculated and depend on temperature and pressure. Since in the risk analysis, reservoir depth D is a primary parameter and its range includes depths where CO2 occurs in the gaseous, liquid, and supercritical states, the calculation of CO2 fluid properties needs to cover a broad range of possible pressures and temperatures. A difficulty occurs when conditions change during the simulation, so

872

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879 Table 2 Definition of primary and secondary model input parameters, dependencies, and sources. Parameter type

Symbol Unit

Primary parameters Porosity φ

[–]

G.grad

dTdz

[K/m]

Depth

D

[m]

Permeability anisotropy

AnIso

[–]

Secondary parameters (variable) [m2] Horizontal Kh permeability Vertical Kv [m2] permeability [K] Temperature TBC + Init

Fig. 4. Interfacial tension between CO2 and water dependent on pressure and temperature [42]. Accuracy is reported to be greater than 95% within the experimental range. The experimental range includes temperatures between 278 and 335 K and pressures between 0.1 and 20 MPa. The range included in this work is expanded beyond the experimental range; however interfacial tension does not vary significantly at high pressures and temperatures (reflecting great depth). At shallower depth, pressure and temperature conditions are covered by the experimental range.

that the vapor pressure curve is crossed, and gas and liquid can coexist. In these cases CO2 properties change discontinuously by orders of magnitude. If the conditions change from sub-critical to supercritical, or vice versa, without crossing the vapor pressure curve, CO2 properties change continuously. Other secondary parameters, like brine fluid properties, solubility of CO2 in brine etc. are calculated by functions given in literature (as cited in Table 2) or are defined as constant values since results are insensitive to variation of these parameters. All parameters given in the NPC database are tested for mutual interrelations. Apart from the considered dependencies (between porosity and absolute permeability and between depth and temperature) no significant correlations have been found.

Dependency/ Type

Source/Remark

Statistical Distribution Statistical Distribution Statistical Distribution Statistical Distribution

[54]

f(φ) f(Kh,AnIso) f(D,dTdz]

pBC

[Pa]

f(D,ϱ⁎b )

+ Init XCO2 w

[kg/kg]

f(T,p)

ϱCO2 hCO2 μCO2 ϱb

[kg/m3] [J/kg] [Pa⋅ s] [kg/m3]

f(T,p) f(T,p) f(T,p) f(T,p,S,XCO2 w )

Brine enthalpy Brine viscosity CO2-brine interfacial tension Capillary pressure

hb μb σ

[J/kg] [Pa⋅ s] [N/m]

CO2 f(T,p,S,Xw ) f(T,S) f(T,p)

pc

[Pa]

f(φ,σ,Θ,Sw⁎⁎)

Relative permeability

kr

[–]

f(Sw,r,SCO2,r,λ, Sw⁎⁎)

Pressure CO2 solubility in brine CO2 density CO2 enthalpy CO2 viscosity Brine density

Secondary parameters (constant)

Value

Salinity Water solubility in CO2 Contact angle Residual water saturation Residual CO2 saturation Sorting factor Reservoir dip

[54] Layered reservoir model

Carman-Kozeny model ([64]) Definition Geothermal boundary and initial condition Hydrostatic boundary and initial condition [14] [65] [65] [18] [1,21] [29] [12,14,29,47] [1] [42]

scaled by a Leverett J-function based on a pc relation given in [58] [7]

S Xw CO2

[kg/kg] [kg/kg]

0.048 0.001

[54] median value Negligible

Θ Sw,r

[∘ ] [–]

49.46 0.3

[58] –

SCO2,r

[–]

0.05

–

λ

[–]

5.87

dip

[∘ ]

0

[58] [7] –

Injection CO2 mass influx Heat influx

[54]

Value/ Dependency qCO2 qh

[MtCO2/ yr] [J/yr]

1 f(p,T,qCO2)

Medium-sized coal fired power plant Resembles injection at actual reservoir conditions

.⁎ ϱb is constant for the calculation of initial and boundary conditions (1160 [kg/m3]). ⁎⁎ Sw indicates water phase saturation [–].

3.3. Procedure of defining a simulation case

Fig. 5. Capillary pressure dependence on water saturation and porosity. Here σ1 = σ2, i.e. relations are at constant temperature and pressure.

The four independent primary parameters are sampled by generating four random numbers between zero and one. These random numbers reflect percentiles of the parameter distributions, for example if the random number is 0.5, the median of the parameter distribution is obtained. Primary parameters are constant throughout a simulation. Following the generation of four random primary parameters, the secondary parameters are calculated. All secondary

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

parameters, which are dependent on the unknowns in the simulation, i.e. pressure, temperature, saturation or CO2 mass fraction dissolved in brine, are continuously updated during model simulations. This means, for example, that the capillary pressure relation could be different in each gridblock and that it changes with time. Other secondary parameters are fixed during the simulation. 4. Numerical simulation In the simulations, a two-phase two-component non-isothermal model approach is used. Phases include a CO2-rich phase and a waterrich-phase (brine). Components are CO2 and water. Mass transfer between the phases is considered, depending on pressure, temperature and salinity. An energy balance is solved (assuming local thermodynamic equilibrium), which accounts for temperature changes in the reservoir, e.g. due to gas expansion and resulting Joule-Thompson cooling. Salt is not treated as a component, but salinity as a (constant) secondary parameter. An overview of dependencies of CO2 and brine fluid properties and the solubility model used to calculate CO2 dissolution in the water-rich phase is given in Table 2 where literature references are cited. While the fluids are compressible, compressibility of the rock matrix is not considered, since the pressure-induced reduction of porosity and permeability is negligible at great depth if the plume evolution is of interest (not if resulting overpressures are of interest). The numerical simulator used for the study is the research code MUFTE_UG. It has been extensively tested in code intercomparison studies and provides results in good agreement with other commercial and non-commercial codes ([9,35,60]). For a detailed review of the simulation environment MUFTE_UG see Refs. [26] and [10]. Detailed information on the capabilities of MUFTE_UG for simulation of CO2 storage in geological reservoirs is given in [6]. CO2 is injected in a radially symmetric model domain as shown in Fig. 1. The reservoir has a lateral extent (radius) of 100 km and a height of 100 m. A radially symmetric grid is used with 5750 elements. The smallest radius is 1 m. Grid spacing increases towards the outer edge of the domain with 100 elements used for r b 5000 m and an additional 15 elements between 5000 m and 100.000 m radius. A constant vertical spacing of 2 m per grid cell is used. The numerical performance of the simulator is documented in the intercomparison study of Ref. [9]. CO2 is injected at the center boundary at a constant rate of 1 MtCO2 per year. This rate resembles the CO2 produced by a medium-sized coal-fired power plant. Since a heat balance is solved, a heat influx at the center boundary also has to be defined. We assume that CO2 is injected at reservoir conditions. Note that the reservoir conditions change with time and heat influx is adapted to these changing conditions. Top, bottom and side boundaries are closed to any flux (Neumann). At the lateral boundary, hydrostatic pressure, geothermal temperature and fully water saturated conditions are assumed. This resembles a laterally infinite aquifer with impermeable cap- and baserock. To conduct the computationally very expensive simulations, the CO2 Community Grid was used [32]. The CO2 Community Grid is a virtual research environment (VRE) for massively parallel simulations related to CO2 storage in porous media. The VRE allows for the usage of a number of supercomputers located in Norway. Table 3 lists some details on the supercomputers used. The grid infrastructure is controlled by a kind of meta-computing system, based on the Nordugrid [50] software ARC (Advanced resource connector [17]), a middleware for lightweight computational grids, which provides a unified access to supercomputers. The gateway to the VRE is through a service computer. This computer gives access to users and provides them a simple, unified access to supercomputer platforms. The system allows the execution of multiple parallel simulations simultaneously. During development, testing, and production of the results, statistics

873

Table 3 Details on the supercomputers used in this study. Frontend name

Location

ce02.titan.uio.no

Univ. of Oslo

No. of No. of CPU nodes cores type 304

2432

norgrid.uit.no Univ. of Tromso 704 hexgrid.bccs.uib.no Univ. of Bergen 1388

5632 5552

Theoretical total peak [Tflops]

Xeon/ 21.5 Opteron Xeon 2 60 Opteron 51

on the number of CPUs used per job have been collected. The maximum number of CPUs used for an individual job was 512 while most jobs used 64 CPUs. 5. Results Typical results are shown in Fig. 6 for two cases. In both cases, the amount of injected CO2 is 20 MtCO2. Randomly sampled primary parameters for Case 11 are φ11 = 0.19, dTdz11 = 0.020364 K/m, D11 = 1920.24 m, AnIso11 = 0.00898, and for Case 48 are φ48 = 0.26, dTdz48 = 0.053162 K/m, D48 = 640.08 m, AnIso11 = 0.02255. These settings lead to CO2 density around 790 kg/m3 in Case 11 and 330 kg/m3 in Case 48. In Case 11 no leakage occurs up to the time represented by the figure, whereas for Case 48 leakage occurs. Given the definition of risk in Eq. (1), the two terms of the equation are evaluated separately and then combined to yield risk. In Fig. 7 the likelihood of failure (LOF) is shown. Each point on the surface shows the combined results of all simulations conducted. The range of LOF is between zero (no case has failed) and one (all cases have failed). Initially no CO2 is injected and consequently no failure has occurred. After 7000 days all cases below rleakywell = 1100 m have failed and LOF is equal to one. For distances larger than rleakywell = 1100 m LOF reduces, since less cases fail within creasing distances. For rleakywell = 2000 m only 20% of cases fail after 7000 days of CO2 injection. Generally LOF increases with time and decreases with distance rleakywell. In Fig. 8 the consequence of failure (COF) is shown. As before, each point on the surface shows the combined results of all simulations conducted. The range of COF is between zero (no leakage occurred in any case) and 2.72·107 kg, where all cases produced considerable damage after 7000 days at rleakywell = 100 m. COF increases with time and decreases with distance rleakywell. For any given distance rleakywell COF increases quickly once damage has started to occur. For example, the COF is larger than zero for a leaky well at 200 m distance after 91 days of injection, then increases to 106 kg at 792 days, and to 107 kg after 5401 days. Another example is a leaky well in 2000 m distance, where the COF is larger than zero after 2409 days of injection, increases to 104 kg at 2806 days, and to 105 kg after 4819 days. Risk is the product of the likelihood of failure and the consequence of failure and is shown in Fig. 9. The range of risk is identical to the range of COF because LOF has a value between zero and one. Risk increases with time and decreases with distance rleakywell, as one would expect. Risk contour lines (e.g., the red line in Fig. 9 identifies risk equal to 0.001 kg) are fitted by the power-law

A⋅logðrÞ + B

t½d = e

;

ð7Þ

where r is leaky well distance in metres, and A and B are power-fitted coefficients given in Table 4. The accuracy of fit of the power-function is shown in the last column of the table.

874

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

Fig. 6. Slice of the radially symmetric domain showing CO2 saturation after 20 years of continuous injection for Case 11 (left) and Case 48 (right). Shown is only the center (inner) part of the domain and the vertical axis is exaggerated. No leakage can occur in Case 11, whereas for Case 48, CO2 has spread beyond the leaky well distance and the potential for leakage is high.

where r is leaky well distance in metres, and Ai are coefficients given in Table 5. In Fig. 10 risk is shown for a few selected leaky well distances and for various primary parameters. In general, risk for all cases increases with time of injection and with smaller leaky well radii. For a leaky well with a radius of 200 m, risk is basically identical in all cases. This result is expected given the short time and spatial scale, because the results are independent of the multiphase flow and transport regime in the reservoir and consequently independent of reservoir parameters. The effect and sensitivity of primary parameters is discussed below: Fig. 10 (top-left) shows risk as a function of porosity. No clear trend of risk as a function of porosity can be seen. Cases with a porosity of 0.106 produce almost the same risk as those with a porosity of 0.19, and 0.32 respectively. With increasing porosity one would expect a retardation of the plume evolution velocity, simply

due to increased pore space, and consequently a later increase of risk, i.e. lower risk at any time. But the permeability also increases with porosity, since it is coupled by a Karman-Kozeny functionality (c.f. Section 3.2). The increased permeability leads to stronger gravity segregation and hence to increased plume evolution velocity below the caprock. Note that the vertical permeability is coupled to horizontal permeability via the Anisotropy parameter. The effects of the retarded plume evolution velocity due to increased porosity and stronger gravity segregation due to increased horizontal and vertical permeability cancel each other out, i.e. risk appears to be independent of porosity. Small differences in risk are due to the other (random) primary parameters (e.g. depth). Fig. 10 (top-right) shows risk for a given geothermal gradient. It can be seen that risk increases earlier with increasing geothermal gradient. This is because of the lower CO2 density for a higher reservoir temperature. Gravity segregation is stronger at higher temperatures and fast lateral plume evolution occurs below the caprock. For rleakywell = 1000 m, risk increases after 884 days for a thermal gradient of 0.0532 K/m, after 2287 days for a thermal gradient of 0.0323 K/m, and after 3538 days for a thermal gradient of 0.0204 K/m. Risk is also quite different after 7000 days, i.e. 4.31·105 kg, 8.96·105 kg, and 1.2·106 kg (Factor 3 in between).

Fig. 7. Likelihood of failure [–] calculated as nf/N. nf is the number of cases that failed for distance rleakywell [m] between injection well and leaky well at time t [d]), N [–] is the total number of cases simulated.

i Fig. 8. Consequence of failure [log kg] calculated as ∑N i = 1 n versus distance rleakywell [m] f between injection well and leaky well and time t [d]. Di [kg] is damage in case i, nf [–] is the number of cases that failed (produced damage), N [–] is the total number of cases simulated.

The time contour lines (c.f. black lines in Fig. 9) can be used to calculate risk dependent on leaky well radii r. They are fitted by an exponential function to the risk surface given by 9

i

Risk½log kg = ∑ Ai ⋅r ; i=0

ð8Þ

D

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

n

D

i Fig. 9. Risk [log kg] calculated as f ⋅∑N i = 1 n versus distance rleakywell [m] between N f injection well and leaky well and time t [d]. Di [kg] is damage in case i, nf [–] is the number of cases that failed (produced damage), N [–] is the total number of cases simulated.

Fig. 10 (bottom-left) shows risk as a function of depth below ground surface and shows that risk increases earlier with decreasing depth. This is also because of CO2 density: at lower density gravity segregation is stronger, leading to fast lateral plume evolution below the caprock. For rleakywell = 1000 m, risk increases after 884 days for a depth of 640 m, after 1433 days for a depth of 1403 m, and after 2684 days for a depth of 2943 m. However, after 7000 days, risk can be viewed as being identical for all depths. Fig. 10 (bottom-right) shows risk as a function of permeability anisotropy. Risk increases earlier with increasing anisotropy factor. The gravity segregation effect is stronger for larger values of the anisotropy factor, i.e. less anisotropic permeability. Parameter sensitivity is however lower, compared with the sensitivity to geothermal gradient and depth. For rleakywell = 1000 m, risk increases after 1952 days for an anisotropy of 0.04274, after 2165 days for an anisotropy of 0.01282, and after 2318 days for an anisotropy of 0.00898. This difference in time gets larger with increasing rleakywell. Results show that the risk of leakage is primarily dependent on the plume development and the breakthrough of the formed gravity tongue. The influence of buoyant flow is well known in reservoir engineering [11,44]. In a previous study [37,38] developed estimates for CO2 storage capacity by dimensional analysis. There the gravitational number Gr and the capillary number Ca are used to describe the influence of buoyant forces and capillary forces on the development of the CO2 plume. The gravitational number is defined as Gr =

ðϱw −ϱCO2 Þ⋅g⋅k ; μCO2 ⋅vcr

ð9Þ

Table 4 Power-fitted coefficients A and B in Eq. (7) to calculate risk contour lines. R2 indicates goodness of fit (that is the ratio of the sum of the squares of the regression to the total sum of the squares) between {0,1}, where a value closer to one indicates better fit. Risk contour line [log kg]

A

B

R2

0.001 1 2 3 4 5 6 7

1.5487 1.4654 1.4024 1.3565 1.4615 1.6048 1.3875 1.0798

−3.9691 −3.2771 −2.7152 −2.1370 −2.3706 −2.7717 −0.6535 2.8553

0.9953 0.9967 0.9980 0.9976 0.9887 0.9956 0.9991 0.9981

875

where vcr is a characteristic velocity in the system, for example, the average plume velocity. High values of Gr indicate a large influence of buoyant forces in comparison to viscous forces, thus a tendency towards gravity segregation. It can be shown that low Gr leads to high CO2 storage capacities, mainly because of the plume shape (less gravity segregation). Since high storage capacities should correlate with a low risk of leakage, risk should also be lower for smaller Gr. The gravitational number Gr can be related to the primary parameters considered in this risk study: porosity, depth, geothermal gradient, and an anisotropy factor kv/kh. Porosity does not appear in the definition of Gr. However, it is coupled to the permeability k. Low φ means low k, thus also a low Gr number. On the other hand, low porosity requires more matrix volume for a given amount of CO2 to be filled, thus a larger extent of the plume. As outlined above, the results of this study showed that the porosity does not strongly influence risk. The influence of depth on Gr is rather obvious. With depth the density difference ϱw − ϱCO2 decreases resulting in low Gr and in a lower risk. High geothermal gradients imply more risky scenarios leading to high density differences and high Gr numbers. Anisotropy is not considered in the definition of Gr but if kv is employed instead of the isotropic k in the definition of Gr, then low Gr is achieved by low kv which is in agreement with less tendency towards gravity segregation. Ref. [37] defined the capillary number as Ca =

k⋅pcr ; μCO2 ⋅vcr ⋅lcr

ð10Þ

where pcr and lcr represent characteristic values for the capillary pressure and a corresponding length scale in the reservoir. Low values of Ca indicate that viscous forces dominate over capillary effects, while high Ca numbers are favorable for good storage capacity and, presumably, also for low risk of leakage. High geothermal gradients lead to lower capillary pressures. If capillary pressure was the only parameter that depends on temperature, this would result in low Ca numbers. However, the viscosity and the characteristic velocities are also temperature dependent. Thus, Ca (and Gr) are not strongly dependent on temperature. Low porosities (and thereby low permeabilities) also result in low Ca numbers. In general it has been shown that low Gr and high Ca lead to low risk, but the results of this risk study lead to the conclusion that Gr dominates both the estimation of risk and storage capacity. Summarizing this discussion of the results, it can be concluded that high risk is produced by small leaky well distances, long injection time, high geothermal gradients, high permeability anisotropy, and by low depth. Risk, as it is defined here, is almost not dependent on porosity. 6. Qualitative sensitivity analysis The assumptions made for the secondary parameters, together with other simplifications, are discussed in the following with respect to their potential to lower (−) or increase (+) leakage (lead to later (−) or earlier (+) leakage): • Reservoir geometry +/−: An increased thickness of the reservoir, a negative dip angle towards the leaky well, and a sealing fracture in the reservoir in between the injection and the leaky well leads to lower/later leakage at the leaky well. If assumptions were vice versa, i.e. decreased thickness, positive dip angle and a sealing fracture on the opposite side of the injection well, this would lead to higher/ earlier leakage. Ref. [49] studied the influence of varying dipping angles and described the influence using dimensionless numbers. Ref. [27] suggests that an increasing slope of an aquifer accelerates residual trapping and that lateral migration of the injected CO2 traps the CO2 relatively quickly as residual saturation. The results were, however, obtained with simplified 2D models neglecting (amongst other simplifications) density and viscosity changes.

876

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

Table 5 Coefficients Ai fitted to Eq. (8) to calculate time contour lines. R2 indicates goodness of fit (that is the ratio of the sum of the squares of the regression to the total sum of the squares) between {0,1}, where a value closer to one indicates better fit. Time contour line [days]

A0

A1

A2

A3

A4

A5

A6

A7

A8

A9

R2

1000 2000 3000 4000 5000 6000 7000

1.0685E+00 5.0088E+00 7.4770E+00 8.4316E+00 8.0576E+00 8.1624E+00 8.3008E+00

1.4859E−01 4.9038E−02 −3.9704E−03 −1.9807E−02 −1.0807E−02 −1.0943E−02 −1.1898E−02

−1.5375E−03 −4.6774E−04 −7.6240E−06 9.8945E−05 3.9636E−05 3.8833E−05 4.2956E−05

8.01300E−06 2.0879E−06 7.2986E−08 −2.7483E−07 −8.7095E−08 −8.0753E−08 −8.8591E−08

−2.3983E−08 −5.2227E−09 −1.9675E−10 4.3676E−10 1.1024E−10 9.6220E−11 1.0433E−10

4.2862E−11 7.7607E−12 2.5317E−13 −4.2399E−13 −8.6559E−14 −7.0288E−14 −7.4669E−14

−4.5849E−14 −6.9974E−15 −1.7582E−16 2.5517E−16 4.2855E−17 3.1944E−17 3.2931E−17

2.8163E−17 3.7556E−18 6.6951E−20 −9.2601E−20 −1.2960E−20 −8.7698E−21 −8.7074E−21

−8.8019E−21 −1.1029E−21 −1.2949E−23 1.8511E−23 2.1768E−24 1.3262E−24 1.2621E−24

9.8655E−25 1.3634E−25 9.6473E−28 −1.5627E−27 −1.5516E−28 −8.4625E−29 −7.6946E−29

0.9997 0.9998 0.9995 0.9987 0.9973 0.9974 0.9964

• Diffuse leakage through caprock -: Diffusive leakage into shallower reservoirs could reduce leakage at the leaky well. The effect on risk depends of course on the judgment whether this diffusive leakage is acceptable (e.g. because leaking CO2 is trapped in the shallower aquifer) or is not acceptable (leads also to damage/risk). However,

Ref. [45] showed that if the sealing pressure of the caprock is not exceeded, leakage of CO2 by molecular diffusion is negligible during the short-term injection stage. If the sealing pressure is exceeded though, leakage rates can become (very) large. Ref. [24] shows that when the transport of chemicals primarily occurs by molecular

Fig. 10. Risk [log kg] versus time [d] for selected rleakywell distances (shown in the boxes) and various values of primary parameters including porosity (top-left), geothermal gradient (top-right), depth (bottom-left), and absolute permeability anisotropy (bottom-right).

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

•

•

•

•

•

•

diffusion in the liquid phase, CO2 leakage becomes self-limiting, i.e. pores become clogged after a very short time. In case of a fractured caprock, however, gaseous CO2 penetrates into the caprock and induces some enhancement in porosity and permeability, which reduces the sealing efficiency of the caprock. Permeability: +/−: A different functionality of permeability on porosity (e.g. [28,57], or a Karman-Kozeny model with modified parameters), could lead to reduction or increase of leakage. Considering the concept of risk developed here, leakage is lower/ occurs later for lower permeabilities, i.e. the CO2 plume is more compact [37]. Viscous fingering: +/−: The radially symmetric and homogeneous setup prevents viscous fingering from occurring. Since the injected CO2 is less viscous than the displaced brine, viscous fingering may lead to an earlier breakthrough than predicted when radial symmetry is assumed. Leakage events can become more likely to occur through viscous fingers. However, the resulting damage and the total amounts of leaking CO2 might be even less since capillary pressure and dissolution of CO2 into brine are acting against the protrusion of fingers. Thus, a higher probability for leakage events does not necessarily imply higher risk. In the literature [22,68] show that the distribution of permeability dominates the fluid displacement and development of fingers rather than hydrodynamic instability. Heterogeneity +/−: This study neglects reservoir heterogeneity even though it has dominant influence on potential finger development, as explained in the previous paragraph. For example, a high permeable channel structure towards the leaky well, leads to earlier leakage. On the other hand, low permeability structures in the reservoir, can reduce gravitational segregation, and consequently leakage is retarded. Ref. [20] studied the impact of different netsand-to-gross-shale ratios on the development of the CO2 plume shape and on dissolved, residually trapped, and mobile CO2 mass fractions. They conclude that with increasing amounts of shale in the reservoir, vertical movement of the plume is restricted and lateral movement encouraged. An increasing amount of shale and consequently a restricted vertical movement leads to reduced leakage. Consequently, for such cases an assumption of homogeneity is presumably a conservative approach on risk. In general, however, reservoir heterogeneity can lead to either increased/ earlier or decreased/later leakage leaving this point open to be addressed in more detail. Relative Permeability +/−: Considering a different shape of the relative permeability relation, or different irreducible saturations, could lead to either increased/earlier or decreased/later leakage. Ref. [3] measured relative permeability relations for CO2-brine systems in the Alberta basin in Canada, and these do influence the plume evolution behavior (as shown by Ref. [38]) in a similar way to the reservoir properties considered here (depth, geothermal gradient, etc.). Capillary Pressure +/−: The approach employed to scale a measured capillary pressure relation to actual reservoir conditions at a given pressure, temperature, and (constant) porosity is quite sophisticated. However, such a scaling assumes the same rock-type in all reservoirs. As this is not generally the case, a different capillary pressure relation might lead to an increase or decrease of leakage. Generally, stronger capillary forces lead to less leakage, since more CO2 is stored by capillary trapping ([30,38]). In a heterogeneous reservoir, the effect of high capillary entry pressures of low permeable shale structures amplifies the effects discussed in bullet Heterogeneity. Hysteresis -: Hysteretic behavior in either the relative permeability relation or in capillary pressure influences plume evolution, and thus influences leakage. Various authors ([19,33,66], [13], [67]) investigated hysteretic behavior by experimental setups or by numerical investigations. Incorporating hysteresis of any kind will lead to decreased/later leakage, since additional CO2 is trapped.

877

• Anisotropy model +/−: A different anisotropy distribution (here a primary parameter derived by conceptual modeling) could lead to increased/earlier or decreased/later leakage. Generally, an increase in anisotropy leads to decreased/later leakage due to reduced gravity segregation. • Salinity +/−: As salinity influences brine density, viscosity and the solubility of CO2 in brine, a higher or lower salinity than assumed influences leakage. Generally, lower salinity leads to an increase in CO2 solubility in brine [14], hence leakage is decreased/occurs later. • Leakage simulation -: In this study, the leakage process itself through the leaky well is not modeled. It is assumed that all of the mass passing by the leakage point will flow through the leaky well out of the reservoir, which should be expected to be a worst case scenario. To justify this approach, a comparison is made with the model of Ref. [15] who investigated CO2 leakage rates through an abandoned well. The considered reservoir has an open leaky well of 0.30 m diameter at a 100 m distance to the injection well, which connects to a shallower aquifer. Results of the model of Ref. [15] are compared with the leakage rates predicted using the simplified approach of this study in Fig. 11 for the reservoir parameters and the setup presented by Ref. [15]. The time when significant leakage starts to occur is almost identical. The difference is due to slightly different boundary conditions and the different reservoir geometry (single radially symmetric reservoir in this study versus a full 3D simulation with two reservoirs, connected by the leaky well in the center in the study of Ref. [15]). Leakage rates start to increase at a slower rate in the study of Ref. [15], reach a peak leakage rate, and decrease on the long term. The gradual increase of the leakage rate is explained by up-coning of the brine into the well [53]. The higher peak rate and later decrease is explained by thermodynamic and hydraulic processes in the leaky well. These results show that despite neglecting the leakage process the predicted leakage agrees satisfactory with a model that describes the leakage process in detail. It is therefore reasonable to employ a simplified model without simulating leakage to gain the benefit of faster simulation and greater flexibility in the post-simulation evaluation of risk. • Additional assumptions on leaky well properties -: Here the leaky well is assumed to penetrate the entire formation thickness and is assumed to be completely open to the atmosphere. In practice this is not the case. For example partial penetration of the formation, well plugs, incorporation of detailed leakage pathways in the leaky well,

Fig. 11. Comparison of leakage rates versus time obtained by Ref. [15] and calculated by the methodology presented in this study. The well diameter is 0.30 m. Leakage rates are given as a fraction of the CO2 rate injected.

878

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879

knowledge about material behavior etc. will always reduce leakage rates. • Injection scheme +/−: In this work a continuous injection of CO2 is assumed. By a modified injection scheme, i.e. alternating injection of gas and fresh brine (WAG [61]), additional CO2 could be residually trapped [30]. Ref. [38] showed that the injection rate also has an considerable influence on the shape of the plume and therefore on the potential leakage. High injection rates lead to cylindrical plume shapes, which leads to less lateral extent of the plume. • Injection well design -: By intelligent design of the injection well (horizontal setup, screening in lower regions, etc) leakage could be reduced (could be deferred to later times), see also Ref. [40]. • Geochemistry -: Geochemical processes can lead to mineral trapping in the reservoir ([48,63,69]), which could lead to lower/ later leakage. These processes are important in the long term, i.e. later than 10-100 years, but are not considered in this study. If the caprock is subject to chemical processes, this can lead to very large leakage rates as stated in the Diffuse leakage through caprock bullet. 7. Conclusions The leakage risk concept developed utilizes a well tested simulation code to derive risk of potential CO2 leakage through preexisting leaky wells, subject to multiphase transport processes occurring in a reservoir. Leakage is not explicitly modeled, but calculated as the accumulated CO2 mass that has spread within the leaky well sector beyond the leaky well distance of interest at time t. The reservoir is assumed to have a simplified geometry and is homogeneous. The study takes into account four independent primary parameters randomly sampled from probability density functions derived from a large reservoir properties database. The risk surface derived is thus an average risk over the considered primary parameter ranges. Additional leaky wells in the surrounding of the injection well can be easily incorporated by summing up the individual risk for each well provided that leaky wells do not interact or significantly affect the CO2 distribution. Since other than primary reservoir properties are fixed or are dependent on primary parameters, different assumptions made for these parameters or for the reservoir could lead to lower/higher or earlier/later leakage at the leaky well, and hence to lower/larger risk. This is discussed in Section 6. This study adopts a conservative approach to estimating the parameters for determining the leakage of CO2 through abandoned wells. However, risk can be underestimated if reservoir geometry is very different (large dip, small height) or if the relative permeability shows high CO2 permeability together with large irreducible water saturation (significantly larger than 0.3). The main findings and conclusions are: • The first quantitative approach for the evaluation of risk with respect to CO2 leakage through pre-existing wells based on comprehensive reservoir properties statistics has been presented. • Within the given framework, a range of possible risks is defined. This can be used to determine whether an individual site is relatively good or not. Hence, assistance is given to experts when rating storage sites with unknown/uncertain reservoir properties. The cumulative risk for any site with given leaky well distances can be calculated, and sites can be compared to each other. Thus, the relative risk might assist in making decisions on where to conduct further investigations or to help experts when utilizing more comprehensive screening and ranking frameworks. • The study identifies reservoir parameters of importance to risk assessment. Among the four selected independent primary parameters, the depth of the reservoir and the geothermal gradient are shown to have the greatest influence on risk. Anisotropy is influential only for short distances from the injection well. In this

study, risk is independent of porosity (due to coupling of permeability to porosity). An ideal reservoir should thus be located at great depth, should have a low geothermal gradient, and should have a high anisotropy. • The placement of the injection well can be optimized with respect to risk arising from abandoned wells in the surroundings. For a given injection well location, the combined risk for any number of leaky wells in the surroundings can be calculated analytically for the time of interest. Thus, it is possible to compare risk for several injection well locations and pre-select the location yielding the lowest risk. Acknowledgements The authors gratefully acknowledge funding and support from the CO2SINK project (SES6-CT-2004-5025599) sponsored by the Commission of the European Communities and the industry (Statoil, RWE, Shell, Schlumberger, Vattenfall, and VNG), the GEOTECHNOLOGIEN program of the German Federal Ministry of Education and Research and the Deutsche Forschungsgemeinschaft, the Nordic Data Grid Facility (NDGF) and the Norwegian national infrastructure for computational science (NOTUR). References [1] Batzle M, Wang Z. Seismic properties of pore fluids. Geophysics 1992;57: 1396–408. [2] Bennion B, Bachu S. Relative Permeability Characteristics for Supercritical CO2 Displacing Water in a Variety of Potential Sequestration Zones in the Western Canada Sedimentary Basin. Society of Petroleum Engineers SPE 95547; 2005. [3] Bennion B, Bachu S. Drainage and Imbibition Relative Permeability Relationships for Supercritical CO2/Brine and H2S/Brine Systems in Intergranular Sandstone, Carbonate, Shale, and Anhydrite Rocks. SPE Reservoir Evaluation & Engineering 2008;11(3). [4] Benson S. Lessons learned from industrial and natural analogs for health, safety and environmental risk assessment for geologic storage of carbon dioxide. Carbon Dioxide Capture for Storage in Deep Geologic Formations - Results from the CO2 Capture Project, Elsevier Publishing, UK 2; 2005. [5] Benson S, Hepple R, Apps J, Tsang C, Lippmann M. Lessons Learned from Natural and Industrial Analogues for Storage of Carbon Dioxide in Deep Geological Formations. Energy Citations Database, Report Number LBNL–51170, Lawrence Berkeley National Laboratory, Berkeley, CA, U.S; 2005. [6] Bielinski, A., 2006. Numerical Simulation of CO2 Sequestration in Geological Formations. Ph.D. thesis, Universität Stuttgart, Mitteilungsheft Nr. 155, Institut für Wasserbau, Universität Stuttgart. [7] Brooks A, Corey A. Hydraulic Properties of Porous Media. Colorado State University Hydrology Paper; 1964. [8] Celia M, Bachu S, Nordbotten J, Gasda S, Dahle H. Quantitative Estimation of CO2 Leakage from Geological Storage: Analytical Models, Numerical Models, and Data Needs. In: Rubin ES, Keith DW, Gilboy CF, editors. 7th International Conference on Greenhouse Gas Control Technologies (GHGT-7). Vancouver, Canada, Sept. 5-9; 2004. p. 9. [9] Class H, Ebigbo A, Helmig R, Dahle H, Nordbotten J, Celia M, et al. A benchmarkstudy on problems related to CO2 storage in geologic formations — summary and discussion of the results. Computers and Geosciences 2009;13(4):409–34. [10] Class H, Helmig R, Bastian P. Numerical simulation of non-isothermal multiphase multicomponent processes in porous media. — 1. An efficient solution technique. Advances in Water Resources 2002;25:533–50. [11] Dake LP. Fundamentals of Reservoir Engineering. Elsevier Science B.V; 1978. [12] Daubert T, Danner R. Physical and thermodynamic properties of pure chemicals: data compilation. Design Institute for Physical Property Data; 1989. [13] Doughty C. Modeling geologic storage of carbon dioxide: comparison of nonhysteretic and hysteretic characteristic curves. Energy Conversion and Management 2007;48(6):1768–81. [14] Duan Z, Sun R. An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chemical Geology 2003;193:257–71. [15] Ebigbo A, Class H, Helmig R. CO2 leakage through an abandoned well: problemoriented benchmarks. Computers and Geosciences 2007;11(2):103–15. [16] Ennis-King J, Paterson L. The role of convective mixing in the long-term storage of carbon dioxide in deep saline formations. Society of Petroleum Engineers Journal 2005;10(3):349–56. [17] Ellert M, Gronager M, Konstantinov A, Konya B, Lindemann J, Livenson I, et al. Advanced resource connector middleware for lightweight computational grids. Future Generation Computer Systems 2007;23:2191–240. [18] Fenghour A, Wakeham W, Vesovic V. The viscosity of carbon dioxide. Journal of Physical and Chemical Reference Data 1998;27(1):31–44. [19] Flett M, Gurton R, Taggart I. The function of gas–water relative permeability hysteresis in the sequestration of carbon dioxide in saline formations. Society of Petroleum Engineers SPE 88485-MS; 2005.

A. Kopp et al. / Advances in Water Resources 33 (2010) 867–879 [20] Flett M, Gurton R, Weir G. Heterogeneous saline formations for carbon dioxide disposal: impact of varying heterogeneity on containment and trapping. Journal of Petroleum Science and Engineering 2007;57(1–2):106–18. [21] Garcia J. Density of Aqueous Solutions of CO2. Tech. rep., LBNL Report 49023, Lawrence Berkeley National Laboratory, Berkeley, CA, U.S.A; 2001. [22] Garcia J, Pruess K. Flow instabilities during injection of CO2 into saline aquifers. Proceedings TOUGH Symposium 2003, Lawrence Berkeley National Laboratory, Berkeley, CA, U.S.A; 2003. [23] Gasda S, Bachu S, Celia M. The potential for CO2 leakage from storage sites in geological media: analysis of well distribution in mature sedimentary basins. Environmental Geology 2004;46(67):707–20. [24] Gherardi F, Xu T, Pruess K. Numerical modeling of self-limiting and self-enhancing caprock alteration induced by CO2 storage in a depleted gas reservoir. Chemical Geology 2007;244:103–29. [25] Haszeldine RS, Quinn O, England G, Wilkinson M, Shipton ZK, Evans JP, et al. Natural geochemical analogues for carbon dioxide storage in deep geological porous reservoirs, a United Kingdom perspective. Oil & Gas Science and Technology 2005;60(1):33–49. [26] Helmig R, Class H, Huber R, Sheta H, Ewing R, Hinkelmann R, et al. Architecture of the Modular Program System MUFTE_UG for simulating multiphase flow and transport processes in heterogeneous porous media. Mathematische Geologie 1998;2:123–31. [27] Hesse M, Orr JF, Tchelepi H. Gravity currents with residual trapping. Journal of Fluid Mechanics 2008;611:35–60. [28] Holtz M. Residual Gas Saturation to Aquifer Influx: A Calculation Method for 3-D Computer Reservoir Model Construction. SPE 75503, paper presented at the SPE Gas Technology Symposium held in Calgary, Alberta, Canada, 30 April - 2 May 2002; 2002. [29] IAPWS. Release on the International Association for the Properties of Water and Steam Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam; 1997. (http://www.iapws.org/). [30] Ide S, Jessen K, F.M., O. J.. Storage of CO2 in saline aquifers: effects of gravity, viscous, and capillary forces on amount and timing of trapping. International Journal of Greenhouse Gas Control 2007;1:481–91. [31] IPCC. Carbon dioxide capture and storage. Intergovernmental Panel on Climate Change. uk: Cambridge university press; 2005. p. 431. (http://www.ipcc.ch/ ipccreports/special-reports.htm). [32] Johannsen K, Tourunen O, Kleist J. The CO2-Community Grid. http://wiki.ndgf.org/ display/ndgfwiki/CO22008. [33] Juanes R, Spiteri E, Orr F, Blunt M. Impact of relative permeability hysteresis on geological CO2 storage. Water Resources Research 2006;42. [34] Kaplan S. The words of risk analysis. Risk Analysis 1997;17(4):407–8. [35] Kopp, A., Bielinski, A., Class, H., Helmig, R., Hurter, S., Bech, N., 2006. Deliverable 6.1–1 Interim Progress Report- Concept for numerical modelling (Intercomparison Study). Progress report, not published, EU CO2SINK Project (Contract Number SES6-CT-2004–502599). [36] Kopp A, Class H, Helmig R. Sensitivity analysis of CO2 injection processes in brine aquifers. Proceedings European Geosciences Union (EGU) General Assembly, 18th April 2007, Vienna, Austria; 2007. [37] Kopp A, Class H, Helmig R. Investigations on CO2 storage capacity in saline aquifers — part 1: dimensional analysis of flow processes and reservoir characteristics. International Journal of Greenhouse Gas Control 2009;3:263–76. [38] Kopp A, Class H, Helmig R. Investigations on CO2 storage capacity in saline aquifers — part 2: estimation of storage capacity coefficients. International Journal of Greenhouse Gas Control 2009;3:277–87. [39] Kopp, A., 2009. Evaluation of CO2 injection processes in geological formations for site screening. Ph.D. thesis, Universität Stuttgart, Mitteilungsheft Nr. 182, Institut für Wasserbau, Universität Stuttgart. [40] Kovscek AR, Cakici CM. Geological storage of carbon dioxide and enhanced oil recovery II: cooptimization of storage and recovery. Energy Conversion and Management 2005;46:1941–56. [41] Kumar A, Ozah RC, Noh M, Pope GA, Bryant SL, Sepehrnoori K, et al. Reservoir simulation of CO2 storage in deep saline aquifers. Society of Petroleum Engineers Journal 2005;10(3):336–48. [42] Kvamme B, Kuznetsova T, Hebach A, Oberhof A, Lunde E. Measurements and modelling of interfacial tension for water + carbon dioxide systems at elevated pressures. Journal of Computational Materials Science 2007;38:506–13. [43] Laha IMRG, Chakravarti JR. Handbook of Methods of Applied Statistics, Volume I. John Wiley and Sons; 1967.

879

[44] Lake L. Enhanced Oil Recovery. Englewood Cliffs, New Jersey: Prentice-Hall, Inc.; 1989. [45] Li Z, Dong M, Li S, Huang S. CO2 sequestration in depleted oil and gas reservoirs caprock characterization and storage capacity. Energy and Conversion Management 2006;47:1372–82. [46] Maul P, Metcalf R, Pearce J, Savage D, West J. Performance assessment for the geologic storage of carbon dioxide: learning from the radioactive waste disposal experience. International Journal of Greenhouse Gas Control 2007;1:444–55. [47] Michaelides E. Thermodynamic properties of geothermal fluids. Geothermal Resources Council Transactions 1981;5:361–4. [48] Mito S, Xue Z, Ohsumi T. Case study of geochemical reactions at the Nagaoka CO2 injection site, Japan. International Journal of Greenhouse Gas Control 2008;2:309–18. [49] Mosthaf, K., 2007. CO2 Storage into Dipped Saline Aquifers Including Ambient Water Flow. Diploma Thesis, Institute of Hydraulic Engineering, University Stuttgart, Germany, http://www.iws.uni-stuttgart.de/publikationen/ausgabe. php?person=1&autor[0]=1195. [50] NDGF. Nordic Data Grid Facility. http://www.ndgf.org2008. [51] Nordbotten J, Celia M, Bachu S. Analytical solutions for leakage rates through abandoned wells. Water Resources Research 2004;40(4):W04204. [52] Nordbotten J, Celia M, Bachu S. Injection and storage of CO2 in deep saline aquifers: analytical solution for CO2 plume evolution during injection. Transport in Porous Media 2005;58(3):339–60. [53] Nordbotten J, Celia M, Bachu S, Dahle H. Semi-analytical solution for CO2 leakage through an abandoned well. Environmental Science & Technology 2005;39(2):602–11. [54] NPC. U.S. National Petroleum Council Public Database. http://www.netl.doe.gov. 1984. [55] Oldenburg C. Screening and ranking framework for geologic CO2 storage site selection on the basis of health, safety, and environmental risk. Environmental Geology 2007. [56] Pacala S. Global constraints on reservoir leakage. In: Gale J, Kaya Y, editors. Proceedings of the Sixth International Conference on Greenhouse Gas Control Technologies; 2003. p. 267–72. [57] Pape H, Clasuer C, Iffland J. Permeability prediction based on fractal pore-space geometry. Geophysics 1999;64(5):1447–60. [58] Plug W-J, Bruining J. Capillary pressure for the sand–CO2–water system under various pressure conditions. Application to CO2 sequestration. Advances in Water Resources 2007;30(11):2339–53. [59] Pruess K. On CO2 fluid flow and heat transfer behavior in the subsurface, following leakage from a geologic storage reservoir. Environmental Geology 2008;54:1677–86. [60] Pruess K, Bielinski A, Ennis-King J, Fabriol R, Le Gallo Y, Garcia J, et al. Code intercomparison builds confidence in numerical models for geologic disposal of CO2. In: Gale J, Kaya Y, editors. GHGT-6 Conference Proceedings: Greenhouse Gas Control Technologies; 2003. p. 463–70. [61] Qi R, LaForce TC, Blunt MJ. Design of carbon dioxide storage in aquifers. International Journal of Greenhouse Gas Control 2009;3:195–205. [62] Riaz A, Hesse M, Tchelepi H, Orr FM. Onset of convection in a graviatationally unstable, diffusive boundary layer in porous media. Journal of Fluid Mechanics 2006;548:87–111. [63] Rosenbauer R, Koksalan T, Palandri J. Experimental investigation of CO2–brine–rock interactions at elevated temperature and pressure: implications for CO2 sequestration in deep-saline aquifers. Fuel Processing Technology 2005;86:1581–97. [64] Scheidegger AE. The physics of flow through porous media. New York, NY: Macmillan Co.; 1960. [65] Span R, Wagner W. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa. Journal of Physical and Chemical Reference Data 1996;25(6):1509–96. [66] Spiteri E, Juanes R, Blunt M, Orr Jr F. Relative-Permeability Hysteresis: Trapping Models and Application to Geological CO2 Sequestration. SPE Annual Technical Conference and Exhibition, 9-12 October 2005, Dallas, U.S. SPE 96448-MS; 2005. [67] Spiteri E, Juanes R, Blunt M, Orr Jr F. A new model of trapping and relative permeability hysteresis for all wettability characteristics. Society of Petroleum Engineers 2008;13(3):277–88 (SPE 96448-PA). [68] Tchelepi HA, Orr FM. Interaction of viscous fingering, permeability heterogeneity, and gravity segregation in three dimensions. SPE Reservoir Engineering 1994;SPE 25235:266–71. [69] Xu T, Apps J, Pruess K, Yamamoto H. Numerical modeling of injection and mineral trapping of CO2 with H2S and SO2 in a sandstone formation. Chemical Geology 2007;242:319–46.