Modeling of time-lapse seismic monitoring using CO2 leakage simulations for a model CO2 storage site with realistic geology: Application in assessment of early leak-detection capabilities

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Contents lists available at ScienceDirect

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

Modeling of time-lapse seismic monitoring using CO2 leakage simulations for a model CO2 storage site with realistic geology: Application in assessment of early leak-detection capabilities

T

⁎

Zan Wanga, , William P. Harberta,b, Robert M. Dilmorec, Lianjie Huangd a

ORISE, National Energy Technology Laboratory, United States Department of Energy, Pittsburgh, PA 15236, USA Department of Geology and Environmental Science, University of Pittsburgh, Pittsburgh, PA 15260, USA c National Energy Technology Laboratory, United States Department of Energy, Pittsburgh, PA 15236, USA d Los Alamos National Laboratory, MS D452, Los Alamos, NM 87545, USA b

A R T I C LE I N FO

A B S T R A C T

Keywords: Carbon sequestration Leak detection Seismic monitoring Seismic modeling Rock physics modeling

Time-lapse surface seismic surveys have been widely used at carbon sequestration sites for site characterization, monitoring subsurface CO2 plume migration, and detecting potential CO2 leakage from a storage reservoir. Monitoring in the first permeable unit directly above the primary seal is important for early detection of CO2 leakage. Forward modeling of time-lapse seismic data can be used to assess the utility of the seismic method for CO2 leakage detection. We develop a workflow for forward modeling of time-lapse seismic data, including constructing seismic velocity models using flow simulation outputs, modeling of pre-stack and post-stack synthetic seismic data following seismic data processing sequence and analysis of processed synthetic time-lapse seismic data. We apply the forward modeling and analysis workflow to assessing the detectability and the earliest detection time of seismic monitoring using the hypothetical CO2 leakage scenarios for a model geologic storage site with realistic geology. We derive the detection thresholds using the simulated normalized rootmean-square (NRMS) differences for the no-leakage case at a range of signal-to-noise ratios, representing the background noise levels in seismic data. We then compare NRMS differences triggered by the CO2 leakage to the detection thresholds at each time step to quantify the detectability and the earliest detection time of seismic monitoring. We analyze the effects of the acquisition parameters and elastic parameters on the produced synthetic seismic data and earliest detection time. Our modeling results indicate that high signal-to-noise ratio is needed to detect the CO2 leakage at the model site. Minimizing the background noise in seismic data is crucial for improving the detectability of the seismic method. Increasing the shot density or increasing the dominant frequency of the source wavelet is likely to increase the possibility of the leakage detection and reduce the time needed for the detection. The elastic parameters used in the rock physics modeling have significant effects on the resultant seismic velocity models and synthetic seismic data, highlighting their critical roles in understanding and interpreting time-lapse seismic reflection data associated with CO2 monitoring, verification and accounting activities.

1. Introduction Carbon sequestration in deep subsurface formations, such as saline aquifers, oil and gas reservoirs, and unmineable coal seams, is considered to be a promising method to reduce carbon dioxide (CO2) emissions and combat global warming. Monitoring for possible leakage is required to assure long-term secure containment of CO2 (U.S. EPA, 2016). A variety of monitoring techniques have been deployed and tested at CO2 storage sites around the world (Chadwick et al., 2006; Daley et al., 2013, 2007; Ivandic et al., 2015; Meckel et al., 2008).

⁎

Corresponding author. E-mail address: [email protected] (Z. Wang).

https://doi.org/10.1016/j.ijggc.2018.06.011 Received 12 January 2018; Received in revised form 25 May 2018; Accepted 14 June 2018 1750-5836/ © 2018 Elsevier Ltd. All rights reserved.

Different monitoring technologies are used for different applications and can provide either direct or indirect measurements of the injected CO2. For example, soil CO2 flux and groundwater sampling are direct measurements focused on the shallow subsurface (Oldenburg et al., 2003; Yang et al., 2011). Seismic monitoring methods, including surface seismic surveys, vertical seismic profiling and cross-well seismic monitoring, can provide indirect measurements of geophysical parameters in the deep subsurface (Chadwick et al., 2004; Cole et al., 2002; Silver et al., 2007; Souza and Lumley, 2015; Vasco et al., 2014; Zhou et al., 2010). A number of geophysical monitoring technologies that are

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

difference of 7%, after applying acquisition noise reduction and several noise-reduction processing steps. The objective of this study is to evaluate the early leak-detection capabilities of seismic monitoring using forward seismic modeling and analysis of synthetic time-lapse seismic data. For this forward modeling and analysis, we use sets of simulations of hypothetical, pre-supposed fluid migration out of the storage reservoir. Those simulations are based on a hypothetical point source CO2 leakage into the first permeable unit above the primary seal. We present a workflow for forward modeling of seismic data, including constructing seismic velocity models using multi-phase flow simulation results, generating synthetic common-shot gathers followed by a basic processing sequence and analysis of stacked/migrated synthetic seismic data. The time-lapse NRMS differences caused by CO2 leakage are compared against the seismic noise thresholds to assess the detectability of seismic signal changes and the earliest detection time. We derive the detection thresholds using simulated NRMS differences for the no-leakage case at a range of signalto-noise ratios, representing the background noise levels in processed seismic data. In addition, we study the effects of the acquisition parameters and elastic parameters on the time-lapse seismic data and the estimated earliest detection time.

applicable to assessing regulatory compliance and detecting potential leakage are reviewed in Harbert et al. (2016). Early detection of CO2 leakage, should it occur, is important in both the injection and the post-injection phases of a carbon sequestration project. Leakage detection in the first permeable unit directly above the primary seal has drawn attention for providing early warning of possible leaks (Chabora and Benson, 2009; Jung et al., 2013; Namhata et al., 2017; Wang and Small, 2014). Leakage simulations at the FutureGen 2.0 site – a site near Jacksonville, Illinois that has been well characterized as a potential deep saline storage – focus on the first permeable unit (Ironton Sandstone) above the primary seal (Vermeul et al., 2016; Williams et al., 2014). The model setting of those multiphase flow simulations was designed for evaluating the early leak-detection capabilities of the monitoring system (Williams et al., 2014). Time-lapse surface seismic surveys have been widely used at carbon sequestration sites as an essential method for site characterization, imaging subsurface CO2 plume migration and detecting potential CO2 leakage from the storage reservoir (Chadwick et al., 2009; Ivanova et al., 2012; White, 2013). Arts et al. conducted forward acoustic and elastic seismic modeling of synthetic seismic data at the Sleipner CO2 storage site and compared with field time-lapse seismic data to derive CO2 plume distributions in the reservoir (Arts et al., 2009). Seismic velocities used in forward seismic modeling are calculated as a function of CO2 saturation using the Gassmann equations (Gassmann, 1951), and pressure effects on the seismic velocities are often neglected (Arts et al., 2009). The Gassmann-Biot fluid substitution theory (Biot, 1941; Gassmann, 1951) has been commonly employed in estimating changes in the seismic response caused by fluid injection or leakage (Roach et al., 2015; Souza and Lumley, 2015; Vasco et al., 2014; Wang and Small, 2015) and is widely used in seismic reflection data analyses in hydrocarbon exploration, including using CO2 for enhanced oil recovery (Harbert and Lipinski, 2010). These quantitatively estimated changes in rock properties directly result in changes in synthetic time-lapse seismic data computed using forward seismic modeling. The uncertainties in the derived seismic velocities using the Gassmann equations are largely related to uncertainties in elastic parameters (Chadwick et al., 2004). Roach et al. (2015) simulated changes in seismic response of rock caused by replacement of brine by CO2 at the Aquistore site, and found that analysis of reflection amplitude differences between the baseline and the repeat surface 3D seismic surveys is a powerful tool for detecting small amounts of free-phase CO2 in deep formations. CO2-induced changes in seismic responses of reservoir rocks have been commonly inferred using time-lapse seismic data. To make this method effective, the background noise level in a repeat time-lapse seismic dataset that is not due to changes in the subsurface geology should be quantified (Kragh and Christie, 2002). The actual size of the CO2 plume that can be detected using seismic monitoring largely depends on the noise level in the time-lapse seismic data and the thickness of the perturbed geologic interval contributing to changes in measured seismic signals (Lüth et al., 2015). The noise in seismic data can be attributed to the natural background random noise at the site, which is irreducible, and noise introduced in acquisition and processing steps, which can be suppressed through various techniques such as precise positioning of sources and receivers, and dedicated processing methods (Lumley et al., 2003; Meunier et al., 2001; Roach et al., 2015). Specific criteria and estimates of the risk of failure of 4-D seismic reflection surveys to properly image subsurface targets have been presented in Lumley et al. (1997). In a conformance study at the Ketzin pilot site (Lüth et al., 2015), the noise threshold is calculated as the 95th percentile value of the cumulative distribution of normalized noise amplitude differences, obtained from noise sections in time-lapse 3D seismic data. Roach et al. (2015) determined the background noise level at the Aquistore CO2 storage site using two surface 3D seismic surveys acquired before the start of injection. The noise threshold of real seismic data is found to be a normalized root-mean-square (NRMS)

