Dynamics of geologic CO2 storage and plume motion revealed by seismic coda waves

Significance The sequestration of CO2 into geological formations as a strategy for reducing atmospheric greenhouse gases requires accurate monitoring of CO2 plume for safe long-term storage. However, imaging CO2 plume migration in the subsurface is a challenging problem for existing seismic monitoring methods. Here, we present an approach for CO2 monitoring based on seismic coda waves and apply it to the 3-d period of CO2 injection in a pilot experiment. We find that velocity reduction is nonlinearly correlated with the cumulative CO2 mass inside the reservoir, and that coda waves can effectively monitor the spatiotemporal evolution of CO2 plumes in the subsurface. Our findings suggest approaches for using seismic data in long-term monitoring of subsurface CO2 plume evolution. Quantifying the dynamics of sequestered CO2 plumes is critical for safe long-term storage, providing guidance on plume extent, and detecting stratigraphic seal failure. However, existing seismic monitoring methods based on wave reflection or transmission probe a limited rock volume and their sensitivity decreases as CO2 saturation increases, decreasing their utility in quantitative plume mass estimation. Here we show that seismic scattering coda waves, acquired during continuous borehole monitoring, are able to illuminate details of the CO2 plume during a 74-h CO2 injection experiment at the Frio-II well Dayton, TX. Our study reveals a continuous velocity reduction during the dynamic injection of CO2, a result that augments and dramatically improves upon prior analyses based on P-wave arrival times. We show that velocity reduction is nonlinearly correlated with the injected cumulative CO2 mass and attribute this correlation to the fact that coda waves repeatedly sample the heterogeneous distribution of cumulative CO2 in the reservoir zone. Lastly, because our approach does not depend on P-wave arrival times or require well-constrained wave reflections it can be used with many source–receiver geometries including those external to the reservoir, which reduces the risk introduced by in-reservoir monitoring wells. Our results provide an approach for quantitative CO2 monitoring and plume evolution that increases safety and long-term planning for CO2 injection and storage.

Quantifying the dynamics of sequestered CO 2 plumes is critical for safe long-term storage, providing guidance on plume extent, and detecting stratigraphic seal failure. However, existing seismic monitoring methods based on wave reflection or transmission probe a limited rock volume and their sensitivity decreases as CO 2 saturation increases, decreasing their utility in quantitative plume mass estimation. Here we show that seismic scattering coda waves, acquired during continuous borehole monitoring, are able to illuminate details of the CO 2 plume during a 74-h CO 2 injection experiment at the Frio-II well Dayton, TX. Our study reveals a continuous velocity reduction during the dynamic injection of CO 2 , a result that augments and dramatically improves upon prior analyses based on P-wave arrival times. We show that velocity reduction is nonlinearly correlated with the injected cumulative CO 2 mass and attribute this correlation to the fact that coda waves repeatedly sample the heterogeneous distribution of cumulative CO 2 in the reservoir zone. Lastly, because our approach does not depend on P-wave arrival times or require well-constrained wave reflections it can be used with many source-receiver geometries including those external to the reservoir, which reduces the risk introduced by inreservoir monitoring wells. Our results provide an approach for quantitative CO 2 monitoring and plume evolution that increases safety and long-term planning for CO 2 injection and storage. CO 2 sequestration | geologic monitoring | coda wave | flow rate | mass quantification C onstraints on the spatial distribution and evolution of subsurface CO 2 plumes are critical for the safety and efficiency of large-scale geological carbon storage, which is a promising approach for mitigating anthropogenic climate change. Stored CO 2 perturbs the geophysical equilibria of reservoir rocks by invading fractures and pore spaces and replacing naturally occurring subsurface fluids (e.g., brines), causing spatiotemporal changes in rock properties. The concomitant changes in seismic velocities of the reservoir rocks can provide a basis for detecting and imaging CO 2 plumes (1). However, current approaches rely on active-source, time-lapse seismic techniques based on single-scattering reflected and transmitted waves (2)(3)(4)(5), and a combination of theoretical rock physics predictions (6), laboratory measurements (7,8), and field experiments (9,10) indicate that these can be insensitive to small changes in rock properties depending on CO 2 saturation state and other factors. Correspondingly, the derivation of CO 2 flow dynamics, such as local flow rate, CO 2 saturation, and CO 2 mass, from single-scattering seismic methods remains challenging.
Here, we demonstrate the use of seismic coda waves from continuous, active-source seismic source monitoring (CASSM) to reveal supercritical CO 2 flow dynamics during the Frio-II brine experiment. Coda waves (11) represent scattered and multiply reflected P-and S-waves that traverse the reservoir and sample the spatiotemporal evolution of elastic properties during CO 2 injection. We use a coda wave interferometric technique (12,13) originally developed to detect velocity changes associated with earthquakes (14), local precipitation (15), and mining activities (16). Seismic coda waves bounce repeatedly in the medium, which enables detection of minor variations in physical properties, such as precursory changes in elastic wave velocity before laboratory earthquakes (17,18).
We use the coda wave technique to show that seismic wave velocity in the vicinity of the CO 2 injection decreases systematically with injection rate and cumulated CO 2 mass during the 74-h Frio-II injection. The temporal velocity changes we observe during CO 2 injection are well described by a 3D multiphase flow model based on seismic coda waves. Correspondingly, we ascribe the temporal velocity reductions to the progressive replacement of brine with CO 2 in the injecting reservoir zone and show a robust correlation between injected CO 2 mass and reservoir seismic wave speed. Our findings provide perspectives on CO 2 plume detection, imaging of subsurface flow, and long-term monitoring of CO 2 storage.

