- Email: [email protected]
Toward quantitative CO2 storage estimates from time-lapse 3D seismic travel times: An example from the IEA GHG Weyburn–Midale CO2 monitoring and storage project
International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Contents lists available at SciVerse ScienceDirect
International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc
Toward quantitative CO2 storage estimates from time-lapse 3D seismic travel times: An example from the IEA GHG Weyburn–Midale CO2 monitoring and storage project D.J. White Geological Survey of Canada, Ottawa, Canada
a r t i c l e
i n f o
Article history: Received 12 October 2012 Received in revised form 29 December 2012 Accepted 31 January 2013 Available online 13 March 2013 Keywords: 4D Seismic Time-lapse CO2 Storage Monitoring
a b s t r a c t Time-lapse seismic interval travel time differences are determined during 7 years of CO2 injection within the Weyburn oil field, Saskatchewan, Canada. Travel time difference maps are used in conjunction with geophysical and geological logs, and depth dependant fluid properties to provide upper bounds on the amount of CO2 that may reside within various geological intervals. Calculated travel time sensitivity values increase from 0.01 ms/m to ∼0.25 ms/m moving upward from the reservoir depth (1450 m) to the shallowest aquifer considered (600 m). The seismic-based apportionment of CO2 shows the fraction of CO2 within or below the reservoir increasing from at least 62–70% after 1 year of injection to at least 92–94% after 7 years of injection, whereas the proportion of the total injected CO2 in the zone immediately overlying the reservoir decreases from a maximum of 30–36% to ∼5–6% over the same time period. The estimates of CO2 quantities residing within this zone are significantly overestimated due to clearly identified pressure effects in this zone. The maximum estimated proportion of CO2 residing in either interval above the regional sealing formation is ≤1%. The maximum amount of CO2 potentially residing above the regional seal by 2007 is <56, 000 tonnes, with a maximum of 0.31 Mtonnes just above the reservoir, and 4.6–4.7 Mtonnes in the reservoir-containing interval. Crown Copyright © 2013 Published by Elsevier B.V. All rights reserved.
1. Introduction The primary objective of time-lapse seismic monitoring in CO2 storage projects is to determine the subsurface distribution of CO2 . This includes mapping the location of CO2 within the primary injection zone (the reservoir), but also monitoring for CO2 that is “out-of-zone”. The latter may occur either as a result of CO2 migrating from the reservoir or alternatively as a result of out-of-zone injection. Time-lapse interval travel times determined from the seismic data provide a useful means of monitoring the overburden above the reservoir. Although they lack the resolving power of amplitude differences, they provide an integrated response of the entire interval. Seismic amplitudes (e.g., White, 2013), though very sensitive to small changes in impedance, are also more prone to noise. In this paper, P-wave interval travel time differences are determined for several depth intervals including the reservoir and overlying saline formations. Travel time difference maps are used in conjunction with geophysical and geological logs, and depth (pressure) dependant fluid properties to provide an upper bound on
E-mail address: [email protected]
the amount of CO2 that may reside within the geological column above the reservoir. White et al. (2011) used a crude version of this approach to produce an upper bound estimate. Here, this calculation is refined to account for the increased seismic sensitivity to CO2 residing at shallower depths and the higher compliance of clastic rocks in the overburden relative to the carbonate reservoir rocks. 2. Weyburn–Midale geology The Weyburn–Midale field is located within the Williston basin in southeastern Saskatchewan, Canada (see White, 2013 for location map). Since 2000, CO2 has been injected into the oil-bearing Midale beds as a means of tertiary oil recovery. The stratigraphic column within the Weyburn–Midale field has been characterized as a series of alternating aquifers and aquitards extending from the Precambrian basement to the surface (e.g., Fig. 1.3 of Whittaker et al., 2004). As shown in Fig. 1, lying above the Midale reservoir (M/V in Fig. 1) is a thick sequence of Mesozoic strata constituting primarily clastic rocks of which shale is the predominant lithology. These units include several regional saline aquifers (Shaunavon and Gravelbourg, which collectively are know as the Jurassic, Mannville and Newcastle formations) and aquitards (Watrous Formation, Vanguard Group, Joli Fou Formation, Colorado Group and Bearpaw Formation). The Watrous Formation is
1750-5836/$ – see front matter. Crown Copyright © 2013 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ijggc.2013.01.047
S96
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Fig. 1. Geophysical logs (gamma ray, density, and Vp ), formation tops and log-based zero-offset synthetic seismic traces. The intervals used for the travel time analysis are annotated on the right. Abbreviations: SSPK, 2nd White Specks; USH, Upper Shaunavon; GVBG, Gravelbourg; LGVBG, Lower Gravelbourg; WATR, Watrous; MISS, Mississippian; M, Midale Marly; V, Midale Vuggy.
considered to be the primary regional seal that will inhibit upward migration of CO2 from the reservoir. The Mannville Formation is a major regional aquifer. The aquifers provide the most likely locations for accumulation of upward migrating CO2 as they have significant porosity (mean values of 15–23%) and permeability, and are capped by low porosity, low permeability aquitard units (mostly shales). 3. Seismic travel time sensitivity to CO2 The goal of this study was to use seismic interval travel time changes to estimate the relative proportions of injected CO2 that reside within various strata within the subsurface. Ultimately, we would like to obtain quantitative estimates of the subsurface CO2 distribution. The ability to make these calculations is dependant on linking changes in the subsurface seismic properties with the presence of CO2 . The change in seismic properties associated with the injection of CO2 will be governed by several factors including host-rock lithology (e.g., clastic or carbonate), porosity, pressure and temperature, and pore fluid composition. A rock-fluid physics model is required to map changes in seismic properties associated with the presence of CO2 . For this purpose, we adopted commonly used formulae for determining the composite seismic properties (P-wave velocity (Vp ) and density) of fluid saturated rocks. In the absence of a calibrated Vp versus CO2 saturation relation (e.g., from in situ calibrations or flow-simulation constrained models), use of a more sophisticated rock physics model is unwarranted. Wood’s (1955) formula was used to determine the acoustic properties of mixed pore fluids and Gassmann’s (1951) relation was used to calculate the bulk seismic properties for a fluid-saturated rock. The properties of the saline fluids were determined using equations from Batzle and Wang (1992) and the CO2 properties using the empirical formulae from Han et al. (2010). Details of these formulae and how they are utilized can be found in Ramirez et al. (in
press), Smith et al. (2003) or Mavko et al. (2003). The effects of rock lithology and porosity are accounted for in these calculations. Pore pressure and temperature effects were included for the constituent pore fluids, but not for the rock matrix. The rock matrix properties used for each geological interval are shown in Table 1. Mean porosities for the aquifers of interest, as determined from regional well log compilations for the area, are shown there along with the layer depths and thicknesses taken from the interpreted well log shown in Fig. 1. Table 2 lists the fluid properties used for each of the aquifers. Fig. 2 illustrates the calculated depth-dependence of properties for pore fluids and a generic rock column with variable CO2 saturation. The Vp values of brine (Fig. 2a) have values that exceed those of CO2 by factors of 3–5, whereas for density, the brine values are ∼1.5 times greater. Below the vapour-to-supercritical fluid transition depth (∼800 m depth), the difference in Vp of supercritical CO2 relative to brine decreases with depth and thus the travel time sensitivity to the presence of CO2 will generally decrease with depth. Fig. 2b illustrates the effect of substituting CO2 for brine at various saturations for a sandstone. A similar depth-dependence of Vp is observed as for the fluids. The Wood–Gassmann formulation results in a very non-linear relationship between CO2 saturation and the change in seismic properties as compared to other possible substitution models. As illustrated in Fig. 3 for the reservoir units, the seismic P-wave velocity is much more sensitive to changes in CO2 saturation for low saturation levels as compared to high saturations. This will be true for the overlying saline formations as well. Practically, this means that in situations where high CO2 saturations occur, there will be ambiguity in estimating resident CO2 quantities from the seismic travel times. In this study, for the saline intervals above the reservoir it was assumed that the pores were originally 100% brine-filled and that the brine was completely replaced by CO2 in the zones where fluid substitution is considered. This results in an estimate of the maximum amount of CO2 resident within a given depth interval for an observed seismic change. Hence, the CO2 mass estimates for these intervals obtained from the observed seismic data represent maximum values. Seismic-based estimation of CO2 quantities in the oil reservoir is complicated by the additional presence of oil, and the significantly different properties of the lithologic units (Marly and Vuggy) that constitute the reservoir. Pre-CO2 injection oil/brine ratios for the Midale Marly and Vuggy units are approximately 50/50 and 30/70, respectively, as estimated from production history in the field. The difference between the two units is the result of more efficient production from the Vuggy unit during the long history of water-flooding in the field. Fig. 3 illustrates the changes in Vp that are expected for various modes of fluid substitution within these two units. The following can be observed. (1) The change in Vp versus CO2 saturation is most acute for low CO2 saturations. (2) The presence of oil in either unit significantly decreases the change in Vp due to introduction of CO2 in the pore space. (3) The fractional decrease in Vp for a given CO2 saturation is larger by a factor of ∼2 for the Marly as opposed to the Vuggy unit. (4) There is little difference in the dependence of Vp on CO2 saturation for the scenarios where CO2 replaces brine and oil in equal proportions, or CO2 preferentially replaces oil. Although, in the latter case the maximum attainable CO2 saturation is limited to that of the pre-injection oil saturation. To account for the ambiguities in these Vp -CO2 relationships, in this study CO2 quantities were calculated assuming that all of the reservoir CO2 resides in either the Marly or the Vuggy unit, and the oil is completely replaced by CO2 (i.e., CO2 saturations reach 50% and 30% in the Marly and Vuggy units in zones where CO2 substitution occurs). These scenarios provide a lower and upper estimate on the amount of CO2 residing in the reservoir based on reasonable assumptions about the mode of fluid
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
S97
Table 1 Geological properties for Mesozoic aquifer units as determined from regional compilations of log data. The salinity values were arbitrarily chosen at the low end of the range of typical values (10–50 g/L; Whittaker et al., 2004) for Mesozoic aquifers in this area. Abbreviations: ss, sandstone; sh, shale; l, limestone; d, dolomite. Geological unit
Rock type
Thickness (m)
Mean porosity (%)
Porosity RMS deviation (%)
Belly River Newcastle Mannville Shaunavon Gravelbourg Watrous Midale Marly/Vuggy
ss ss ss ss ss sh d/l
40 10 140 20 20 90 8/22
18 20 22.6 15 15 5 24/10
10 10 25 10 10 10 10/10
Salinity (ppm) 10,000 10,000 10,000 10,000 10,000 10,000 10,000
Matrix bulk and shear moduli (GPa)
Matrix bulk density (kg/m3 )
40, 44 40, 44 40, 44 40, 44 40, 44 40, 44 60, 30
2650 2650 2650 2650 2650 2650 2824/2740
Table 2 Fluid properties assigned to each geological unit (after Nader, 2012). Temperatures at depth were calculated using a constant vertical gradient of 0.025 ◦ C/m and an average surface temperature of 20 ◦ C. Pore pressures at depth were calculated by assuming a constant hydrostatic gradient of 0.0105 MPa/m. Geological unit
Depth (m)
Pore pressure (MPa)
Temperature (◦ C)
Saline water density (kg/m3 )
Saline water bulk modulus (GPa)
CO2 density (kg/m3 )
CO2 bulk modulus (GPa)
Belly River Newcastle Mannville Shaunavon/Gravelbourg Watrous Midale Marly Midale Vuggy
600 900 930 1175 1290 1400 1406
6.3 9.5 9.8 12.3 13.5 14.7 14.8
35 42 40 44 47 55 55
1000 1002 1003 1000 1000 1060 1060
2.354 2.427 2.423 2.470 2.487 2.900 2.900
267 407 473 540 568 590 590
0.024 0.031 0.050 0.064 0.074 0.083 0.083
substitution. However, in that higher CO2 saturations may be attained in zones where brine is partially displaced, they do not formally represent limiting estimates of CO2 quantities as was the case for the saline formations. Pressure effects that alter the seismic properties of the porous rock framework are not accounted for in either the saline formations or the reservoir units. It is very likely that pressure effects are only potentially significant at depths close to the reservoir. This clearly appears to be the case in the interval immediately above the reservoir as further discussed below. An example of the log-based fluid substitution results are shown for the Mannville aquifer in Fig. 4. As can be seen, substitution of brine with CO2 within this interval results in a reduction in Vp of between 12 and 17% which for the entire interval would produce a two-way vertical travel time delay of 12 ms. Results for the various aquifers and the Watrous Formation are summarized in Table 3. As can be seen in Fig. 5, the sensitivity of the seismic travel time delay increases greatly for each successively shallower aquifer unit
indicating that travel time monitoring should be very effective for detecting the presence of CO2 above the reservoir. 4. Interval travel time differences The stratigraphic column was divided into 4 depth intervals (Fig. 6) for the purposes of apportioning the injected CO2 within the subsurface: (1) surface to 2nd White Specks (SSPK), (2) SSPK to Lower Gravelbourg, (3) Watrous to Mississippian, and (4) Mississippian to Bakken. The 4th interval contains the reservoir where the CO2 is injected. The data are less reliable at shallower depths, and thus SSPK was the shallowest horizon considered. Interval 1 contains the Belly River aquifer, interval 2 contains the major aquifers above the reservoir (Newcastle, Mannville, Shaunavon and the upper part of the Gravelbourg), interval 3 constitutes the Watrous Formation, whereas interval 4 contains the reservoir and immediately underlying zones. These intervals were chosen as they are demarcated by strong coherent reflections from their
Fig. 2. (A) Density and Vp versus depth for brine and CO2 . Depth dependant pressure and temperature were used in determining the fluid/gas properties. (B) Vp versus CO2 saturation for a hypothetical sandstone with 20% porosity. Figures are modified from White (2012).
S98
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Fig. 3. (A) Vp and (B) fractional change (dVp ) in Vp versus CO2 saturation for different fluid substitution scenarios in the Midale Marly unit. (C) Vp and (D) fractional change (dVp ) in Vp versus CO2 saturation for different fluid substitution scenarios in the Midale Vuggy unit. The fluid substitution scenarios referred to in the legend are as follows. (1) Brine to CO2 : brine replaces CO2 . (2) Brine/oil to CO2 : CO2 replaces the brine/oil mixture in equal proportions. (3) Brine/oil to brine/CO2 : CO2 replaces only the oil in a brine/oil mixture.
boundaries that are well suited to the determination of maximum cross-correlation time shift analysis. Time-lapse interval travel time changes were determined using the 3D time-lapse seismic data from White (2013). Details of the data processing sequence applied to the data can be found there. The time-lapse data volumes utilized in this study had all post-stack calibration steps applied (as described by White, 2013) with the exception of the final step of time-variant cross-correlation time stretching. The latter step was eliminated in order to preserve the interval travel time changes that were analyzed here. The Lower Gravelbourg horizon was used as a reference horizon for determining relative time shifts during data processing. White (2013) notes the high repeatability achieved amongst the different vintages of the time-lapse data volumes with normalized-root-mean-square differences of 0.30–0.34 relative to the baseline data. Repeatability denotes the degree of similarity between a baseline seismic survey and subsequent monitor surveys except for differences induced by changes in the subsurface geology. High repeatability is achieved by reproducing survey parameters as closely as possible, including source-receiver geometry, source and receiver characteristics, data
processing and near-surface ground conditions. High repeatability is a fundamental requirement for robust time-lapse analysis. Interval travel time changes have been determined in the following manner. For each seismic horizon of interest, a time-shift was calculated that produced the maximum correlation value between the monitor and the baseline seismic volumes for a 20–40 ms window centered on the horizon. Then difference times were calculated for specific depth intervals by subtracting the maximum correlation time shifts (determined in the previous step) corresponding to the seismic horizons bounding the depth interval. Maximum correlation time shifts determined for the Lower Gravelbourg reference horizon had a mean value of 0.0 ms and root-mean-square deviations of 0.3–0.4 ms. In that this was the reference horizon used during processing, it is expected that this variance represents the best that can be achieved for the other horizons. Subsequently, interval differences exceeding 0.4–0.6 are expected to be significant. Based on this, a value of 0.5 ms was chosen as a threshold value for travel time differences that exceed the estimated uncertainty level. Values that had a magnitude less than a chosen threshold value were set to zero. In addition, a 3 × 3
Table 3 Summary of fluid substitution results by geological interval. For each interval the mean Vp determined from the log prior to CO2 injection is provided as well as the change in Vp , change in 2-way travel time (TWT) per meter, and the total time delay over each interval for replacement of brine to CO2 . Geological unit
Vp (m/s)
dVp (fraction)
dT/dz (ms/m)
Porosity (fraction)
Max. TWT delay (ms)
Belly River Viking/Newcastle Mannville Jurassic Watrous Midale Marly Midale Vuggy
2500 3000 3300 3300 3800 3700 5200
0.24 0.17 0.15 0.15 0.10 0.06 0.03
0.25 0.14 0.11 0.11 0.058 0.035 0.012
0.38 0.25 0.25 0.21 0.12 0.24 0.10
15.2 2.0 11.8 11.2 5.3 0.21 0.29
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
S99
Fig. 4. Pore fluid substitution modelling for the Mannville aquifer interval. The original log values (black) are measured in situ and correspond to brine-filled pores, whereas the modified properties (red) correspond to complete replacement of brine with CO2 . The effects of fluid substitution were calculated using the Wood–Gassmann formulation described in the text.
Fig. 7. Seismic interval travel time differences relative to the 2000 baseline data for (a) 2001, (b) 2002, (c) 2004 and (d) 2007, determined for stratigraphic interval 4 (Mississippian–Bakken) which contains the reservoir. The area shown is the Phase 1A CO2 flood area of the Weyburn field. Superposed are the locations of vertical wells (dots), horizontal CO2 injection wells (red lines) and horizontal production wells (black lines). The individual seismic bins (corresponding to pixels in the image) have dimensions of 40 m × 40 m.
Fig. 5. Log-based seismic travel time sensitivity for chosen geological units. Shown are the travel time changes per meter thickness of the formation when brine is completely replaced by CO2 in the pore space, except in the case of the Midale units where CO2 completely replaces oil. See the text for details.
median filter was applied to the maps to eliminate zones where non-zero travel time differences were of limited spatial extent (1–4 pixels). Although it is difficult to directly assess the validity of these assumptions, the good correlation of the travel time anomalies with reservoir simulations (e.g., White, 2013) and with well information as discussed below, provides supporting corroboration of the results. Interval travel time differences for the Mississippian–Bakken (reservoir containing) and the Watrous–Mississippian intervals are shown in Figs. 7 and 8, respectively. Fig. 9 shows the interval travel time differences for each depth interval for the 2004 versus the 2000 baseline data. As can be seen in Fig. 7, the interval which includes the reservoir has pronounced negative travel time anomalies of up to ∼2 ms that are clustered about the trajectories of the horizontal CO2 injection wells. In the immediately overlying interval (Fig. 7), there are a few smaller magnitude, but significant, travel time delays. In contrast (compare Fig. 9a and b with Fig. 9c and d), there are few (if any) significant travel time anomalies in the intervals above the regional sealing Watrous Formation. 5. CO2 quantity estimation
Fig. 6. Weyburn seismic data. Purple intervals represent aquitards. The reservoir zone is indicated in orange. Interval travel time differences have been determined for the depth intervals 1–4 as noted along the right of the plot with chosen formation tops. Note that the top of the seismic section corresponds to a travel time of 500 ms or approximately 750 m depth. Abbreviations: TWT, two-way travel time; SSPK, 2nd White Specks; LGVB, Lower Gravelbourg; WATR, Watrous; and MISS, Mississippian.
The cumulative travel time delays portrayed in Figs. 7–9 and their associated areal extent have been used as a proxy for the volume of CO2 present within the individual intervals. Whereas White et al. (2011) produced a crude estimate of CO2 proportions assuming equal significance of the anomalies regardless of the depth or geological interval, this study has used the log-based seismic
S100
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Fig. 8. Seismic interval travel time differences relative to the 2000 baseline data for (a) 2001, (b) 2002, (c) 2004 and (d) 2007, determined for stratigraphic interval 3 (Watrous–Mississippian) which resides immediately above the reservoir. The labels A and B identify vertical water injection wells that are discussed in the text. The area shown is the same as in Fig. 7. Superposed are the locations of vertical wells (dots), horizontal CO2 injection wells (red lines) and horizontal production wells (black lines). The individual seismic bins (corresponding to pixels in the image) have dimensions of 40 m × 40 m.
Fig. 9. Seismic interval travel time differences relative to the 2000 baseline data for 2004 for stratigraphic intervals 1–4. The area shown is the same as in Fig. 7. Superposed are the locations of vertical wells (dots), horizontal CO2 injection wells (red lines) and horizontal production wells (black lines). The individual seismic bins (corresponding to pixels in the image) have dimensions of 40 m × 40 m.
sensitivities (dT/dz) determined by fluid substitution for the various geological intervals (see Fig. 5). For an observed interval delay time T, the corresponding thickness of the zone where CO2 substitution has occurred is zsat = T/(dT/dz), from which an estimate of the CO2 mass MCO2 is given by MCO2 = A · zsat · ϕ · SCO2 · CO2 where A is the area over which the travel time delay is observed, ϕ is the fractional porosity of the geological interval, SCO2 is the CO2 saturation within the zone, and CO2 is the CO2 density determined for the appropriate pressure and temperature conditions. The resultant proportions of CO2 determined in this manner are shown in Fig. 9. These can be used along with the known CO2 net injection quantities (i.e., total CO2 injected minus CO2 recycled from oil production) to estimate the CO2 quantities within the various stratigraphic intervals as shown in Fig. 10. The seismic-based apportionment of CO2 (Fig. 9) shows 62–70% of the net injected CO2 residing within the reservoir-containing interval after one year of injection, increasing to 92–94% after 7 years. The estimated proportion of CO2 in the zone that immediately overlies the reservoir decreases from 30–36% after 1 year of injection to 5–6% after 7 years of injection. As discussed below, the CO2 estimates for this zone are likely significantly inflated by the contribution of pressure effects in the seismic travel time maps. The maximum estimated proportion of CO2 residing in either interval above the regional seal is ≤1%. Applying these relative amounts to
Fig. 10. Seismic-based relative mass estimates by stratigraphic interval. The range of values indicated for intervals 3 and 4 correspond to estimates for the reservoir interval in which it is assumed that CO2 replaces the oil in either the Marly or Vuggy units. As noted in the text, pressure effects substantially inflate the CO2 estimates for interval 3 which lies immediately above the reservoir.
the known net injected quantitites of CO2 (Fig. 10) indicates that the maximum amount of CO2 that could reside above the regional seal by 2007 is <56,000 tonnes. 4.6–4.7 Mtonnes are attributed to the reservoir and below, and a maximum of 0.31 Mtonnes potentially resides just above the reservoir.
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
S101
6. Discussion The interval travel time differences provide a systematic means of apportioning quantities of CO2 to the various geological intervals. However, there are limitations that must be appreciated. These limitations include uncertainties in the rock-physics model used to calculate the seismic sensitivities to CO2 , the difficulty of assessing very low levels of CO2 from seismic data that is inherently nonrepeatable, the uncertainties of how CO2 is distributed within a particular geological interval, the subjectivity in choosing a threshold level for the data, characterization of the subsurface using logs from a limited number of wells, and the assumption that all travel time delays are related to pore fluid substitution (as opposed to pressure changes, for example). As noted, a threshold was applied to the interval travel time differences in an attempt to eliminate small travel time variations that result from “noise” in the data. As a result, bona fide CO2 -related travel time changes that are less than 0.5 ms may be missed in the analysis. This will correspond to CO2 that is unaccounted for by the seismic data. The estimation of CO2 quantities demonstrated in this study provides a practical means of establishing an inventory of stored CO2 subject to the limitations listed above. For example, the ability to state that 1.0% is a maximum bound for CO2 potentially residing above the regional seal is useful for this purpose. However, in regard to detecting whether CO2 is leaking from the reservoir, this result is of limited value as it does not completely rule out the possibility of leakage. This is a consequence of the time-lapse seismic data being inherently non-repeatable and having limited sensitivity. An implicit assumption of the calculation is that all of the travel time delays are ascribed to fluid substitution effects. Decreases in effective pressure associated with increased pore pressures can have a similar effect. However, there are two mitigating factors here. First, significant pore pressure increases (e.g., 5-10 MPa at the reservoir level) are required to produce significant travel time delays. In most injection scenarios, such pressure changes would be restricted to the immediate vicinity of the injection reservoir, and thus would be expected to be generally insignificant at shallower levels. Thus, pressure contributions to the travel time delays should be most significant at the reservoir level and there they are generally expected to be smaller than fluid substitution effects (e.g., White, 2013). In the Weyburn example shown here, an exception to this almost certainly occurs in the zone immediately overlying the reservoir (see Fig. 8). There, several pronounced local travel time delays are observed which are centered on vertical wells that have been used exclusively for water injection since 2000, indicating that the observed seismic anomalies are non-CO2 related. For example, consider the two vertical water injectors A and B in Fig. 8. Water injection in well B occurred at a relatively high rate (∼1.26 × 105 m3 /year) from 2000 to 2002, was then shut-in until the middle of 2004, after which it continued at a reduced rate (∼0.5 × 105 m3 /year). The related travel time anomalies in Fig. 8 increase in intensity in 2001 and 2002 during the high injection rate period, and then fade in the 2004 and 2007 maps which correspond to the period of no water injection or reduced injection. Similarly, in well A, relatively high water injection rates to 2001 were followed by no injection till 2003 when injection began at a reduced level. In this case, the most pronounced travel time anomaly occurs in the 2004 map after 2 years of continuous water injection. The prominence of these non-CO2 related travel time delays artificially increases the estimate of CO2 for this layer. Thus, the recognition that the CO2 estimates presented here are maximum bounds is especially true for geological interval 3. The inferred pressure effect (described above) assigned to the interval immediately above the reservoir explains the apparent decrease in CO2 within this zone (interval 3) after 2002 that is
Fig. 11. Mass estimates calculated by combining the relative mass estimates from Fig. 10 with known net injection quantities. The range of values indicated for intervals 3 and 4 correspond estimates for the reservoir interval in which it is assumed that CO2 replaces the oil in either the Marly or Vuggy units.
depicted in Fig. 11. This apparent decrease of CO2 might be interpreted as due to upward migration into the overlying intervals (1 and 2). However, it most likely represents depressurization of this zone over time. Furthermore, the interval travel time delays for intervals 1 and 2 show no evidence for upward migrating CO2 . 7. Conclusions Time-lapse seismic interval travel time changes have been used to estimate the proportions of injected CO2 within four stratigraphic intervals including the reservoir interval and 3 overburden intervals. Log-based fluid substitution modeling formed the basis for calculating the seismic sensitivity to the presence of CO2 within the various intervals, accounting for the effects of rock lithology, porosity, pore pressure and temperature. Travel time sensitivity values increase by an order of magnitude (0.012 ms/m to ∼0.25 ms/m) moving upward from the reservoir depth (1450 m) to the shallowest aquifer considered (600 m), indicating that travel time monitoring should be very effective for detecting the presence of CO2 above the reservoir. The seismic-based relative apportionment of CO2 shows the fraction of net injected CO2 within the reservoir and below increasing from at least 62–70% after 1 year of injection to 92–94% after 7 years of injection, whereas the proportion in the immediately overlying zone decreases from 30–36% to ∼5–6% over the same time period. The estimate for CO2 residing in the interval just above the reservoir is significantly overestimated due to the contribution of pressure effects associated with water injection. The maximum estimated proportion of CO2 residing in either interval above the regional seal is ≤1%. The maximum amount of CO2 potentially residing above the regional seal by 2007 is <56,000 tonnes, with a maximum of 0.31 Mtonnes just above the reservoir, and 4.6-4.7 Mtonnes in the reservoir and below. In the Weyburn example presented here, pressure effects are considered of secondary importance in contributing to travel time delays at the reservoir level where CO2 quantities are expected to be large. However, pressure effects are very significant in the zone just above the reservoir and greatly affect the corresponding estimation of CO2 quantities for that zone. Thus, improved seismic-based CO2 quantity estimates will require the inclusion of field-based knowledge, incorporation of constraints from flow simulations, and adoption of rock-fluid physics models with site-specific calibration.
S102
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Acknowledgements Cenovus provided the wells logs and the seismic data for the 2004 and 2007 surveys. Initial seismic processing was done by Dave Secord at Sensor Geophysical. Discussions with Dave Cooper provided much insight based on his extensive experience in the Weyburn field operations. Barbara Dietiker contributed to the time-lapse seismic analysis and produced the log-based property compilations. Hampson-Russell Pro4D software was used for some of the time-lapse seismic analysis and fluid substitution modelling. This is contribution number 20120441 of the Geological Survey of Canada. References Batzle, M., Wang, Z., 1992. Seismic properties of pore fluids. Geophysics 57, 1396–1408. Gassmann, F., 1951. Uber die Elastizitat poroser Medien. Veirteljahrsschrift der Naturforschenden Gesellschaft in Zurich 96, 1–23. Han, D., Sun, M., Batzle, M., 2010. CO2 velocity measurements and models for temperatures up to 200 ◦ C and pressures up to 100 MPa. Geophysics 75, E123–E129. Mavko, G., Mukerji, T., Dvorkin, J., 2003. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media. Cambridge University Press, Cambridge.
Nader, R., 2012. Effectiveness of seismic detection modelling for CO2 migation scenarios near the Aquistore CO2 injection site, Estevan, Saskatchewan. B.Sc. Thesis. Carleton University, Ottawa. Ramirez, A., White, D., Hao, Y., Dyer, K., Johnson, J. Estimating reservoir permeabilities using the seismic response to CO2 injection and stochastic inversion. International Journal of Greenhouse Gas Control Technologies, http://dx.doi.org/10.1016/j.ijggc.2012.11.031, in press. Smith, T.M., Sondergeld, C.H., Rai, C.S., 2003. Gassmann fluid substitutions: a tutorial. Geophysics 68 (2), 430–440. White, D., 2012. Geophysical monitoring. In: Hitchon, B. (Ed.), Best Practices for Validating CO2 Geological Storage. Geoscience Publishing, Sherwood Park, pp. 155–210. White, D., 2013. Seismic characterization and time-lapse imaging during 7 years of CO2 flood in the Weyburn Field, Saskatchewan, Canada. International Journal of Greenhouse Gas Control, http://dx.doi.org/10.1016/j.ijggc.2013.02.006, in press. White, D.J., Meadows, M., Cole, S., Ramirez, A., Hao, Y., Carle, S., Duxbury, A., Samson, C., Kendall, J.-M., Verdon, J.P., Dietiker, B., Johnson, J., Morozov, I., 2011. Geophysical monitoring of the Weyburn CO2 flood: results during 10 years of injection. Energy Procedia 4, 3628–3635. Whittaker, S., et al., 2004. Theme 1: geological characterization. In: Wilson, M., Monea, M. (Eds.), IEA GHG Weyburn CO2 Monitoring and Storage Project Summary Report 2000–2004. From the Proceedings of the 7th International Conference on Greenhouse Gas Control Technologies. Petroleum Technology Research Centre, Regina, pp. 1–72 (Chapter 1). Wood, A.W., 1955. A Textbook of Sound. McMillan Co., New York.
Contents lists available at SciVerse ScienceDirect
International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc
Toward quantitative CO2 storage estimates from time-lapse 3D seismic travel times: An example from the IEA GHG Weyburn–Midale CO2 monitoring and storage project D.J. White Geological Survey of Canada, Ottawa, Canada
a r t i c l e
i n f o
Article history: Received 12 October 2012 Received in revised form 29 December 2012 Accepted 31 January 2013 Available online 13 March 2013 Keywords: 4D Seismic Time-lapse CO2 Storage Monitoring
a b s t r a c t Time-lapse seismic interval travel time differences are determined during 7 years of CO2 injection within the Weyburn oil field, Saskatchewan, Canada. Travel time difference maps are used in conjunction with geophysical and geological logs, and depth dependant fluid properties to provide upper bounds on the amount of CO2 that may reside within various geological intervals. Calculated travel time sensitivity values increase from 0.01 ms/m to ∼0.25 ms/m moving upward from the reservoir depth (1450 m) to the shallowest aquifer considered (600 m). The seismic-based apportionment of CO2 shows the fraction of CO2 within or below the reservoir increasing from at least 62–70% after 1 year of injection to at least 92–94% after 7 years of injection, whereas the proportion of the total injected CO2 in the zone immediately overlying the reservoir decreases from a maximum of 30–36% to ∼5–6% over the same time period. The estimates of CO2 quantities residing within this zone are significantly overestimated due to clearly identified pressure effects in this zone. The maximum estimated proportion of CO2 residing in either interval above the regional sealing formation is ≤1%. The maximum amount of CO2 potentially residing above the regional seal by 2007 is <56, 000 tonnes, with a maximum of 0.31 Mtonnes just above the reservoir, and 4.6–4.7 Mtonnes in the reservoir-containing interval. Crown Copyright © 2013 Published by Elsevier B.V. All rights reserved.
1. Introduction The primary objective of time-lapse seismic monitoring in CO2 storage projects is to determine the subsurface distribution of CO2 . This includes mapping the location of CO2 within the primary injection zone (the reservoir), but also monitoring for CO2 that is “out-of-zone”. The latter may occur either as a result of CO2 migrating from the reservoir or alternatively as a result of out-of-zone injection. Time-lapse interval travel times determined from the seismic data provide a useful means of monitoring the overburden above the reservoir. Although they lack the resolving power of amplitude differences, they provide an integrated response of the entire interval. Seismic amplitudes (e.g., White, 2013), though very sensitive to small changes in impedance, are also more prone to noise. In this paper, P-wave interval travel time differences are determined for several depth intervals including the reservoir and overlying saline formations. Travel time difference maps are used in conjunction with geophysical and geological logs, and depth (pressure) dependant fluid properties to provide an upper bound on
E-mail address: [email protected]
the amount of CO2 that may reside within the geological column above the reservoir. White et al. (2011) used a crude version of this approach to produce an upper bound estimate. Here, this calculation is refined to account for the increased seismic sensitivity to CO2 residing at shallower depths and the higher compliance of clastic rocks in the overburden relative to the carbonate reservoir rocks. 2. Weyburn–Midale geology The Weyburn–Midale field is located within the Williston basin in southeastern Saskatchewan, Canada (see White, 2013 for location map). Since 2000, CO2 has been injected into the oil-bearing Midale beds as a means of tertiary oil recovery. The stratigraphic column within the Weyburn–Midale field has been characterized as a series of alternating aquifers and aquitards extending from the Precambrian basement to the surface (e.g., Fig. 1.3 of Whittaker et al., 2004). As shown in Fig. 1, lying above the Midale reservoir (M/V in Fig. 1) is a thick sequence of Mesozoic strata constituting primarily clastic rocks of which shale is the predominant lithology. These units include several regional saline aquifers (Shaunavon and Gravelbourg, which collectively are know as the Jurassic, Mannville and Newcastle formations) and aquitards (Watrous Formation, Vanguard Group, Joli Fou Formation, Colorado Group and Bearpaw Formation). The Watrous Formation is
1750-5836/$ – see front matter. Crown Copyright © 2013 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ijggc.2013.01.047
S96
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Fig. 1. Geophysical logs (gamma ray, density, and Vp ), formation tops and log-based zero-offset synthetic seismic traces. The intervals used for the travel time analysis are annotated on the right. Abbreviations: SSPK, 2nd White Specks; USH, Upper Shaunavon; GVBG, Gravelbourg; LGVBG, Lower Gravelbourg; WATR, Watrous; MISS, Mississippian; M, Midale Marly; V, Midale Vuggy.
considered to be the primary regional seal that will inhibit upward migration of CO2 from the reservoir. The Mannville Formation is a major regional aquifer. The aquifers provide the most likely locations for accumulation of upward migrating CO2 as they have significant porosity (mean values of 15–23%) and permeability, and are capped by low porosity, low permeability aquitard units (mostly shales). 3. Seismic travel time sensitivity to CO2 The goal of this study was to use seismic interval travel time changes to estimate the relative proportions of injected CO2 that reside within various strata within the subsurface. Ultimately, we would like to obtain quantitative estimates of the subsurface CO2 distribution. The ability to make these calculations is dependant on linking changes in the subsurface seismic properties with the presence of CO2 . The change in seismic properties associated with the injection of CO2 will be governed by several factors including host-rock lithology (e.g., clastic or carbonate), porosity, pressure and temperature, and pore fluid composition. A rock-fluid physics model is required to map changes in seismic properties associated with the presence of CO2 . For this purpose, we adopted commonly used formulae for determining the composite seismic properties (P-wave velocity (Vp ) and density) of fluid saturated rocks. In the absence of a calibrated Vp versus CO2 saturation relation (e.g., from in situ calibrations or flow-simulation constrained models), use of a more sophisticated rock physics model is unwarranted. Wood’s (1955) formula was used to determine the acoustic properties of mixed pore fluids and Gassmann’s (1951) relation was used to calculate the bulk seismic properties for a fluid-saturated rock. The properties of the saline fluids were determined using equations from Batzle and Wang (1992) and the CO2 properties using the empirical formulae from Han et al. (2010). Details of these formulae and how they are utilized can be found in Ramirez et al. (in
press), Smith et al. (2003) or Mavko et al. (2003). The effects of rock lithology and porosity are accounted for in these calculations. Pore pressure and temperature effects were included for the constituent pore fluids, but not for the rock matrix. The rock matrix properties used for each geological interval are shown in Table 1. Mean porosities for the aquifers of interest, as determined from regional well log compilations for the area, are shown there along with the layer depths and thicknesses taken from the interpreted well log shown in Fig. 1. Table 2 lists the fluid properties used for each of the aquifers. Fig. 2 illustrates the calculated depth-dependence of properties for pore fluids and a generic rock column with variable CO2 saturation. The Vp values of brine (Fig. 2a) have values that exceed those of CO2 by factors of 3–5, whereas for density, the brine values are ∼1.5 times greater. Below the vapour-to-supercritical fluid transition depth (∼800 m depth), the difference in Vp of supercritical CO2 relative to brine decreases with depth and thus the travel time sensitivity to the presence of CO2 will generally decrease with depth. Fig. 2b illustrates the effect of substituting CO2 for brine at various saturations for a sandstone. A similar depth-dependence of Vp is observed as for the fluids. The Wood–Gassmann formulation results in a very non-linear relationship between CO2 saturation and the change in seismic properties as compared to other possible substitution models. As illustrated in Fig. 3 for the reservoir units, the seismic P-wave velocity is much more sensitive to changes in CO2 saturation for low saturation levels as compared to high saturations. This will be true for the overlying saline formations as well. Practically, this means that in situations where high CO2 saturations occur, there will be ambiguity in estimating resident CO2 quantities from the seismic travel times. In this study, for the saline intervals above the reservoir it was assumed that the pores were originally 100% brine-filled and that the brine was completely replaced by CO2 in the zones where fluid substitution is considered. This results in an estimate of the maximum amount of CO2 resident within a given depth interval for an observed seismic change. Hence, the CO2 mass estimates for these intervals obtained from the observed seismic data represent maximum values. Seismic-based estimation of CO2 quantities in the oil reservoir is complicated by the additional presence of oil, and the significantly different properties of the lithologic units (Marly and Vuggy) that constitute the reservoir. Pre-CO2 injection oil/brine ratios for the Midale Marly and Vuggy units are approximately 50/50 and 30/70, respectively, as estimated from production history in the field. The difference between the two units is the result of more efficient production from the Vuggy unit during the long history of water-flooding in the field. Fig. 3 illustrates the changes in Vp that are expected for various modes of fluid substitution within these two units. The following can be observed. (1) The change in Vp versus CO2 saturation is most acute for low CO2 saturations. (2) The presence of oil in either unit significantly decreases the change in Vp due to introduction of CO2 in the pore space. (3) The fractional decrease in Vp for a given CO2 saturation is larger by a factor of ∼2 for the Marly as opposed to the Vuggy unit. (4) There is little difference in the dependence of Vp on CO2 saturation for the scenarios where CO2 replaces brine and oil in equal proportions, or CO2 preferentially replaces oil. Although, in the latter case the maximum attainable CO2 saturation is limited to that of the pre-injection oil saturation. To account for the ambiguities in these Vp -CO2 relationships, in this study CO2 quantities were calculated assuming that all of the reservoir CO2 resides in either the Marly or the Vuggy unit, and the oil is completely replaced by CO2 (i.e., CO2 saturations reach 50% and 30% in the Marly and Vuggy units in zones where CO2 substitution occurs). These scenarios provide a lower and upper estimate on the amount of CO2 residing in the reservoir based on reasonable assumptions about the mode of fluid
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
S97
Table 1 Geological properties for Mesozoic aquifer units as determined from regional compilations of log data. The salinity values were arbitrarily chosen at the low end of the range of typical values (10–50 g/L; Whittaker et al., 2004) for Mesozoic aquifers in this area. Abbreviations: ss, sandstone; sh, shale; l, limestone; d, dolomite. Geological unit
Rock type
Thickness (m)
Mean porosity (%)
Porosity RMS deviation (%)
Belly River Newcastle Mannville Shaunavon Gravelbourg Watrous Midale Marly/Vuggy
ss ss ss ss ss sh d/l
40 10 140 20 20 90 8/22
18 20 22.6 15 15 5 24/10
10 10 25 10 10 10 10/10
Salinity (ppm) 10,000 10,000 10,000 10,000 10,000 10,000 10,000
Matrix bulk and shear moduli (GPa)
Matrix bulk density (kg/m3 )
40, 44 40, 44 40, 44 40, 44 40, 44 40, 44 60, 30
2650 2650 2650 2650 2650 2650 2824/2740
Table 2 Fluid properties assigned to each geological unit (after Nader, 2012). Temperatures at depth were calculated using a constant vertical gradient of 0.025 ◦ C/m and an average surface temperature of 20 ◦ C. Pore pressures at depth were calculated by assuming a constant hydrostatic gradient of 0.0105 MPa/m. Geological unit
Depth (m)
Pore pressure (MPa)
Temperature (◦ C)
Saline water density (kg/m3 )
Saline water bulk modulus (GPa)
CO2 density (kg/m3 )
CO2 bulk modulus (GPa)
Belly River Newcastle Mannville Shaunavon/Gravelbourg Watrous Midale Marly Midale Vuggy
600 900 930 1175 1290 1400 1406
6.3 9.5 9.8 12.3 13.5 14.7 14.8
35 42 40 44 47 55 55
1000 1002 1003 1000 1000 1060 1060
2.354 2.427 2.423 2.470 2.487 2.900 2.900
267 407 473 540 568 590 590
0.024 0.031 0.050 0.064 0.074 0.083 0.083
substitution. However, in that higher CO2 saturations may be attained in zones where brine is partially displaced, they do not formally represent limiting estimates of CO2 quantities as was the case for the saline formations. Pressure effects that alter the seismic properties of the porous rock framework are not accounted for in either the saline formations or the reservoir units. It is very likely that pressure effects are only potentially significant at depths close to the reservoir. This clearly appears to be the case in the interval immediately above the reservoir as further discussed below. An example of the log-based fluid substitution results are shown for the Mannville aquifer in Fig. 4. As can be seen, substitution of brine with CO2 within this interval results in a reduction in Vp of between 12 and 17% which for the entire interval would produce a two-way vertical travel time delay of 12 ms. Results for the various aquifers and the Watrous Formation are summarized in Table 3. As can be seen in Fig. 5, the sensitivity of the seismic travel time delay increases greatly for each successively shallower aquifer unit
indicating that travel time monitoring should be very effective for detecting the presence of CO2 above the reservoir. 4. Interval travel time differences The stratigraphic column was divided into 4 depth intervals (Fig. 6) for the purposes of apportioning the injected CO2 within the subsurface: (1) surface to 2nd White Specks (SSPK), (2) SSPK to Lower Gravelbourg, (3) Watrous to Mississippian, and (4) Mississippian to Bakken. The 4th interval contains the reservoir where the CO2 is injected. The data are less reliable at shallower depths, and thus SSPK was the shallowest horizon considered. Interval 1 contains the Belly River aquifer, interval 2 contains the major aquifers above the reservoir (Newcastle, Mannville, Shaunavon and the upper part of the Gravelbourg), interval 3 constitutes the Watrous Formation, whereas interval 4 contains the reservoir and immediately underlying zones. These intervals were chosen as they are demarcated by strong coherent reflections from their
Fig. 2. (A) Density and Vp versus depth for brine and CO2 . Depth dependant pressure and temperature were used in determining the fluid/gas properties. (B) Vp versus CO2 saturation for a hypothetical sandstone with 20% porosity. Figures are modified from White (2012).
S98
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Fig. 3. (A) Vp and (B) fractional change (dVp ) in Vp versus CO2 saturation for different fluid substitution scenarios in the Midale Marly unit. (C) Vp and (D) fractional change (dVp ) in Vp versus CO2 saturation for different fluid substitution scenarios in the Midale Vuggy unit. The fluid substitution scenarios referred to in the legend are as follows. (1) Brine to CO2 : brine replaces CO2 . (2) Brine/oil to CO2 : CO2 replaces the brine/oil mixture in equal proportions. (3) Brine/oil to brine/CO2 : CO2 replaces only the oil in a brine/oil mixture.
boundaries that are well suited to the determination of maximum cross-correlation time shift analysis. Time-lapse interval travel time changes were determined using the 3D time-lapse seismic data from White (2013). Details of the data processing sequence applied to the data can be found there. The time-lapse data volumes utilized in this study had all post-stack calibration steps applied (as described by White, 2013) with the exception of the final step of time-variant cross-correlation time stretching. The latter step was eliminated in order to preserve the interval travel time changes that were analyzed here. The Lower Gravelbourg horizon was used as a reference horizon for determining relative time shifts during data processing. White (2013) notes the high repeatability achieved amongst the different vintages of the time-lapse data volumes with normalized-root-mean-square differences of 0.30–0.34 relative to the baseline data. Repeatability denotes the degree of similarity between a baseline seismic survey and subsequent monitor surveys except for differences induced by changes in the subsurface geology. High repeatability is achieved by reproducing survey parameters as closely as possible, including source-receiver geometry, source and receiver characteristics, data
processing and near-surface ground conditions. High repeatability is a fundamental requirement for robust time-lapse analysis. Interval travel time changes have been determined in the following manner. For each seismic horizon of interest, a time-shift was calculated that produced the maximum correlation value between the monitor and the baseline seismic volumes for a 20–40 ms window centered on the horizon. Then difference times were calculated for specific depth intervals by subtracting the maximum correlation time shifts (determined in the previous step) corresponding to the seismic horizons bounding the depth interval. Maximum correlation time shifts determined for the Lower Gravelbourg reference horizon had a mean value of 0.0 ms and root-mean-square deviations of 0.3–0.4 ms. In that this was the reference horizon used during processing, it is expected that this variance represents the best that can be achieved for the other horizons. Subsequently, interval differences exceeding 0.4–0.6 are expected to be significant. Based on this, a value of 0.5 ms was chosen as a threshold value for travel time differences that exceed the estimated uncertainty level. Values that had a magnitude less than a chosen threshold value were set to zero. In addition, a 3 × 3
Table 3 Summary of fluid substitution results by geological interval. For each interval the mean Vp determined from the log prior to CO2 injection is provided as well as the change in Vp , change in 2-way travel time (TWT) per meter, and the total time delay over each interval for replacement of brine to CO2 . Geological unit
Vp (m/s)
dVp (fraction)
dT/dz (ms/m)
Porosity (fraction)
Max. TWT delay (ms)
Belly River Viking/Newcastle Mannville Jurassic Watrous Midale Marly Midale Vuggy
2500 3000 3300 3300 3800 3700 5200
0.24 0.17 0.15 0.15 0.10 0.06 0.03
0.25 0.14 0.11 0.11 0.058 0.035 0.012
0.38 0.25 0.25 0.21 0.12 0.24 0.10
15.2 2.0 11.8 11.2 5.3 0.21 0.29
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
S99
Fig. 4. Pore fluid substitution modelling for the Mannville aquifer interval. The original log values (black) are measured in situ and correspond to brine-filled pores, whereas the modified properties (red) correspond to complete replacement of brine with CO2 . The effects of fluid substitution were calculated using the Wood–Gassmann formulation described in the text.
Fig. 7. Seismic interval travel time differences relative to the 2000 baseline data for (a) 2001, (b) 2002, (c) 2004 and (d) 2007, determined for stratigraphic interval 4 (Mississippian–Bakken) which contains the reservoir. The area shown is the Phase 1A CO2 flood area of the Weyburn field. Superposed are the locations of vertical wells (dots), horizontal CO2 injection wells (red lines) and horizontal production wells (black lines). The individual seismic bins (corresponding to pixels in the image) have dimensions of 40 m × 40 m.
Fig. 5. Log-based seismic travel time sensitivity for chosen geological units. Shown are the travel time changes per meter thickness of the formation when brine is completely replaced by CO2 in the pore space, except in the case of the Midale units where CO2 completely replaces oil. See the text for details.
median filter was applied to the maps to eliminate zones where non-zero travel time differences were of limited spatial extent (1–4 pixels). Although it is difficult to directly assess the validity of these assumptions, the good correlation of the travel time anomalies with reservoir simulations (e.g., White, 2013) and with well information as discussed below, provides supporting corroboration of the results. Interval travel time differences for the Mississippian–Bakken (reservoir containing) and the Watrous–Mississippian intervals are shown in Figs. 7 and 8, respectively. Fig. 9 shows the interval travel time differences for each depth interval for the 2004 versus the 2000 baseline data. As can be seen in Fig. 7, the interval which includes the reservoir has pronounced negative travel time anomalies of up to ∼2 ms that are clustered about the trajectories of the horizontal CO2 injection wells. In the immediately overlying interval (Fig. 7), there are a few smaller magnitude, but significant, travel time delays. In contrast (compare Fig. 9a and b with Fig. 9c and d), there are few (if any) significant travel time anomalies in the intervals above the regional sealing Watrous Formation. 5. CO2 quantity estimation
Fig. 6. Weyburn seismic data. Purple intervals represent aquitards. The reservoir zone is indicated in orange. Interval travel time differences have been determined for the depth intervals 1–4 as noted along the right of the plot with chosen formation tops. Note that the top of the seismic section corresponds to a travel time of 500 ms or approximately 750 m depth. Abbreviations: TWT, two-way travel time; SSPK, 2nd White Specks; LGVB, Lower Gravelbourg; WATR, Watrous; and MISS, Mississippian.
The cumulative travel time delays portrayed in Figs. 7–9 and their associated areal extent have been used as a proxy for the volume of CO2 present within the individual intervals. Whereas White et al. (2011) produced a crude estimate of CO2 proportions assuming equal significance of the anomalies regardless of the depth or geological interval, this study has used the log-based seismic
S100
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Fig. 8. Seismic interval travel time differences relative to the 2000 baseline data for (a) 2001, (b) 2002, (c) 2004 and (d) 2007, determined for stratigraphic interval 3 (Watrous–Mississippian) which resides immediately above the reservoir. The labels A and B identify vertical water injection wells that are discussed in the text. The area shown is the same as in Fig. 7. Superposed are the locations of vertical wells (dots), horizontal CO2 injection wells (red lines) and horizontal production wells (black lines). The individual seismic bins (corresponding to pixels in the image) have dimensions of 40 m × 40 m.
Fig. 9. Seismic interval travel time differences relative to the 2000 baseline data for 2004 for stratigraphic intervals 1–4. The area shown is the same as in Fig. 7. Superposed are the locations of vertical wells (dots), horizontal CO2 injection wells (red lines) and horizontal production wells (black lines). The individual seismic bins (corresponding to pixels in the image) have dimensions of 40 m × 40 m.
sensitivities (dT/dz) determined by fluid substitution for the various geological intervals (see Fig. 5). For an observed interval delay time T, the corresponding thickness of the zone where CO2 substitution has occurred is zsat = T/(dT/dz), from which an estimate of the CO2 mass MCO2 is given by MCO2 = A · zsat · ϕ · SCO2 · CO2 where A is the area over which the travel time delay is observed, ϕ is the fractional porosity of the geological interval, SCO2 is the CO2 saturation within the zone, and CO2 is the CO2 density determined for the appropriate pressure and temperature conditions. The resultant proportions of CO2 determined in this manner are shown in Fig. 9. These can be used along with the known CO2 net injection quantities (i.e., total CO2 injected minus CO2 recycled from oil production) to estimate the CO2 quantities within the various stratigraphic intervals as shown in Fig. 10. The seismic-based apportionment of CO2 (Fig. 9) shows 62–70% of the net injected CO2 residing within the reservoir-containing interval after one year of injection, increasing to 92–94% after 7 years. The estimated proportion of CO2 in the zone that immediately overlies the reservoir decreases from 30–36% after 1 year of injection to 5–6% after 7 years of injection. As discussed below, the CO2 estimates for this zone are likely significantly inflated by the contribution of pressure effects in the seismic travel time maps. The maximum estimated proportion of CO2 residing in either interval above the regional seal is ≤1%. Applying these relative amounts to
Fig. 10. Seismic-based relative mass estimates by stratigraphic interval. The range of values indicated for intervals 3 and 4 correspond to estimates for the reservoir interval in which it is assumed that CO2 replaces the oil in either the Marly or Vuggy units. As noted in the text, pressure effects substantially inflate the CO2 estimates for interval 3 which lies immediately above the reservoir.
the known net injected quantitites of CO2 (Fig. 10) indicates that the maximum amount of CO2 that could reside above the regional seal by 2007 is <56,000 tonnes. 4.6–4.7 Mtonnes are attributed to the reservoir and below, and a maximum of 0.31 Mtonnes potentially resides just above the reservoir.
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
S101
6. Discussion The interval travel time differences provide a systematic means of apportioning quantities of CO2 to the various geological intervals. However, there are limitations that must be appreciated. These limitations include uncertainties in the rock-physics model used to calculate the seismic sensitivities to CO2 , the difficulty of assessing very low levels of CO2 from seismic data that is inherently nonrepeatable, the uncertainties of how CO2 is distributed within a particular geological interval, the subjectivity in choosing a threshold level for the data, characterization of the subsurface using logs from a limited number of wells, and the assumption that all travel time delays are related to pore fluid substitution (as opposed to pressure changes, for example). As noted, a threshold was applied to the interval travel time differences in an attempt to eliminate small travel time variations that result from “noise” in the data. As a result, bona fide CO2 -related travel time changes that are less than 0.5 ms may be missed in the analysis. This will correspond to CO2 that is unaccounted for by the seismic data. The estimation of CO2 quantities demonstrated in this study provides a practical means of establishing an inventory of stored CO2 subject to the limitations listed above. For example, the ability to state that 1.0% is a maximum bound for CO2 potentially residing above the regional seal is useful for this purpose. However, in regard to detecting whether CO2 is leaking from the reservoir, this result is of limited value as it does not completely rule out the possibility of leakage. This is a consequence of the time-lapse seismic data being inherently non-repeatable and having limited sensitivity. An implicit assumption of the calculation is that all of the travel time delays are ascribed to fluid substitution effects. Decreases in effective pressure associated with increased pore pressures can have a similar effect. However, there are two mitigating factors here. First, significant pore pressure increases (e.g., 5-10 MPa at the reservoir level) are required to produce significant travel time delays. In most injection scenarios, such pressure changes would be restricted to the immediate vicinity of the injection reservoir, and thus would be expected to be generally insignificant at shallower levels. Thus, pressure contributions to the travel time delays should be most significant at the reservoir level and there they are generally expected to be smaller than fluid substitution effects (e.g., White, 2013). In the Weyburn example shown here, an exception to this almost certainly occurs in the zone immediately overlying the reservoir (see Fig. 8). There, several pronounced local travel time delays are observed which are centered on vertical wells that have been used exclusively for water injection since 2000, indicating that the observed seismic anomalies are non-CO2 related. For example, consider the two vertical water injectors A and B in Fig. 8. Water injection in well B occurred at a relatively high rate (∼1.26 × 105 m3 /year) from 2000 to 2002, was then shut-in until the middle of 2004, after which it continued at a reduced rate (∼0.5 × 105 m3 /year). The related travel time anomalies in Fig. 8 increase in intensity in 2001 and 2002 during the high injection rate period, and then fade in the 2004 and 2007 maps which correspond to the period of no water injection or reduced injection. Similarly, in well A, relatively high water injection rates to 2001 were followed by no injection till 2003 when injection began at a reduced level. In this case, the most pronounced travel time anomaly occurs in the 2004 map after 2 years of continuous water injection. The prominence of these non-CO2 related travel time delays artificially increases the estimate of CO2 for this layer. Thus, the recognition that the CO2 estimates presented here are maximum bounds is especially true for geological interval 3. The inferred pressure effect (described above) assigned to the interval immediately above the reservoir explains the apparent decrease in CO2 within this zone (interval 3) after 2002 that is
Fig. 11. Mass estimates calculated by combining the relative mass estimates from Fig. 10 with known net injection quantities. The range of values indicated for intervals 3 and 4 correspond estimates for the reservoir interval in which it is assumed that CO2 replaces the oil in either the Marly or Vuggy units.
depicted in Fig. 11. This apparent decrease of CO2 might be interpreted as due to upward migration into the overlying intervals (1 and 2). However, it most likely represents depressurization of this zone over time. Furthermore, the interval travel time delays for intervals 1 and 2 show no evidence for upward migrating CO2 . 7. Conclusions Time-lapse seismic interval travel time changes have been used to estimate the proportions of injected CO2 within four stratigraphic intervals including the reservoir interval and 3 overburden intervals. Log-based fluid substitution modeling formed the basis for calculating the seismic sensitivity to the presence of CO2 within the various intervals, accounting for the effects of rock lithology, porosity, pore pressure and temperature. Travel time sensitivity values increase by an order of magnitude (0.012 ms/m to ∼0.25 ms/m) moving upward from the reservoir depth (1450 m) to the shallowest aquifer considered (600 m), indicating that travel time monitoring should be very effective for detecting the presence of CO2 above the reservoir. The seismic-based relative apportionment of CO2 shows the fraction of net injected CO2 within the reservoir and below increasing from at least 62–70% after 1 year of injection to 92–94% after 7 years of injection, whereas the proportion in the immediately overlying zone decreases from 30–36% to ∼5–6% over the same time period. The estimate for CO2 residing in the interval just above the reservoir is significantly overestimated due to the contribution of pressure effects associated with water injection. The maximum estimated proportion of CO2 residing in either interval above the regional seal is ≤1%. The maximum amount of CO2 potentially residing above the regional seal by 2007 is <56,000 tonnes, with a maximum of 0.31 Mtonnes just above the reservoir, and 4.6-4.7 Mtonnes in the reservoir and below. In the Weyburn example presented here, pressure effects are considered of secondary importance in contributing to travel time delays at the reservoir level where CO2 quantities are expected to be large. However, pressure effects are very significant in the zone just above the reservoir and greatly affect the corresponding estimation of CO2 quantities for that zone. Thus, improved seismic-based CO2 quantity estimates will require the inclusion of field-based knowledge, incorporation of constraints from flow simulations, and adoption of rock-fluid physics models with site-specific calibration.
S102
D.J. White / International Journal of Greenhouse Gas Control 16S (2013) S95–S102
Acknowledgements Cenovus provided the wells logs and the seismic data for the 2004 and 2007 surveys. Initial seismic processing was done by Dave Secord at Sensor Geophysical. Discussions with Dave Cooper provided much insight based on his extensive experience in the Weyburn field operations. Barbara Dietiker contributed to the time-lapse seismic analysis and produced the log-based property compilations. Hampson-Russell Pro4D software was used for some of the time-lapse seismic analysis and fluid substitution modelling. This is contribution number 20120441 of the Geological Survey of Canada. References Batzle, M., Wang, Z., 1992. Seismic properties of pore fluids. Geophysics 57, 1396–1408. Gassmann, F., 1951. Uber die Elastizitat poroser Medien. Veirteljahrsschrift der Naturforschenden Gesellschaft in Zurich 96, 1–23. Han, D., Sun, M., Batzle, M., 2010. CO2 velocity measurements and models for temperatures up to 200 ◦ C and pressures up to 100 MPa. Geophysics 75, E123–E129. Mavko, G., Mukerji, T., Dvorkin, J., 2003. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media. Cambridge University Press, Cambridge.
Nader, R., 2012. Effectiveness of seismic detection modelling for CO2 migation scenarios near the Aquistore CO2 injection site, Estevan, Saskatchewan. B.Sc. Thesis. Carleton University, Ottawa. Ramirez, A., White, D., Hao, Y., Dyer, K., Johnson, J. Estimating reservoir permeabilities using the seismic response to CO2 injection and stochastic inversion. International Journal of Greenhouse Gas Control Technologies, http://dx.doi.org/10.1016/j.ijggc.2012.11.031, in press. Smith, T.M., Sondergeld, C.H., Rai, C.S., 2003. Gassmann fluid substitutions: a tutorial. Geophysics 68 (2), 430–440. White, D., 2012. Geophysical monitoring. In: Hitchon, B. (Ed.), Best Practices for Validating CO2 Geological Storage. Geoscience Publishing, Sherwood Park, pp. 155–210. White, D., 2013. Seismic characterization and time-lapse imaging during 7 years of CO2 flood in the Weyburn Field, Saskatchewan, Canada. International Journal of Greenhouse Gas Control, http://dx.doi.org/10.1016/j.ijggc.2013.02.006, in press. White, D.J., Meadows, M., Cole, S., Ramirez, A., Hao, Y., Carle, S., Duxbury, A., Samson, C., Kendall, J.-M., Verdon, J.P., Dietiker, B., Johnson, J., Morozov, I., 2011. Geophysical monitoring of the Weyburn CO2 flood: results during 10 years of injection. Energy Procedia 4, 3628–3635. Whittaker, S., et al., 2004. Theme 1: geological characterization. In: Wilson, M., Monea, M. (Eds.), IEA GHG Weyburn CO2 Monitoring and Storage Project Summary Report 2000–2004. From the Proceedings of the 7th International Conference on Greenhouse Gas Control Technologies. Petroleum Technology Research Centre, Regina, pp. 1–72 (Chapter 1). Wood, A.W., 1955. A Textbook of Sound. McMillan Co., New York.