2. Workflow of forward modeling and analysis of synthetic seismic data A workflow is detailed and demonstrated using a hypothetical CO2 storage site with realistic geology based loosely on the FutureGen 2.0 site. Previously-published simulations (Williams et al., 2014) are used to illustrate the process of forward modeling and analysis of synthetic seismic data for assessing the early leak-detection capabilities of seismic monitoring. 2.1. Description of the hypothetical CO2 leakage simulations A set of hypothetical leakage simulations were developed to support leak detection capability and risk assessment at the FutureGen 2.0 site. These simulations were based on a hypothetical point source CO2 leakage into the first permeable unit (Ironton Sandstone) above the caprock. The bottom of the model domain is the bottom of the Ironton Sandstone and the top of the model domain is within the unit above the lowermost groundwater aquifer. A CO2 point leakage source is applied at the center bottom of the model domain. The storage reservoir, the cap rock and the leakage pathway are not included in the leakage model. The simulated leakage scenarios assume that 1% of the total planned injected supercritical CO2 mass (22 MMT) has leaked over 20, 100 and 500 years, respectively. More detailed description of the FutureGen 2.0 leakage model can be found in Williams et al. (2014). In this study, we use the simulation outputs for the 20-year leakage scenario in which 1% of total injected CO2 leaked over 20 years (0.011 MMT/yr). The output files for the 20-year leakage scenario contain 31 time steps from t = 0 to t = 20 years. CO2 migration through the primary seal of this magnitude is not expected for the FutureGen 2.0 site or other well characterized geologic carbon storage sites. It is emphasized that this simulation presupposes a large leak event only for purposes of evaluating whether or not such leakage can be detected. 2.2. Rock physics modeling We first construct seismic velocity models at each simulation time step for computing synthetic time-lapse seismograms. We build the initial seismic velocity model (i.e., before CO2 leakage) using the processed wireline log data (measurements of the P- and S-wave slowness) from the initial stratigraphic borehole (Alliance, 2012). The slowness measurements are averaged over the interval of each stratigraphic unit using the arithmetic averaging method. The depth intervals of the stratigraphic units at the model site are obtained from the site 40

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 1. Background P- and S-wave velocity models constructed based on reported stratigraphic borehole data (Alliance, 2012).

mineral grains and effective pore fluid (mixed pore fluid phases), respectively. The bulk modulus and density of effective pore fluid (Kfl and ρfl ) can be estimated using inverse bulk modulus averaging, a parallel law or Wood’s law (Wood, 1955) and arithmetic averaging of densities of the separate fluid phases (brine phase and supercritical CO2 phase), respectively (Kumar, 2006):

characterization report (Alliance, 2015). The beginning depth of the wireline log data is 172 m below the ground surface. For the intervals above the beginning depth, we assume that the slowness is constant at the slowness value of the first interval where wireline log data are available. The seismic velocity for each stratigraphic unit is calculated as the reciprocal of the averaged slowness. The grid and zonation definition, including the inactive nodes specification, in the leakage simulations by Williams et al. (2014) is used to characterize the dip along the bottom of the stratigraphic unit. Fig. 1 shows the constructed background P- and S-wave velocity models. For those simulated CO2 leakage scenarios, we calculate changes in the seismic velocity caused by CO2 leakage from the storage reservoir at each time step using the Gassmann-Biot fluid substitution theory (Biot, 1941; Gassmann, 1951) and the Hertz-Mindlin contact theory (Mindlin, 1949). The instantaneous velocity field in the model domain is related to the moduli and density of the rock by:

Vp =

Ksat + 4μ/3 , ρsat

Vs =

μ , ρsat

Sg 1 Sw = + and ρfl = Sw ρbrine + Sg ρco2 , Kfl Kbrine K co2 Where Sg is the CO2 saturation (as a volume fraction) and Sw (= 1−Sg ) is the brine saturation. The bulk moduli and densities of pure brine and supercritical CO2 (Kbrine, K co2 and ρbrine , ρco2 ) are calculated as a function of temperature, pore pressure and salinity using relationships developed by Batzle and Wang (1992). The input parameters (i.e., temperature, pore pressure, salinity and CO2 saturation) at each node in the model domain at each time step are obtained from the leakage simulation outputs (Williams et al., 2014). Example outputs at t = 20 years in a 20-year leakage scenario are shown in Fig. 2. The bulk and shear moduli of mineral grains (Kmineral and μmineral ) are estimated by taking Voigt-Reuss-Hill averaging (Hill, 1952) of the mineral constituents:

Where Vp and Vs are the P-wave and S-wave velocities; Ksat is the bulk modulus of the rock after fluid substitution; μ is the shear modulus of the rock, which is constant during the fluid substitution in the Gassmann-Biot theory. The density of rock, ρsat , after fluid substitution is calculated using:

n

k voigt =

Ksat = Kframe +

φ Kfl

+

Kframe

1−φ Kmineral

)

1 ; Kmineral (or μmineral ) n ∑i = 1 vi/ ki

= 0.5 *(k voigt + kreuss ), Where vi is the volume fraction and ki is the bulk (or shear) modulus of mineral constituent i; and n is the number of mineral constituents, which is 7 in this study. We obtain the mineral composition of each stratigraphic unit (volume fraction) from the processed wireline log data acquired during site characterization. Elemental analysis of the lithology log indicates that the mineral constituents at the model site are quartz, dolomite, calcite, K-feldspar, kaolinite, illite and pyrite. The moduli of the mineral constituents, which can be obtained from published literature and textbooks such as Mondol et al. (2008); Mavko et al. (2009) and Roach et al. (2015), are summarized in Table 1. The bulk and shear moduli of mineral grains remain constant during the Gassmann fluid substitution. initial The initial bulk modulus of the dry rock frame (K frame ) can be derived from the wireline log data by rewriting the Gassmann equation for Kframe (Zhu and McMechan, 1990) as

where ρinitial is the initial bulk density of the rock, which can be obtained from the density log acquired during site characterization; ρbrine and ρfl are the densities of the fluid phase before (i.e., pure brine) and after CO2 invasion (i.e., mixture of brine and CO2), respectively; φ is the porosity of the rock as a volume fraction, which can be obtained from the wireline log data. In our study, we use the same porosity values as used in the flow simulations by Williams et al. (2014), which were derived based on the characterization data from the initial stratigraphic borehole at the model site (Vermeul et al., 2016). We estimate the bulk modulus of the rock after fluid substitution (Ksat ) at each node in the model domain using the low frequency Gassmann equation (Gassmann, 1951): Kmineral