Results
Continuous active-source seismic data were collected for about 24 h before and 74 h during the injection of supercritical CO 2 into a high permeability reservoir in the Frio-II pilot site in southeast Texas ( Fig. 1A; see details in Methods). We analyzed data from two subsurface receivers installed at depths of 1,634 and 1,640 m within the Anahuac shale, which is the primary seal above the reservoir sand (Fig. 1B). CO 2 was injected at 1,657 m. Previous results (9, 10) verify that CO 2 was trapped in the reservoir and did not breach the shale (Fig. 1C, Top). Here, we focus on coda waves of the seismic records ( Fig. 1B and SI Appendix, Fig. S1). To improve the signal-to-noise ratio (SNR) of

Significance
The sequestration of CO 2 into geological formations as a strategy for reducing atmospheric greenhouse gases requires accurate monitoring of CO 2 plume for safe long-term storage. However, imaging CO 2 plume migration in the subsurface is a challenging problem for existing seismic monitoring methods. Here, we present an approach for CO 2 monitoring based on seismic coda waves and apply it to the 3-d period of CO 2 injection in a pilot experiment. We find that velocity reduction is nonlinearly correlated with the cumulative CO 2 mass inside the reservoir, and that coda waves can effectively monitor the spatiotemporal evolution of CO 2 plumes in the subsurface. Our findings suggest approaches for using seismic data in longterm monitoring of subsurface CO 2 plume evolution.
individual records we stack the pulse recordings in sets of 3,600 (originally four pulses per second), which leads to a series of full seismic data gathers sampled at ∼15 min intervals during injection. The 15-min stacked records are then processed with a bandpass filter of 200-1,000 Hz before coda wave analysis.
Our data include seismic records obtained with CASSM before and after CO 2 injection ( Fig. 2 and SI Appendix, Fig. S2). Before CO 2 injection, the seismograms show consistent waveforms as a function of elapsed time, including coda waves (Fig.  2C), which confirms that the controlled piezoelectric source signature is highly reproducible. After CO 2 injection, the waveforms of the first arrivals (direct waves) are stable and quite reproducible, but coda waveforms show an evolution in phase and amplitude as a function of injection time. Waveforms at the beginning of injection (black line in Fig. 2 B and C) provide a reference for changes with injection. As shown, first arrivals still correlate well with the reference wave because direct waves traveling from the source to the sensor via the shortest path are not affected by the CO 2 plume. However, the coda waves of even these early seismic traces differ significantly from the reference (Fig. 2C), showing gradual but systematic phase shifts with injection time. Amplitude differences in the later coda ( Fig. 2A), after 70 ms, are also apparent. We interpret these observations as an indication that the coda waves observed above the reservoir are sensitive to CO 2 migration within the heterogeneous reservoir. The time shift fluctuations between coda waves indicate multiple scattering due to different local slowness perturbations in the heterogeneous reservoir. Coda waves passing through local slowness perturbations will reflect the time shifts while others, not passing through, are expected to show small or no time shifts. Next we quantify the coda waves using coda wave interferometry (Methods) following previous methods (12).
Seismic velocity changes are negligible before CO 2 injection, but evolve systematically after injection (Fig. 3). During CO 2 pumping, seismic velocity in the vicinity of the CO 2 injection decreases with time and net injection mass (Fig. 3A). The time derivative of velocity changes d À dv v Á =dt (Fig. 3B) compared with injection rate (Fig. 3C) indicates that a large reduction in elastic wave speed (orange arrows) occurs shortly after CO 2 injection begins (first dashed line in Fig. 3), followed by multiple changes with time. We distinguish three stages of CO 2 injection based on pumping and the reservoir response: stage 1 from 0 to 39 h, a transition stage from 39 to 47 h when pumping was stopped but fluid continued to migrate within the reservoir, and stage 2 when injection occurred between 47 and 74 h.
After injection began we hypothesize that the first significant reduction in seismic wave speed occurred via rapid emplacement of CO 2 in the reservoir ( Fig. 3 A and B). Our data suggest that as the mass of injected CO 2 increased (Fig. 3E), the CO 2 plume migrated up-dip and replaced brine that was originally in the sand reservoir. During this time of increasing CO 2 saturation, our results indicate that elastic wave speed decreased steadily (Fig. 3A). Following a halt in injection due to an unexpected pump failure at about 39 h ( Fig. 3 =dt occurs at about 40 h (orange arrows in Fig. 3B). We attribute this to ongoing CO 2 plume migration within the reservoir. The beginning of the last, and largest magnitude, reduction in velocity coincided approximately with an increase in injection rate (green arrows in Fig. 3 B and C). We hypothesize that faster injection of CO 2 causes the CO 2 plume to spread more widely and that a wider CO 2 plume produces larger coda wave time shifts; we should also note that buoyancy-driven flow is occurring during all phases of injection.
Compared with stage 1 injection, stage 2 shows larger dynamic variations (Fig. 3B). The largest P-wave velocity reduction (1.1 ± 0.23%) (see green arrows in Fig. 3 A and B) occurs approximately 3 h after the CO 2 injection rate was increased to ∼1.0 kg/s (see graydashed line at about 63 h in Fig. 3 A and B). This velocity reduction, and the corresponding time shift of the coda waves, is also clear in the seismograms at about 63 h ( Fig. 2A and SI Appendix, Fig. S2). The results from the top sensor at 1,634 m, shown in SI Appendix, Fig. S3, demonstrate a similar magnitude of velocity reduction in comparison, confirming that the velocity reduction is the result of increasing CO 2 in the reservoir zone. The largest velocity reduction, up to 1.8 ± 0.44% (see purple arrows in SI Appendix, Fig. S3), is seen at the upper sensor, at 1,634 m, consistent with the expectation that scattering was greatest within the CO 2 reservoir, and scattered waves arriving at the topmost sensor will have sampled the reservoir more thoroughly than other waves.
Our data show a correlation between the mass of CO 2 injected and the properties of scattered seismic waves (Fig. 4). As the reservoir CO 2 concentration increases during pumping, the plume causes increased layer reflectivity, intralayer scattering, and intrinsic attenuation, all of which tend to enhance scattering coda waves in the plume layer. Therefore, we attribute the velocity reduction observed with coda waves to the increasing volume of supercritical CO 2 , as CO 2 migrates and replaces brine in the plume layer, reducing the bulk modulus of heterogeneous porous rocks. Pore pressure variations in the Frio-II experiment (as measured by a downhole P/T gauge in the injection well) were relatively small (maximum 0.24 MPa in stage 1) due to the high porosity and permeability of the formation, whereas the overburden pressure was about 16 MPa (Fig. 3D). In stage 2, the pore pressure variations are nearly zero. Such small pore pressure changes are unlikely to cause velocity reduction given the effective stress state of the formations (16). Another observation also supports a connection with injected CO 2 volume. That is, note that the magnitudes of the estimated velocity reductions from two receivers depend nonlinearly on the cumulative CO 2 mass (Fig. 4). The sharp variations in slope are related to rapid injection in stage 2 (blue lines in Fig. 4). This is a unique observation from coda waves. Previous work based on direct P-wave arrivals did not observe such changes (9). We postulate that direct CO 2 monitoring travel time measurements saturate at a fixed value after the plume advances past a specular ray path.
To better interpret our findings, we conducted seismic forward modeling tests to simulate seismic coda wave behavior during CO 2 injection (see details in Methods). We constructed CO 2 seismic velocity models based on predictions of a 3D multiphase flow model for a given CO 2 saturation and calculated seismic velocity changes using modeled coda waves. We overlay the model predictions (black crosses) and the observations of seismic velocity changes in Fig. 3A and SI Appendix, Fig. S3A. They show reasonable agreement between calculated velocity changes from modeled coda waves and estimates from both 1,634-and 1,640-m sensors in the period 0-47 h, which is consistent with the velocity-mass relation (red lines in Fig. 4). During stage 2 injection (after hour 47), seismic velocity reduction is nonlinearly proportional to CO 2 mass, which indicates a more heterogeneous CO 2 distribution. We interpret this as a result of newly injected CO 2 adding to and expanding the original CO 2 plume in the reservoir. Modeled velocity reduction seems not to capture the nonlinearity, particularly in the 1,634-m sensor case ( Fig. 4A and SI Appendix, Fig. S3A). Possible reasons for this include: (i) a constant CO 2 injection rate (0.9 kg/s) used in the model as opposed to dynamic CO 2 injection rates in the field operation, which leads to the strong nonlinear velocity-mass relation (blue lines in Fig. 4), and/or (ii) limitations of the model, which is based on travel time results (19) for which higher CO 2 saturation are not predicted.

Discussion
Our results demonstrate that CO 2 migration can be imaged using seismic coda waves. We demonstrate a relation between cumulative CO 2 mass and elastic wave speed reduction within the 74-h injection and show the unique sensitivity of coda waves to cumulative CO 2 mass. Our interpretation is consistent with previous work (9) and provides additional information on plume dynamics. Previous work based on direct waves found that seismic velocity reduction became negligible, likely due to the fact that CO 2 saturation in the direct ray path was either not increasing or that the direct P-wave response to saturation changes was flat. Spetzler et al. (20) showed that finite-frequency kernels based on single scattering can overcome the limitations of direct ray paths by including Fresnel zone effects. Here, we show that multiple scattering coda waves sample a larger volume of the plume and are more sensitive to the cumulative plume mass compared with local effects.
Coda wave observations of CO 2 migration are significant because traditional approaches for monitoring, verification, and accounting (MVA) seek to provide global estimates of plume mass as well as structural information on CO 2 migration. Our results suggest that seismic coda measurements may enable monitoring in the deep subsurface without the use of wells that directly penetrate the reservoir and for cases where cumulative CO 2 mass is the target parameter. This could have important implications for safe, long-term storage because the use of nonpenetrating monitoring wells decreases leakage risk while preserving the advantages of locating repeatable sources below the highly attenuating near-surface regions. Our observations also suggest that coda wave analysis might have a role in broader areal monitoring if a sufficiently stable source could be utilized. Recently developed orbital sources including the Accurately Controlled Routinely Operated Signal System (21) and the  A key question in seismic monitoring of CO 2 is how to quantify injected CO 2 flow properties at the reservoir scale. Our results show a nonlinear correlation between velocity reduction and the cumulative CO 2 mass (Fig. 4) in contrast with the linear correlation between measured time shifts of reflections and injected CO 2 mass observed in Sleipner time-lapse data (5,23). This is likely due to the dramatic differences in wavelength, timescale, and processing; the Sleipner approach is largely based on classical analysis of reflection pull down and is executed at wavelengths close to two orders of magnitude larger than the Frio-II experiment. In Fig. 4, we show that the relationship between wave speed and CO 2 mass has two key stages: weak nonlinearity in stage 1 and strong nonlinearity in stage 2. This might be related to different injection rates (or cumulative CO 2 mass), which would be consistent with previous experimental results showing that CO 2 saturation and local permeability are strongly controlled by flow rate of the invaded phase (24,25). Therefore, the correlation between velocity reductions and the flow rate presented here provides constraints on quantifying CO 2 flow properties in subsurface processes using coda waves. In addition, different velocity CO 2 mass trends seen in the top two sensors (Fig. 4) confirm the heterogeneous distribution of the CO 2 plume that is very important for fluid migration. Such heterogeneity needs to be well understood and characterized using coda wave tomography (26) with multiple source CASSM data. Analysis of coda waves from multisource CASSM deployments (27) might provide a path toward improved constraints on spatial  heterogeneity, even for locations not directly within the aperture of transmission imaging.
In summary, our results show key connections between the dynamics (flow rate and mass) of injected CO 2 plume and temporal seismic velocity changes using coda waves that may be a promising avenue for the long-term monitoring and quantifying of injected CO 2 plume in the future geological CO 2 storage projects.

Methods
Site Description and Seismic Monitoring Data. The Frio-II pilot was a small-scale injection of supercritical CO 2 into a high permeability reservoir to test geologic storage in saline aquifers (28). CO 2 injection began at ∼7:30 PM, Central Daylight Time, on September 25, 2006. About 380 tons of CO 2 were injected into the blue sand of the Frio formation (Fig. 1A) in about 5 d. Injected CO 2 fluid is in a supercritical (or relatively dense phase) state under the condition of in situ pressure (15 MPa) and temperature (55°C). The pilot site had two wells, about 30 m apart, a down-dip injector and a dedicated up-dip observation well. A highly repeatable piezoelectric source was installed in the injection well near the top of the reservoir sand at about 1,657 m, and multiple hydrophones were installed in the observation well spanning the reservoir (Fig. 1B). A full description of the acquisition design and deployment can be found in Daley et al. (9). We processed about 74 h of continuous cross-well seismic monitoring data, which provided information on the spatial and temporal variation of the CO 2 plume as it migrated updip, driven by buoyancy, across different source-receiver pairs (Fig. 1B). First arrival travel time and amplitude changes, measured by the receivers at various depths in the observation well, allowed hour-by-hour monitoring of the growing CO 2 plume via the induced seismic velocity change (9,19) and seismic attenuation change (10). Each ray path had a unique response with those in the reservoir showing a clear delay in arrival time, which tended to stabilize after a few hours of injection (Fig. 1C). This CASSM waveform data are ideal for constraining the spatiotemporal evolution of the CO 2 plume since it captures variations on the same spatial and temporal scales as a typical reservoir model.
Relative Time Shift dt=t by the Local Similarity Method. The local waveform similarity method (29) provides a smooth continuous measure of similarity between two signals and thus quantitative estimation and extraction of variable time shifts between signals in iterative optimization schemes. In the implementation of iterative optimization, we chose the shaping regularization to enforce smoothness and stabilize the results. To adaptively determine the regularization coefficient (i.e., smoothness of the shaping operator), we started with strong smoothing and iteratively decreased it when the results stopped changing and before they became unstable. The first iteration with a Gaussian smoothing is equivalent to the windowed cross-correlation method. Further iterations using relative amplitude normalization can mitigate amplitude effects on the local similarity. Another parameter (the smoothing radius) needs to be specified in the local similarity method.
In our examples, we used the local waveform similarity module in the opensource Madagascar package (30) to compute the time shift, τ = dt=t, between the reference coda wave trace and other time-lapse traces. A reference seismic record is chosen as the first trace of each selected lapse time window to avoid cycle skipping, which may occur after longer injection times (larger velocity reduction). We chose three lapse time window sizes: 1 (15 min), 5 (75 min), and 10 (150 min) traces. In each window i, the relative time shift τ i for the lapse times are derived and added to the last time shift τ i−1 . After all calculation, three derived relative time shifts for three window sizes are averaged (SI Appendix, Fig. S4). We observed that the trend of time shifts is consistent for three window cases. We used three independent estimates of the relative time shift to calculate mean and SD. Assume that the time shift dt=t is caused by a spatially homogeneous relative velocity change dv=v, the relative time shift is therefore the reflection of the relative velocity variations dv=v (12).
Modeling of Seismic Coda Waves During CO 2 Plume Migration in Frio-II Site. The CO 2 flow model used in this study was developed in Daley et al. (19) using the TOUGH2 (Transport of Unsaturated Groundwater and Heat 2) flow simulator (31) for simulating CO 2 flow dynamics in porous and fractured media. TOUGH2 uses the hysteretic formulation for capillary pressure and relative permeability, and an equation of state package (32) to treat a twophase (liquid, gas), three-component (water, salt, CO 2 ) system in pressure/ temperature regimes above the critical point of CO 2 (P = 7.38 MPa, T = 31°C). In rock physics, forward simulations, brine, and supercritical CO 2 properties are calculated for in situ pressures (15 MPa) and temperatures (55°C) (V p values of 1,574 and 333 m/s, respectively, with corresponding densities of 992 and 653 kg/m 3 ) using the Batzle-Wang model (33) for brine properties and the National Institute of Standards and Technology empirical database (34) for CO 2 properties. Because a minimal temperature change is observed in field measurements, the temperature is assumed to remain constant for the simulations. The numerical schemes use an integral finite difference method for space discretization and first-order fully implicit time differencing. The simulation finally generated a time series of simulated saturation distributions for 5 d of CO 2 injection at an average rate of 0.9 kg/s (76 T/d). These time-lapse CO 2 saturation models were further improved iteratively by matching the CASSM field data (19).
Next, we employ the patchy saturation substitution approach (35,36) to transform changes in CO 2 saturation from the TOUGH2 simulation to changes in observable seismic properties (P-wave velocity, density, and attenuation) using the models mentioned previously. The heterogeneity in the seismic models is introduced by a random perturbation (<5% percent velocity perturbation) by von Kármán correlation function (37). SI Appendix, Fig. S5 shows the snapshots of the updated CO 2 seismic velocity models. We finally employ a time-domain finite-difference acoustic wave equation solver (38) to model seismic data for each CO 2 velocity model with the Frio-II source-receiver geometry. SI Appendix, arrivals and coda waves before and after CO 2 injection are consistent with our field observations (Fig. 2).
Data Availability. Seismic field data used in this study were acquired by Lawrence Berkeley National Lab (LBNL) with support from the US Department of Energy. All data can be freely accessed on the Energy Data eXchange (EDX), operated by the National Energy Technology Laboratory (9) and utilized under the terms of the Creative Commons Attribution (CCA) license.