kreuss =

i=1

ρsat = ρinitial + φ (ρfl − ρbrine ),

(1 −

∑ vi * ki;

2

initial ⎛ Ksat ⎝ =

⎜

Kframe

− K2

initial K frame

mineral

φKmineral K initial fluid

φKmineral K initial fluid

Where Kframe , Kmineral , and Kfl are the bulk moduli of the dry rock frame, 41

+

+ 1−φ⎞−Kmineral ⎠ , ⎟

initial Ksat Kmineral

− 1−φ

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 2. Example outputs from hypothetical leakage simulations (Williams et al., 2014): a) temperature; b) pore pressure; c) salt concentration; d) CO2 saturation profiles for a x–z cross section at y = 0 location in the leakage zone and the lower-permeability unit directly above the leakage zone at t = 20 years in a 20-year leakage scenario.

pore pressure and the confining pressure from the leakage simulation outputs and the wireline log data acquired during site characterization, respectively. Parameter pratio denotes the Poisson’s ratio of the minerals; φc is the critical porosity, above which the rock ceases to be in a consolidated, frame supporting state and becomes a fluid-supported liquid suspension. The reported critical porosity of sandstone from Nur et al. (1998) is 0.4. Parameter C represents the average number of contacts per spherical mineral grain. According to published data (Carcione et al., 2006; Murphy, 1982), C = 2.8/φc . The effective bulk and shear moduli of dry rock frame at a different porosity φ can be estimated using the modified Hashin-Strikman lower bound (Dvorkin and Nur, 1996) based on the original Hashin-Strikman lower bound (Hashin and Shtrikman, 1963).

Table 1 Moduli of rock mineral constituents. Mineral

Bulk modulus (GPa)

Shear modulus (GPa)

Quartz Dolomite Calcite K-feldspar Kaolinite Illitea Pyrite

37 94.9 76.8 37.5 55 61.5 147.4

44 45 32 15 31.8 41.1 132.5

a

Illite is modeled as muscovite.

initial Where Ksat is the initial bulk modulus of the wet rock, which can be derived from the wireline log data as:

φ

Kframe

4 initial Ksat = ρinitial ⎡ (V Pinitial )2− *(V sinitial )2] ⎢ 3 ⎣

φ

1/3

C2 (1−φc )2μ 2mineralpd ⎞ Kmc = ⎜⎛ 2 2 ⎟ ⎝ 18π (1−pratio) ⎠

φ

−1

Where Kmineral and μ mineral are the bulk and shear moduli of mineral grains, respectively. Table 2 summarizes the parameters obtained from the wireline log data and from the leakage simulations for constructing the initial background seismic velocity model and the velocity models at each time step in the leakage simulations. The last column in Table 2 shows representative values for the layer where the simulated CO2 leakage

1/3

2 5−4 * pratio ⎛ 3C2 (1−φc )2μ mineralpd ⎞ ⎜ ⎟ 2 5(2−pratio) ⎝ 2π (1−pratio)2 ⎠

−1

1− φ ⎡ ⎤ φc c ⎥ μ=⎢ + μ mc 9Kmc + 8μ mc μ mc 9Kmc + 8μ mc ⎢ μ mc + ( K + 2μ ) μ mineral + 6 ( K + 2μ ) ⎥ 6 mc mc mc ⎦ mc ⎣ μ mc ⎛ 9Kmc + 8μ mc ⎞ − ⎜ ⎟ 6 ⎝ Kmc + 2μ mc ⎠

The effects of pore pressure changes on bulk modulus of dry rock frame are modeled based on the Hertz-Mindlin contact theory (Mindlin, 1949). The effective elastic properties of rock at critical porosity φc are given by:

μ mc =

φ

1− φ ⎤ ⎡ φc 4 c =⎢ + ⎥ − μ mc, 4 4 3 Kmineral + 3 μ mc ⎥ ⎢ Kmc + 3 μ mc ⎦ ⎣

,

Where pd = pc −p is the differential pressure calculated by subtracting the pore pressure ( p ) from the confining pressure ( pc ). We obtain the 42

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Table 2 Parameters needed for constructing seismic velocity models.

Table 3 Acquisition parameters used in the acoustic seismic modeling.

Parameters from log

Units

Representative values

Parameter

Value

Initial P-wave velocity Initial S-wave velocity Initial bulk density of the rock Confining pressure Mineral composition of the rock

m/s m/s g/cm3

4673 2764 2.47

MPa volume fraction

24.35 0.828 quartz, 0.056 dolomite, 0.01 calcite, 0.03 K-feldspar, 0.064 kaolinite, 0.008 illite and 0.004 pyrite

Number of shots Number of receivers Shot point spacing Receiver station spacing Source wavelet Dominant frequency of the source wavelet Sampling rate Record length

54 299 25 m 5m Ricker 30 Hz 1 ms 3s

Parameters from flow simulations

Units

Representative values

Temperature Pore pressure Salinity of brine CO2 saturation Porosity

°C MPa weight fraction volume fraction volume fraction

33.5 10.14 0.000597 0.7 0.12

occurs. Fig. 3 shows the predicted changes in seismic velocities from the rock physics modeling at t = 20 years in a 20-year leakage scenario.

mean-square velocity model (Margrave, 2003) converted from the Pwave instantaneous velocity model (e.g., Fig. 3a). No gain or mute is applied to the synthetic common-shot gathers in creating the CMP stack. 2D depth migration is performed on the stacked section using the P-wave instantaneous velocity model (e.g., Fig. 3a) with the application of the phase shift plus interpolation method implemented in the CREWES MATLAB Toolbox (Ferguson and Margrave, 2005), to produce the migrated 2D seismic section.

2.3. Forward seismic modeling

2.4. Analysis of synthetic seismic data

We perform two-dimensional (2D) acoustic-wave modeling to assess the CO2 leakage detectability. A 2D surface seismic line runs through the constructed P-wave velocity model at each time step. The 2D line is parallel to the x axis, at y = 0 location so that the imaged x–z cross section intersects the leakage point. The extent of the seismic velocity model in x and z directions are from 0 to 1490 m, and 0 to 1190 m, respectively. The bottom of the leakage zone is at a depth of about 1035 m, which is also the bottom of the multi-phase flow model in the leakage simulations. The seismic velocity model is extended at the bottom of the leakage zone by about 155 m to avoid boundary effects. Along the 2D surface seismic line, there are 54 shots covering the x range of 80 m to 1405 m and 299 receivers covering the x range of 0 m to 1490 m. The shot point spacing is 25 m. The receiver spacing is 5 m. The dominant frequency of the source wavelet, which is modeled as a Ricker wavelet, is 30 Hz. The sampling rate is 1 ms and the record length is 3 s. A summary of the acquisition parameters is given in Table 3. We generate synthetic common-shot gathers using a 2D acoustic finite-difference modeling code in the CREWES MATLAB Toolbox (Margrave, 2000). The computation time step in the acoustic finitedifference wave propagation simulation is 0.2 ms and the computation grid is 5 m * 5 m. The generated common-shot gathers are converted to the common mid-point (CMP) gathers and the CMP stack is created after applying the normal move out (NMO) correction using the root-

The post-stack migrated synthetic seismic data at each time step of CO2 leakage are subtracted from the migrated synthetic data at the initial time step (t = 0), resulting in the time-lapse amplitude difference datasets. To make the amplitude difference values easier to interpret, the percent amplitude changes are calculated as the amplitude difference values divided by the peak amplitude value at the initial time step. At each time step, the maximum, mean and standard deviation of the absolute value of the percent amplitude change are calculated. The repeatability metric: normalized root-mean-square (NRMS) difference along the 2D seismic line at each time step is calculated using the equation (Kragh and Christie, 2002):

NRMS =

200 × d2

d2

∑d1 (tracet −tracet 0)2 / N

∑d1 tracet 2/ N +

d2

,

∑d1 tracet 0 2/ N

where tracet and tracet0 are the traces at time step t and at the initial time step in the post-stack depth migrated seismic datasets. Parameters d1 and d2 are the start and end depths of the desired depth window. Here N is the number of samples per trace within the desired window. The calculated NRMS is in percentage, which quantifies the likeness of two traces. The lower the NRMS values are, the more repeatable the two traces are. The NRMS value is calculated at each trace location x at each time step. The maximum NRMS value at each time step is obtained, which is used as the main variable for inferring CO2 leakage in

Fig. 3. Predicted changes in seismic velocities: a) P-wave velocity changes; b) S-wave velocity changes at t = 20 years in a 20-year leakage scenario. 43

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 4. Workflow of forward modeling and analysis of synthetic seismic data.

this study. Fig. 4 is a schematic flow chart showing the workflow of forward modeling and analysis of synthetic seismic data. The outcomes of the forward modeling and analysis workflow are used to infer the CO2 leakage, which is a combined effect of variations in geophysical parameters, including pore fluid type, pore fluid saturation and effective pressure. 2.5. Definition of the background noise level and the earliest detection time Quantifying the background noise level is crucial as the seismic signals induced by CO2 leakage must be larger than the seismic background noise to be detectable. When real seismic data are available, the background noise level (or detection threshold) can be derived by identifying a noise section of seismic traces with low amplitude difference values in the seismic dataset and calculating a practical threshold, such as the 95th percentile from the cumulative distribution of amplitude difference values for the noise area (Lüth et al., 2015). Or if two baseline seismic surveys are available at the site, the background seismic repeatable noise can be determined by calculating the repeatability metric, NRMS difference using the processed seismic data from the two baseline surveys (Roach et al., 2015). In our study, we derive the detection threshold values from synthetic seismic datasets generated at a range of signal-to-noise (S/N) ratios from 30 to 50 for the baseline seismic surveys. First, we simulate the migrated seismic data without noise at the initial condition (t = 0) using the background P-wave velocity model (Fig. 1a) and the approach described in Section 2.3. Then we add the normally distributed random noise to the synthetic seismic traces at a specified signal-to-noise ratio. A repeat data set is generated from the specified distribution for the added noise, representing the second baseline seismic survey. Then the detection threshold is calculated as the maximum NRMS difference of the simulated two baseline seismic datasets, which quantifies the NRMS difference caused by the background variations (noise). For each signal-to-noise ratio, we simulate the detection threshold 1000 times with different samples for the added noise, and use the mean of the simulated detection thresholds as the final detection threshold at the specified signal-to-noise ratio. Fig. 5 shows the

Fig. 5. Detection threshold of the NRMS value (%) as a function of the signal-tonoise ratio in the migrated synthetic seismic data.

calculated detection threshold as a function of the signal-to-noise ratio in the migrated synthetic seismic data generated at t = 0. The detection threshold of NRMS value decreases from 8.2% at S/N = 30 to 4.9% at S/N = 50. The NRMS differences caused by CO2 leakage need to be larger than the NRMS differences caused by background noise to be detectable. For example, if the maximum NRMS differences caused by CO2 leakage is 7%, the CO2 leakage is detectable if S/N ≥ 36. According to Lumley (2010), NRMS values less than 20% are considered excellent for 4D seismic techniques. Roach et al. (2015) indicates that NRMS value of 7% can be achieved using permanent receiver array and applying a sequence of seismic data processing steps. The earliest detection time can be derived by comparing the maximum NRMS value of the time-lapse seismic datasets with the NRMS detection threshold. It is defined as the time needed for the maximum NRMS difference caused by CO2 leakage to exceed the detection threshold at the specified signal-to-noise ratio in migrated seismic data.

44

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 6. Example plots of the post-stack, post-migration synthetic seismic data. a) Migrated depth section at t = 0 (i.e., the no leakage case); b) Percent amplitude changes at t = 20 years (after the leakage initiation); c) Amplitude difference map at the bottom of the model domain at t = 20 years; and d) NRMS in the depth interval of 600–1170 m at t = 20 years.

3. Results

NRMS difference in the depth interval of 600–1170 m is 3.82%, which occurs at x = 745 m and t = 0.9 years. The NRMS difference caused by the seismic repeatable noise should be smaller than 3.82% to detect the CO2 leakage at t = 0.9 years. For the NRMS thresholds considered in this study (Fig. 5), all threshold values are larger than 3.82%, which means that the NRMS differences resulting from the CO2 leakage are within the background noise levels considered. Higher signal-to-noise ratio or lower NRMS detection threshold is needed to be able to detect the CO2 leakage.

3.1. Assessment of early leak-detection capabilities Fig. 6 shows example plots of the synthetic seismic data. Data in Panel a) shows the post-stack, post-migration synthetic seismic data in the whole model domain at t = 0 (before the start of leakage). The percent amplitude changes at t = 20 years (after the leakage initiation) in the whole model domain are plotted in Panel b). The affected area by CO2 leakage can be identified in Panel b) as the region with relatively high amplitude changes. Based on the flow simulation results of the model site, most CO2 plume is trapped within the leakage zone with its top at about 1011 m below the ground surface (Fig. 2d). The depth and the lateral extent of the region with high amplitude changes in the synthetic seismic data are consistent with the CO2 plume area predicted in the leakage flow simulations. The amplitude difference data at the bottom of the model domain where the simulated CO2 leakage occurs at t = 20 years are shown in Panel c). Panel d) shows the calculated NRMS in the depth interval of 600–1170 m along the 2D seismic line. The peak NRMS difference at t = 20 years is about 3.5%, compared with the traces at t = 0 (i.e., the no leakage case). The peak value occurs at the location x = 745 m, where the leakage point is deployed. The mean, standard deviation, maximum of the percent amplitude changes and the maximum of the NRMS difference in the migrated synthetic seismic data are plotted as a function of time in Fig. 7. The solid blue curves denote the results using the nominal parameter values. The mean and the standard deviation of the absolute amplitude changes increase with increasing time. The rate of increase is larger in the first two years than that after the first two years. The maximum value of the absolute amplitude changes reaches the peak value (3.2%) at t = 0.9 years. It gradually decreases after four years of leakage. The maximum

3.2. Effects of acquisition parameters In this section, we evaluate the effects of shot point spacing and dominant frequency on the computed seismic data and the results of the time-lapse amplitude analysis. The nominal setting of the acquisition geometry contains 54 shots, with the shot point interval of 25 m. Here we increase the total number of shots to 135 by reducing the shot point spacing to 10 m. The 135 shot points cover a wider range from 80 m to 1420 m in the x direction. The other acquisition parameters are kept the same as the nominal values defined in Table 3. Fig. 8 shows example plots of the post-stack migrated seismic data using the high-density acquisition geometry. The resolution of the produced seismic images is higher than that resulting from use of the nominal setting of the acquisition geometry. The maximum NRMS difference at t = 20 years is 4.35%, which is higher than the value in the nominal setting. The mean, standard deviation, maximum of the absolute amplitude changes and the maximum of the NRMS values in the time-lapse seismic datasets for the high-density shot acquisition are plotted as dashed orange curves in Fig. 7. The trend of these statistics of the absolute amplitude changes using the high-density acquisition geometry is similar to that obtained using the nominal acquisition geometry. The 45

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 7. Statistics of the percent amplitude changes and the NRMS values in the time-lapse migrated synthetic seismic datasets: a) mean of the absolute amplitude changes; b) standard deviation of the absolute amplitude changes; c) maximum of the absolute amplitude changes; and d) maximum of the NRMS values in the depth interval of 600–1170 m. The solid blue, dashed orange, dash-dot green, and dotted red curves represent the results using the nominal parameter setting, 135 shots, 60 Hz source wavelet and 90 Hz source wavelet, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

cumulative CO2 mass leaked at 0.5 years is 0.0055 MMT. To evaluate the effects of the dominant frequency of the source signal on the modeling results, we increase the dominant frequency from 30 Hz to 60 Hz and 90 Hz, keeping all other acquisition parameters at the nominal values defined in Table 3. Corresponding simulation results are compared with the base case results in Section 3.1. Fig. 10 shows the percent amplitude changes in the whole model domain, the amplitude difference values at the bottom of the model domain and the NRMS differences at t = 20 years for the dominant frequency of 60 Hz and 90 Hz, respectively. The distribution of the amplitude anomalies using the high-frequency source wavelet is different from that in the base case with the dominant frequency of 30 Hz. The large amplitude anomalies can be identified at two ranges of x locations around the leakage point in the high-frequency cases, while the large amplitude anomalies are more concentrated at regions above the leakage point at dominant frequency of 30 Hz. The resolution along the depth in the z direction improves as the dominant frequency increases. At the dominant frequency of 60 Hz and 90 Hz, the maximum of the NRMS differences at t = 20 years are 5.05% and 5.63%, respectively, both occurring at the location x = 570 m. The mean, standard deviation, and maximum of the absolute amplitude changes and the maximum of the NRMS differences as a function of time of leakage at dominant frequencies of 60 and 90 Hz are illustrated as dash-dot green and dotted red curves in Fig. 7. For percent amplitude changes, the mean, the standard deviation and the maximum values at high frequencies (60 Hz and 90 Hz) are larger than those at the nominal frequency (30 Hz), with the 90 Hz case providing the highest mean value and the 60 Hz case having the largest standard deviation and highest percent amplitude change. The behavior of the NRMS value

mean and the maximum values of the absolute amplitude changes at each time step are smaller than the values in the nominal setting. The maximum NRMS in the time-lapse amplitude datasets for the highdensity acquisition geometry is 5.24% (dashed orange curve in Fig. 7d), achieved at t = 0.9 years and x = 745 m. The standard deviation of the absolute amplitude changes using the high-density acquisition geometry is also smaller than that in the nominal setting, meaning that the absolute amplitude changes are more concentrated around their mean values. The NRMS difference using the high-density acquisition geometry is significantly higher than that using the nominal acquisition geometry. The lowest signal-to-noise ratio in the migrated seismic data needed to detect the CO2 leakage and the earliest detection time using the highdensity acquisition geometry can be derived by comparing the dashed orange curve in Fig. 7d with Fig. 5. The NRMS detection threshold needs to be smaller than the peak NRMS value caused by CO2 leakage (5.24%) to be able to detect the leakage at t = 0.9 years. The corresponding signal-to-noise ratio for the critical detection threshold is 47, meaning that the signal-to-noise ratio needs to be larger than or equal to 47 to be able to detect the CO2 leakage. Fig. 9 shows the computed earliest detection time as a function of the signal-to-noise ratio in the migrated synthetic data, starting from the critical signal-to-noise ratio below which the CO2 leakage cannot be detected. The earliest detection time using the high-density acquisition geometry is 0.9 years (i.e., about 11 months) at the signal-to-noise ratio of 47. The corresponding cumulative CO2 mass leaked at 0.9 years is about 0.01 MMT. Below the critical signal-to-noise ratio at 47, the CO2 leakage cannot be detected. If the signal-to-noise ratio is further improved to 50, the earliest detection time will be 0.5 years (i.e., 6 months) and the corresponding 46

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 8. Example outputs from the time-lapse amplitude analysis using the high-density acquisition geometry with 135 shots. a) Migrated depth section at t = 20 years; b) Percent amplitude changes at t = 20 years; c) Amplitude difference map at the bottom of the model domain at t = 20 years; and d) NRMS in the depth interval of 600–1170 m at t = 20 years.

the earliest detection time is t = 16 years (corresponding to a cumulative mass leaked of 0.176 MMT) at frequency of 60 Hz or t = 9 years (corresponding to a cumulative mass leaked of about 0.1 MMT) at frequency of 90 Hz. 3.3. Effects of elastic parameters The elastic properties of mineral constituents are input parameters in the rock physics modeling for constructing the seismic velocity model. For clay minerals, such as Kaolinite, the elastic properties are extremely difficult to measure and the results of different studies vary widely (Mondol et al., 2008). For example, the measured bulk and shear moduli of Kaolinite by Katahara (1996) are 55.0 GPa and 31.8 GPa, respectively, which are the values used in previous sections. Alternative study of Kaolinite moduli by Vanorio et al. (2003) shows much lower values with the measured bulk and shear moduli of Kaolinite being 11.0 GPa and 6.0 GPa, respectively. The effects of the elastic moduli of the clay mineral Kaolinite on the expected magnitudes of the changes in seismic velocities caused by CO2 leakage and the forward seismic modeling results are analyzed in this section by comparing the modeling results using the lower Kaolinite moduli measurements from Vanorio et al. (2003) with the results using the nominal Kaolinite moduli from Katahara (1996). All other input parameters are constant at the nominal values (Tables 1 and 3). Fig. 11 shows the predicted P-wave and S-wave velocity changes at t = 20 years using the lower Kaolinite moduli measurements. Changes in P-wave velocity are smaller using the lower Kaolinite moduli measurements than those using the nominal parameter values (Fig. 3). The maximum of the P-wave velocity changes decreases from 64 m/s to 24 m/s. The distribution of the P-wave velocity changes is also affected by the Kaolinite moduli input. The Kaolinite moduli input has

Fig. 9. The earliest detection time as a function of the signal-to-noise ratio in the migrated seismic data. The dashed orange, dash-dot green, and dotted red lines represent the results using the high-density acquisition geometry with 135 shots, source wavelet of 60 Hz and source wavelet of 90 Hz, respectively. The CO2 leakage could not be detected using the nominal setting of the acquisition parameters at the range of signal-to-noise ratios considered. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

with increasing time is different at high frequencies than that at the nominal frequency. At high frequencies, the NRMS value fluctuates in the first two years, and gradually increases after two years of leakage. The largest NRMS difference is achieved at t = 20 years for both the 60 Hz and the 90 Hz cases. As shown by the dash-dot green line and the dotted red line in Fig. 9, the CO2 leakage can be detected at t = 20 years (corresponding to a cumulative mass leaked of 0.22 MMT) if S/N ≥ 49 at frequency of 60 Hz or S/N ≥ 44 at frequency of 90 Hz. For S/N = 50, 47

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 10. Seismic modeling results at t = 20 years using high frequency source wavelet (dominant frequency = 60 Hz and 90 Hz, respectively): a) percent amplitude changes in the whole model domain; b) amplitude difference values at the bottom of the model domain; and c) the NRMS differences in the depth interval of 600–1170 m.

the same data as the solid blue curves in Fig. 7, representing the results using nominal parameter values. The predicted amplitude changes and the NRMS differences are much smaller using the lower Kaolinite moduli than those using the nominal Kaolinite moduli.

negligible effect on the predicted S-wave velocities. The example outputs from the forward seismic modeling using the lower Kaolinite moduli are shown in Fig. 12. The percent amplitude changes and the NRMS differences are smaller, compared to those using the nominal Kaolinite moduli. The distributions of the amplitude differences and the NRMS differences have changed corresponding to changes in the seismic velocity model. The maximum NRMS difference at t = 20 years is only 0.92%, occurring at the location x = 690 m. The NRMS difference reaches the second peak at the location x = 800 m with the value being 0.90%. The center of the two peak locations is at x = 745 m where the leakage point is located. Results from the time-lapse amplitude analysis are shown in Fig. 13. The mean, standard deviation, maximum of the amplitude changes and the maximum NRMS difference as a function of time for the lower Kaolinite moduli case are plotted as dashed orange curves in Fig. 13. The solid blue curves in Fig. 13 show

3.4. Effects of fluid mixing model Wood’s law (Wood, 1955), the parallel law was used to estimate the bulk modulus of the fluid mixture (brine and supercritical CO2) in the previous sections. Here, we evaluate the effects of the fluid mixing model using Brie’s empirical pore fluid mixing law (Brie et al., 1995):

Kfl = (Kbrine−K co2)* Sw e + K co2 Where e is an empirical cementation constant for tuning. The typical value of e is 3 (Brie et al., 1995), which is used in our study for 48

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 11. Predicted changes in seismic velocities: a) P-wave velocity changes; b) S-wave velocity changes at t = 20 years in a 20-year leakage scenario using the lower Kaolinite moduli measurements from Vanorio et al. (2003).

Fig. 12. Percent amplitude changes and the NRMS differences at t = 20 years using the lower Kaolinite moduli measurements from Vanorio et al. (2003).

differences above noise was 0.9 years at S/N = 47, corresponding to a cumulative mass leaked of about 0.01 MMT. Even once an advantageous signal-to-noise ratio is achieved, it may still be challenging to identify leakage signals in a timely manner, since the current state-ofpractice for seismic monitoring calls for data collection events at a frequency that is on the order of years, and the NRMS differences from leakage events could (as in the case considered herein) increase rapidly within the first year of leakage. This highlights the value of developing high collection frequency (i.e., frequency of conducting the time-lapse seismic surveys), even continuous, seismic monitoring with high signalto-noise performance. Other monitoring techniques, such as pressure monitoring and electromagnetic monitoring, may be applied together with seismic monitoring to increase the confidence in CO2 leakage detection and achieve improvements in time to detection. We analyze both the absolute amplitude changes and the NRMS differences in the time-lapse synthetic seismic datasets. Our modeling results indicate that the NRMS difference could be a more robust measure of detectability than the absolute amplitude change, especially when there are changes in the acquisition parameters. In some cases, the maximum of the NRMS differences gets larger while the maximum of the absolute amplitude changes gets smaller. We observe this phenomena when comparing the results from the high-density acquisition geometry with those from the nominal acquisition setting (Fig. 7). Therefore, the detection threshold and the earliest detection time in our study are computed based on the NRMS difference in the time-lapse synthetic seismic datasets. In our study, the maximum change in the P-wave velocity is only 1.21%, which is calculated based on the rock type and mineral content

comparison with the Wood’s mixing law. Fig. 14a shows the predicted changes in P-wave velocity of rocks saturated with mixtures of CO2 and in-situ brine after 20 years of leakage using Brie’s fluid mixing law. The predicted velocity changes are slightly smaller using Brie’s law compared to those using Wood’s law (Fig. 3a). The predicted percent amplitude changes at t = 20 years using Brie’s law (Fig. 14b) are similar to those obtained using Wood’s law (Fig. 6b). The statistics of percent amplitude changes and the NRMS values as a function of time using Brie’s law are shown as dotted green curves in Fig. 13. The values of percent amplitude changes and NRMS differences using Brie's law are smaller compared to those using Wood’s law. The peak percent amplitude change (2.9%) and the peak NRMS difference (3.5%) occur at t = 9 years and t = 8 years, respectively using Brie’s law, while the peak values are achieved at t = 0.9 years using Wood’s law. 4. Discussion The detectability of the CO2 leakage and the earliest detection time using the time-lapse seismic data analysis largely depend on the background noise level in the processed seismic data. Noise reductions in the acquisition and processing steps are crucial as high signal-to-noise ratio in the final processed seismic data is needed to be able to detect the CO2 leakage considered at the model site. Minimizing the noise in seismic data can improve the likelihood of detecting the CO2 leakage and reduce the time needed for the leakage detection. For the best-performing scenario considered in this study (i.e., the high-density acquisition geometry), the estimated time to detection based on anomalous NRMS 49

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Fig. 13. Comparison of statistics of percent amplitude changes and the NRMS values in the time-lapse migrated synthetic seismic datasets using nominal parameter values with Wood’s law for the fluid mixture (solid blue curves as shown in Fig. 7), lower Kaolinite moduli (dashed orange curves), Brie’s fluid mixing law (dotted green curves), respectively: a) mean of the absolute amplitude changes; b) standard deviation of the absolute amplitude changes; c) maximum of the absolute amplitude changes; and d) maximum of the NRMS values in the depth interval of 600–1170 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

herein to assess the utility of proposed advanced technologies to enhance seismic signals. For example, the addition of acoustic contrast agents, such as nanomaterials, to the CO2 phase has been proposed as a way to improve CO2 plume tracking using seismic methods (Schaef et al., 2017). The workflows developed herein represent a computationally efficient, yet effective means of exploring the benefit of such proposed technologies for storage plume delineation, detection of leakage through seals, and other related applications. For purposes of comparison, we conduct a pre-stack seismic data analysis for the initial condition of 100% brine and brine/CO2 mixture saturation cases, where the CO2 saturation in the pore space is 5% and

specific at the model site. The changes in the P-wave velocity would vary for formations with different rock types or mineral contents. The acoustic properties of supercritical CO2 fluid are estimated using analytical expressions from Batzle and Wang (1992). Alternative expressions from Span and Wagner (1996) could be used to estimate the acoustic properties of CO2. Comparison of the modeling results using these two methods shows that the Span and Wagner equation of state predicts higher density and lower seismic velocity of the supercritical CO2 fluid, resulting in larger P-wave velocity contrast and smaller Swave velocity contrast for rocks before and after CO2 fluid substitution. Future studies could consider application of the method presented

Fig. 14. Results of rock physics modeling and seismic modeling at t = 20 years using Brie’s fluid mixing law: a) P-wave velocity changes; and b) percent amplitude changes. 50

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

Technology Laboratory for their support and constructive review of this work. We thank Dr. Xianjin Yang at Lawrence Livermore National Laboratory and Dr. Kai Gao at Los Alamos National Laboratory for their valuable suggestions.

70%, respectively. Details of that analysis are included in the Supporting Information section. The interface of the Ironton Sandstone and the Davis Dolomite is chosen for the pre-stack analysis. The input parameters, such as the P- and S-wave velocities and densities, are calculated for the cell at the location of x = 600 m and z = 1000 m using the rock physics modeling method described in section 2.2. The analysis shows that using pre-stack gathers, it is possible to distinguish between these pore-saturating fluids/mixtures and that longer offsets during acquisition improved identification of these various saturation conditions using pre-stack analysis methods. The present study suggests that, by neglecting the effects of noise, and by reducing the analysis to one dimension (i.e., considering only the detectability at the interface between the bottom of the dolomite and the top of the sandstone intervals), the pre-stack analysis overstates the potential for practical detection of the modeled saturation changes.

References Alliance (FutureGen Industrial Alliance, Inc.), 2012. Borehole Completion and Characterization Summary Report for the Stratigraphic Well, Morgan County, Illinois. FG-RPT-015 Rev 1. Alliance (FutureGen Industrial Alliance, Inc.), Washington, D.C. Alliance (FutureGen Industrial Alliance, Inc.), 2015. Site Geology and Detailed Stratigraphy of the FGA-1 Stratigraphic Well. FG-02-RPT-0018. Alliance (FutureGen Industrial Alliance, Inc.), Washington, D.C. Arts, R.J., Trani, M., Chadwick, R.A., Eiken, O., Dortland, S., van der meer, L.G.H., 2009. Acoustic and elastic modelling of seismic time-lapse data from the Sleipner CO2 storage operation. Carbon Dioxide Sequestration Geol. Media – State Sci., vol. 59. pp. 391–403 AAPG Stud. Geol. Batzle, M., Wang, Z., 1992. Seismic properties of pore fluids. Geophysics 57, 1396–1408. http://dx.doi.org/10.1190/1.1443207. Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Phys. 12, 155–164. http://dx.doi.org/10.1063/1.1712886. Brie, A., Pampuri, F., Marsala, A., Meazza, O., 1995. Shear sonic interpretation in gasbearing sands. SPE Annu. Tech. Conf. Exhib. 701–710. http://dx.doi.org/10.2118/ 30595-MS. Carcione, J.M., Picotti, S., Gei, D., Rossi, G., 2006. Physics and seismic modeling for monitoring CO2 storage. Pure Appl. Geophys. 163, 175–207. http://dx.doi.org/10. 1007/s00024-005-0002-1. Chabora, E.R., Benson, S.M., 2009. Brine displacement and leakage detection using pressure measurements in aquifers overlying CO2 storage reservoirs. Energy Procedia 2405–2412. http://dx.doi.org/10.1016/j.egypro.2009.01.313. Chadwick, R.A., Arts, R., Eiken, O., Kirby, G.A., Lindeberg, E., Zweigel, P., 2004. 4D seismic imaging of an injected CO2 plume at the Sleipner Field, Central North Sea, vol. 29. Geological Society, London, Memoirs, pp. 311–320. http://dx.doi.org/10. 1144/GSL.MEM.2004.029.01.29. Chadwick, A., Arts, R., Eiken, O., Williamson, P., Williams, G., 2006. Geophysical monitoring of the CO2 plume at Sleipner, North Sea. Adv. Geol. Storage Carbon Dioxide IV. Earth. pp. 303–314. http://dx.doi.org/10.1007/1-4020-4471-2_25. Chadwick, R.A., Noy, D., Arts, R., Eiken, O., 2009. Latest time-lapse seismic data from Sleipner yield new insights into CO2 plume development. Energy Procedia 2103–2110. http://dx.doi.org/10.1016/j.egypro.2009.01.274. Cole, S., Lumley, D., Meadows, M., Tura, A., 2002. Pressure and saturation inversion of 4D seismic data by rock physics forward modeling. SEG Int’l Expo. 72nd Annu. Meet 4–7. http://dx.doi.org/10.1190/1.1817221. Daley, T.M., Solbau, R.D., Ajo-Franklin, J.B., Benson, S.M., 2007. Continuous activesource seismic monitoring of CO2 injection in a brine aquifer. Geophysics 72, A57. http://dx.doi.org/10.1190/1.2754716. Daley, T.M., Freifeld, B.M., Ajo-Franklin, J., Dou, S., Pevzner, R., Shulakova, V., Kashikar, S., Miller, D.E., Goetz, J., Henninges, J., Lueth, S., 2013. Field testing of fiber-optic distributed acoustic sensing (DAS) for subsurface seismic monitoring. Lead. Edge 32, 699–706. http://dx.doi.org/10.1190/tle32060699.1. Dvorkin, J., Nur, A., 1996. Elasticity of high-porosity sandstones: theory for two North Sea data sets. Geophysics 61, 890–893. http://dx.doi.org/10.1190/1.1444059. Ferguson, R.J., Margrave, G.F., 2005. Planned seismic imaging using explicit one-way operators. Geophysics 70, 101–109. http://dx.doi.org/10.1190/1.2073885. Gassmann, F., 1951. Über die elastizität poröser medien. Vierteljahrss-chrift der Naturforschenden Gesellschaft Zurich 96, 1–23. Harbert, W., Lipinski, B., 2010. Technologies Monitor CO2 EOR Floods. The American Oil and Gas Reporter September, p. 119–123. Harbert, W., Daley, T.M., Bromhal, G., Sullivan, C., Huang, L., 2016. Progress in monitoring strategies for risk reduction in geologic CO2 storage. Int. J. Greenh. Gas Control 51, 260–275. http://dx.doi.org/10.1016/j.ijggc.2016.05.007. Hashin, Z., Shtrikman, S., 1963. A variational approach to the theory of the elastic behaviour of multiphase materials. J. Mech. Phys. Solids 11, 127–140. http://dx.doi. org/10.1016/0022-5096(63)90060-7. Hill, R., 1952. The elastic behaviour of a crystalline aggregate. Proc. Phys. Soc. Sect. A 65, 349–354. http://dx.doi.org/10.1088/0370-1298/65/5/307. Ivandic, M., Juhlin, C., Lüth, S., Bergmann, P., Kashubin, A., Sopher, D., Ivanova, A., Baumann, G., Henninges, J., 2015. Geophysical monitoring at the Ketzin pilot site for CO2 storage: new insights into the plume evolution. Int. J. Greenh. Gas Control 32, 90–105. http://dx.doi.org/10.1016/j.ijggc.2014.10.015. Ivanova, A., Kashubin, A., Juhojuntti, N., Kummerow, J., Henninges, J., Juhlin, C., Lüth, S., Ivandic, M., 2012. Monitoring and volumetric estimation of injected CO2 using 4D seismic, petrophysical data, core measurements and well logging: a case study at Ketzin, Germany. Geophys. Prospect. 60, 957–973. http://dx.doi.org/10.1111/j. 1365-2478.2012.01045.x. Jung, Y., Zhou, Q., Birkholzer, J.T., 2013. Early detection of brine and CO2 leakage through abandoned wells using pressure and surface-deformation monitoring data: concept and demonstration. Adv. Water Resour. 62, 555–569. http://dx.doi.org/10. 1016/j.advwatres.2013.06.008. Katahara, K., 1996. Clay mineral elastic properties. 1996 SEG Annu. Meet. pp. 1–4. Kragh, E., Christie, P., 2002. Seismic repeatability, normalized rms, and predictability. Lead. Edge 21, 640–647. http://dx.doi.org/10.1190/1.1497316. Kumar, D., 2006. A tutorial on gassmann fluid substitution: formulation, algorithm and matlab code. Geohorizons 4–12.

5. Conclusions We have developed a workflow for forward modeling of time-lapse seismic data, including constructing seismic velocity models using flow simulation outputs, modeling of synthetic seismic data followed by a basic processing sequence and analysis of processed synthetic seismic data, including determining the detection threshold of changes in timelapse seismic data. We have applied the forward modeling and analysis workflow to the leakage simulations by Williams et al. (2014) for assessing the early leak-detection capabilities of seismic monitoring. The NRMS differences caused by the CO2 leakage are compared to the NRMS detection threshold to evaluate the detectability and the earliest detection time. Our results indicate that high signal-to-noise ratio in the processed seismic data is required to be able to detect the CO2 leakage considered at the model site. Minimizing the background noise in seismic data is crucial for improving the detectability of the seismic method. Increasing the shot density or increasing the dominant frequency of the source wavelet is likely to increase the possibility of the leakage detection and reduce the time needed for the detection. The elastic parameters used in the rock physics modeling have significant effects on the constructed seismic velocity models and the computed synthetic seismic data. Reducing the uncertainties in the elastic parameters could reduce the uncertainties in the predicted seismic responses. The predicted time that the peak amplitude change or the peak NRMS difference occurs can be affected by the fluid mixing model in the rock physics modeling. Disclaimer This paper was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference therein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed therein do not necessarily state or reflect those of the United States Government or any agency thereof. Acknowledgements This technical effort was supported in part by an appointment to the National Energy Technology Laboratory Research Participation Program, sponsored by the U.S. Department of Energy and administered by the Oak Ridge Institute for Science and Education. We would like to thank Dr. Mark Mckoy and Dr. Grant Bromhal at National Energy 51

International Journal of Greenhouse Gas Control 76 (2018) 39–52

Z. Wang et al.

seismic travel time for stress-induced changes. Bull. Seismol. Soc. Am. 97, 281–293. http://dx.doi.org/10.1785/0120060120. Span, R., Wagner, W., 1996. 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. J. Phys. Chem. Ref. Data 25, 1509–1596. http://dx.doi.org/10.1063/1.555991. Souza, R., Lumley, D., 2015. Estimation of reservoir fluid saturation from seismic data: amplitude analysis and impedance inversion as a function of noise. ASEG Ext. Abstr. 2015, 1. http://dx.doi.org/10.1071/ASEG2015ab146. U.S. EPA, 2016. Program class VI well plugging. Post-Inject. Site Care, Site Clos. Guidance 73. Vanorio, T., Prasad, M., Nur, A., 2003. Elastic properties of dry clay mineral aggregates, suspensions and sandstones. Geophys. J. Int. 155, 319–326. http://dx.doi.org/10. 1046/j.1365-246X.2003.02046.x. Vasco, D.W., Daley, T.M., Bakulin, A., 2014. Utilizing the onset of time-lapse changes: a robust basis for reservoir monitoring and characterization. Geophys. J. Int. 197, 542–556. http://dx.doi.org/10.1093/gji/ggt526. Vermeul, V.R., Amonette, J.E., Strickland, C.E., Williams, M.D., Bonneville, A., 2016. An overview of the monitoring program design for the FutureGen 2.0 CO2 storage site. Int. J. Greenh. Gas Control 51, 193–206. http://dx.doi.org/10.1016/j.ijggc.2016.05. 023. Wang, Z., Small, M.J., 2014. A Bayesian approach to CO2 leakage detection at saline sequestration sites using pressure measurements. Int. J. Greenh. Gas Control 30, 188–196. http://dx.doi.org/10.1016/j.ijggc.2014.09.011. Wang, Z., Small, M.J., 2015. Statistical performance of CO2 leakage detection using seismic travel time measurements. Greenh. Gases Sci. Technol. 5, 1–15. White, D., 2013. Seismic characterization and time-lapse imaging during seven years of CO2 flood in the weyburn field, Saskatchewan, Canada. Int. J. Greenh. Gas Control 16. http://dx.doi.org/10.1016/j.ijggc.2013.02.006. Williams, M.D., Vermuel, V.R., Oostrom, M., Porse, S.L., Thorne, P.D., Szecsody, J.E., Horner, J.A., Gilmore, T.J., 2014. Design support of an above cap-rock early detection monitoring system using simulated leakage scenarios at the FutureGen2.0 site. Energy Procedia 4071–4082. http://dx.doi.org/10.1016/j.egypro.2014.11.439. Wood, A.W., 1955. A Textbook of Sound. The MacMillan Co., New York pp.360. Yang, Y.-M., Small, M.J., Junker, B., Bromhal, G.S., Strazisar, B., Wells, A., 2011. Bayesian hierarchical models for soil CO2 flux and leak detection at geologic sequestration sites. Environ. Earth Sci. 64, 787–798. http://dx.doi.org/10.1007/s12665-0110903-5. Zhou, R., Huang, L., Rutledge, J.T., Fehler, M., Daley, T.M., Majer, E.L., 2010. Coda-wave interferometry analysis of time-lapse VSP data for monitoring geological carbon sequestration. Int. J. Greenh. Gas Control 4, 679–686. http://dx.doi.org/10.1016/j. ijggc.2010.01.010. Zhu, X., McMechan, G.A., 1990. Direct estimation of the bulk modulus of the frame in a fluid-saturated elastic medium by Biot theory. SEG Tech. Progr. Expand. Abstr. 787–790.

Lumley, D., 2010. 4D seismic monitoring of CO2 sequestration. Lead. Edge 29, 150–155. http://dx.doi.org/10.1190/1.3304817. Lumley, D.E., Behrens, R.A., Wang, Z., 1997. Assessing the technical risk of a 4-D seismic project. Lead. Edge 16, 1287–1292. http://dx.doi.org/10.1190/1.1437784. Lumley, D., Adams, D.C., Meadows, M., Cole, S., Wright, R., 2003. 4D seismic data processing issues and examples. SEG Technical Program Expanded Abstracts. http:// dx.doi.org/10.1190/1.1817550. Lüth, S., Ivanova, A., Kempka, T., 2015. Conformity assessment of monitoring and simulation of CO2 storage: a case study from the Ketzin pilot site. Int. J. Greenh. Gas Control 42, 329–339. http://dx.doi.org/10.1016/j.ijggc.2015.08.005. Margrave, G.F., 2000. New seismic modelling facilities in Matlab. CREWES Res. Rep. 12, 1–45. Margrave, G.F., 2003. Numerical Methods of Exploration Seismology With Algorithms in MATLAB. Mavko, G., Mukerji, T., Dvorkin, J., 2009. The Rock Physics Handbook: Tools for Seismic Analysis of Porous media. Cambridge University Presshttp://dx.doi.org/10.1017/ CBO9780511626753. Meckel, T.A., Hovorka, S.D., Kalyanaraman, N., 2008. Continuous pressure monitoring for large volume CO2 injections. In: Ninth International Conference on Greenhouse Gas Control Technologies. November 16–20, Washington D.C., U.S.A. Meunier, J., Huguet, F., Meynier, P., 2001. Reservoir monitoring using permanent sources and vertical receiver antennae: the Céré-la-Ronde case study. Lead. Edge 20, 622–629. http://dx.doi.org/10.1190/1.1439008. Mindlin, R.D., 1949. Compliance of elastic bodies in contact. Trans. ASME 71 A-259. Mondol, N.H., Jahren, J., Bjørlykke, K., Brevik, I., 2008. Elastic properties of clay minerals. Lead. Edge 27, 758. http://dx.doi.org/10.1190/1.2944163. Murphy, W.F., 1982. Effects of Microstructure and Pore Fluids on the Acoustic Properties of Granular Sedimentary Materials. Stanford University. Namhata, A., Zhang, L., Dilmore, R.M., Oladyshkin, S., Nakles, D.V., 2017. Modeling changes in pressure due to migration of fluids into the Above zone monitoring interval of a geologic carbon storage site. Int. J. Greenh. Gas Control 56, 30–42. http:// dx.doi.org/10.1016/j.ijggc.2016.11.012. Nur, A.M., Mavko, G., Dvorkin, J., Gal, D., 1998. Critical porosity: a key to relating physical properties to porosity in rocks. Lead. Edge 17, 357–362. Oldenburg, C.M., Lewicki, J.L., Hepple, R.P., 2003. Near-Surface Monitoring Strategies for Geologic Carbon Dioxide Storage Verification. http://dx.doi.org/10.2172/ 840984. Roach, L.A.N., White, D.J., Roberts, B., 2015. Assessment of 4D seismic repeatability and CO2 detection limits using a sparse permanent land array at the aquistore CO2 storage site. Geophysics 80, WA1–WA13. http://dx.doi.org/10.1190/geo2014-0201.1. Schaef, H.T., Strickland, C.E., Jung, K.W., Martin, P.F., Nune, S.K., Loring, J.S., McGrail, B.P., 2017. Injectable contrast agents for enhanced subsurface mapping and monitoring. Energy Procedia 114, 3764–3770. http://dx.doi.org/10.1016/j.egypro.2017. 03.1506. Silver, P.G., Daley, T.M., Niu, F., Majer, E.L., 2007. Active source monitoring of cross-well

52