Research ArticleGEOPHYSICS

Shear-wave anisotropy reveals pore fluid pressure–induced seismicity in the U.S. midcontinent

See allHide authors and affiliations

Science Advances  13 Dec 2017:
Vol. 3, no. 12, e1700443
DOI: 10.1126/sciadv.1700443

Abstract

Seismicity in the U.S. midcontinent has increased by orders of magnitude over the past decade. Spatiotemporal correlations of seismicity to wastewater injection operations have suggested that injection-related pore fluid pressure increases are inducing the earthquakes. We present direct evidence linking earthquake occurrence to pore pressure increase in the U.S. midcontinent through time-lapse shear-wave (S-wave) anisotropy analysis. Since the onset of the observation period in 2010, the orientation of the fast S-wave polarization has flipped from inline with the maximum horizontal stress to inline with the minimum horizontal stress, a change known to be associated with critical pore pressure buildup. The time delay between fast and slow S-wave arrivals exhibits increased variance through time, which is common in critical pore fluid settings. Near-basement borehole fluid pressure measurements indicate pore pressure increase in the region over the earthquake monitoring period.

INTRODUCTION

Seismicity in the midcontinent of the United States has seen a marked increase from the historical average of 21 magnitude (M)≥3 earthquakes per year to 188 M>3 earthquakes in 2011 (1). High levels of seismicity persist to date (2), with 688 M>3 earthquakes in 2014 (3) and uncharacteristically large magnitude events becoming more common. Two of the largest earthquakes in the region occurred in northern Oklahoma in 2016, the Pawnee M5.8 earthquake and the Fairview M5.1 earthquake, with the Pawnee earthquake being the largest in Oklahoma history (4, 5).

Spatiotemporal analysis of earthquakes and wastewater injection point to pore fluid pressure increases as the cause of increased seismicity across the central United States (112), but direct evidence from seismological data has not been documented (6). The time-lapse earthquake shear-wave (S-wave) split analysis presented here demonstrates that pore fluid pressure in the shallow basement has increased over time to a critical pressure and is the cause of the increased seismicity.

S-wave splitting occurs when a wave travels through an anisotropic medium (1318), such as faulted and fractured shallow basement rocks in the U.S. midcontinent (19). Wave propagation oblique to the anisotropy causes the S wave to split into two components: a fast S wave polarized parallel to the fast axis of the anisotropy along the maximum horizontal stress orientation and a slow S wave polarized oblique to the fast axis (1318). The difference in the arrival times between the two S waves is δt, and φ is the azimuthal angle of the fast S-wave polarization orientation (1318).

S-wave splitting analysis of naturally occurring earthquakes has been used to identify critically pore fluid–pressured zones through observation of orthogonal, approximately 90° flips in the fast S-wave orientation, causing it to align with the minimum horizontal stress (13, 14, 20). Flips in φ orientation from maximum to minimum horizontal stress are seen, when the ray path travels a greater distance in elevated pore fluid pressure zones than in normally pressured rock (14). This phenomenon occurs when pore fluid pressure increases sufficiently to cause fluid-saturated fractures that would normally be forced closed by the stress regime to open (14, 15). Change in the travel-path ratio of critical fluid pressure to normal fluid pressure can have significant impact on the range of time delays (δt), with up to an 80% scatter in values (14). The buildup of critical pore fluid pressure allows for natural fault systems to slip (15, 20, 21).

The first observation of induced critical pore fluid pressure change evidenced by S-wave splitting occurred while monitoring the injection of CO2 in a fractured reservoir using time-lapse three-dimensional seismic reflection (22). In this study by Angerer et al. (22), the fast orientation before injection lined up with the maximum horizontal stress; following the injection, which raised reservoir pore fluid pressure by 6.4 MPa (928 psi), the fast S-wave orientation flipped by about 90°, aligning with the minimum horizontal stress.

Here, we examine whether the recent seismicity in northern Oklahoma and southern Kansas exhibits fast S-wave φ flips and increased δt magnitude and scatter, which would constitute direct evidence of critical pore pressure buildup along the ray paths traveled.

RESULTS

The time-lapse earthquake anisotropy analysis spans the period from 2010, when seismicity was sparse in south-central Kansas, to 2016, when earthquake occurrence had become frequent. The region of earthquake activity is constrained to a small area (~10,000 km2) in southern Kansas and northern Oklahoma (Fig. 1) to reduce spatial variability effects and assess temporal changes in seismic anisotropy. Borehole pressure measurements in the Arbuckle Group saline aquifer, approximately 30 m above the basement (~1150 m below ground surface), were obtained from the Kansas Geological Survey (KGS) 1-28 well in the Wellington oil field (Fig. 1).

Fig. 1 Map of the study area in south-central Kansas and northern Oklahoma.

Colored triangles are seismometer station locations, and colored circles are earthquake epicenters. Color identifies the time period of the earthquakes and the corresponding recording stations. Red: 2010–2012, EarthScope Transportable Array (TA); green: 2013–2015, Nanometrics Research Network (NX) and the U.S. Geological Survey (USGS) Networks (GS); blue: 2015–2016 Wellington, Kansas CO2 Sequestration Monitoring network (ZA). Downhole pressure measured in the KGS 1-28 well. Most events used in the study occurred in or near western Sumner County, Kansas. More distant events in northern Oklahoma were incorporated during early time periods, when there was very little seismicity in Kansas. The timing of earthquake occurrence suggests a progression of seismicity from south to north over the 7-year period.

The anisotropy analysis shows a predominant ~90° flip in the φ of events that occurred in 2015–2016 (Fig. 2A3) compared to events from the earlier time windows (2010–2012 and 2013–2015) that exhibit mixed φ responses (Fig. 2, A1 and A2). The early-time solutions of φ (Fig. 2A1) align with the direction of the maximum horizontal stress in the region (70° to 85°), as determined independently from earthquake focal mechanisms and well-bore sonic log data analysis (23). Flipped φ solutions (~330°) are also evident, but the small number of earthquakes that occurred during that period limits the statistical significance of the results in Fig. 2A1. Seismicity increased considerably during the intermediate time period (2013–2015) with solutions of φ exhibiting fast S-wave orientations along the maximum horizontal stress (70° to 85°) and nearly orthogonal to it (Fig. 2A2). In the histogram corresponding to the 2015–2016 period, when most earthquakes occur in Sumner County (Fig. 2A3), the overwhelming majority of fast S-wave φ orientations are offset by approximately 90° from the maximum horizontal stress, causing them to align with the minimum horizontal stress. Although natural stress changes are a possible explanation of the change in fast S-wave azimuth, the similarity in stress orientations from well log data and inversion of recent earthquake focal mechanisms suggests that tectonic stresses over this 7-year period have been stable, as would be expected at an intraplate setting (23). This rotation in φ and the short time frame of its occurrence provide evidence of a change that may be anthropogenic. These changes in φ have been previously identified as an effect of pore fluid pressure increases, where the ray path travels through rock that is critically stressed by pore fluid for a longer distance than rock that is not critically stressed by pore fluid (13, 14, 20, 21). These studies have also identified a large deviation in δt shown to be associated with pore fluid pressure changes.

Fig. 2 Temporal variation in φ and δ.

(A1) Polar histogram of φ from 2010–2012 TA stations (red). Common φ values are near the maximum horizontal stress of 70° to 85° along with flipped values at ~330°. Zero degree values are most often null solutions. (A2) Polar histogram of φ from NX and 2013–2015 GS stations (green) shows the most common solutions in line with maximum horizontal stress and flipped solutions approximately 90° off of maximum horizontal stress. (A3) Polar histogram of φ from 2015–2016 ZA stations (blue) shows the most common solution to be flipped approximately 90° off of the maximum horizontal stress, a direct indicator of critical pore fluid pressure. Polar histograms depict all individual station-earthquake pairs. Black arrows indicate the orientation of maximum horizontal stress at 75°. (B) Average δt/km of earthquakes from 2010 to 2016 showing an increase in variance. Each circle depicts the average δt/km of all stations for an individual earthquake. Yellow stars correspond to average monthly pressure observations in the KGS 1-28 well at the Wellington Oil field. The initial pressure measurement in August 2011 was obtained when the well was drilled. Inset B1 is an expanded view of monthly average downhole pressures from April to November 2016.

The S-wave anisotropy analysis shows an increasing range and scatter in δt estimates (Fig. 2B). The increase in δt suggests increasing anisotropy of the rock, often associated with the fracture density and aperture width (16). It is likely that the basement has become critically stressed by increasing pore fluid pressure. The pore fluid pressure increase reduces the effective stress on the rock, which previously kept fractures that were not parallel to the maximum horizontal stress closed (14). Increasing pore fluid pressure can cause fractures to dilate, increasing the anisotropy and the magnitude of δt. Large scatter in δt estimation is common in critical pore fluid settings because δt is very sensitive to small pressure changes and the length of ray path in critically pressured rock (14, 22). Furthermore, the scatter in δt may be an indicator that the pressure field is nonuniform or “patchy,” with some regions of the shallow basement critically stressed by fluid pressure and other areas not critically stressed (14).

DISCUSSION

The observed flip in φ and the increase in average δt scatter are interpreted as direct evidence of an increase in pore fluid pressure over the time of the investigation. This interpretation is founded on previous studies, where flip in φ and increase in average δt and scatter have been linked to pore fluid increase (1318, 2022, 2426). The seismic anisotropy changes also correlate temporally with changes in downhole pressure data acquired at the KGS 1-28 well in the lower Arbuckle Group saline aquifer at a depth of ~1150 m, approximately 30 m above the basement. Bottom-hole pressure has increased by more than 200 kPa since 2011, when the well was drilled (Fig. 2B). The borehole remained idle until April 2016, when a pressure sensor was installed for continuous monitoring of the lower Arbuckle Group. The high-resolution pressure observations show that downhole pressures have increased at a rate of 3 to 4 kPa per month through the time period of this investigation (Fig. 2B1). The highly permeable Arbuckle Group is believed to be in pressure communication with the underlying fractured basement rocks (6, 7, 23). Increasing pore fluid pressure observed in the Arbuckle Group at KGS 1-28 lends additional support to a regional pore pressure increase in the deeper sedimentary section and shallow basement (depth, 1 to 5 km) resulting from thousands of Arbuckle Group wastewater injection wells throughout southern Kansas and northern Oklahoma (6, 7), although the exact pressure increase at hypocentral depths is still unknown. Previous studies have estimated that changes as small as 0.01 MPa (1.5 psi) (27) to 0.07 MPa (10.5 psi) (6) are sufficient to induce earthquakes, and the region is expected to be very near hydrostatic equilibrium (27, 28).

The TA station (2010–2012) φ results from northern Oklahoma (Figs. 1 and 2A1) are mixed between fast S-wave orientations inline and perpendicular to the maximum horizontal stress. This may indicate that, during that time period, northern Oklahoma and southern Kansas were experiencing the onset of regional elevated pore pressures. We applied the anisotropy analysis to the early-time earthquakes (2010–2013) differentiated by location to two groups: southern group A with events south of 36.5° latitude and northern group B north of 36.5° latitude (Fig. 3). All southern group A events occurred at the earliest time in 2010–2012 and include all the flipped φ solutions (~330°) for the earliest time period and φ along the regional maximum horizontal stress (70° to 85°; Fig. 3A). Flips in φ seen from 2010 to 2012 could signify the onset of pressure increase in the southern portion of the study area because these flips are associated with the southern earthquakes (Fig. 3A). The northern earthquakes of group B exhibit φ solutions that are in line with the maximum horizontal stress and no flipped φ solutions. This observation suggests that the region north of 36.5° latitude during 2010–2013 was not experiencing the same level of elevated pore pressures as the southern region. This spatiotemporal observation of S-wave anisotropy change provides further evidence to a pore pressure increase front moving along Oklahoma and northward to Kansas.

Fig. 3 Spatiotemporal comparison of early time period data.

Separation of early time period earthquakes (2010–2013) between southern group A below the red line at 36.5° latitude and northern group B above the red line at 36.5° latitude. All southern group A earthquakes occurred in 2010–2012 and correspond to the red TA events shown in Figs. 1 and 2. (A) The polar histogram shows the φ solutions of the southern group A earthquakes. (B) The polar histogram shows the φ solutions of the northern group B earthquakes, which include one event from 2010 (red) and the earliest green earthquakes that occurred in 2013. All flipped φ orientations from the red 2010–2012 data shown in Fig. 2 come from the southern group A earthquakes. Corresponding δt values are identified in Fig. 2B to the left of the vertical dashed line at the end of year 2013. The northernmost red event is the earliest earthquake in the data set.

The NX and GS station data covering the time period of 2013–2015 (green events in Figs. 1 to 3) are interpreted as an intermediate stage where pore pressures are increasing, causing a significantly greater number of earthquakes and seismic anisotropy changes becoming more common. This is evident in Fig. 2A2, which shows a mix of φ alignments along the maximum horizontal stress and flipped φ along the minimum horizontal stress. This intermediate time period could be interpreted as a “pressure buildup” period, where some regions are critically stressed, but still, others are not as the pressure front moves across the region. This would produce a mix of φ orientations that are parallel and antiparallel to the maximum horizontal stress. Comparing the early (red) and intermediate (green; Fig. 2, A1 and A2) periods to the φ solutions approximately perpendicular to the direction of maximum horizontal stress in earthquakes from 2015 to 2016 (blue; Fig. 2A3), with ray paths fully contained in southern Kansas (Fig. 1), suggests that the basement has changed to become critically stressed by pore fluid pressure. The predominant flip in φ on the 2015–2016 data suggests that the pore pressure front has reached southern Kansas. This interpretation is in agreement with Arbuckle Group pore pressure changes observed in the KGS 1-28 well.

Although the temporal changes of seismic anisotropy presented in this study are consistent with expected changes induced by pore pressure increase, uncertainties in results may be introduced by the following: (i) Varying earthquake locations: Earthquakes are not colocated, resulting in varying ray paths (azimuth and distance) to the recording stations. Limiting the study to earthquakes that occurred mostly south of stations and in relative proximity (Fig. 1) reduces the azimuthal effect but likely does not eliminate it. Early-time northern Oklahoma earthquake S-wave arrivals (2010–2012 red events) travelled longer paths and may include converted phases. Although converted phases affect δt results, they provide φ information about the rock volume that they pass through. We present the small number of northern Oklahoma earthquakes because they are the only early time period events that occurred near the area of study, recognizing that the overwhelming supporting evidence for this investigation is derived from the large number of recent, proximal earthquakes (2013–2016 green versus blue events). Normalizing δt for distance to milliseconds per kilometer (ms/km) significantly reduces the effect of varying recording distance but adds an additional source of error that affects short ray paths. Proximal earthquakes with a higher ratio of ray path in the shallow sedimentary units could be affected through normalization of δt in ms/km because these shallow units have not experienced fluid pressure change. The similarity in average δt values between the blue and green time periods (Fig. 2B) could be a result of the change in average ray path length. On average, the blue events had shorter ray paths, which are corrected for in calculating δt per kilometer, but these shorter ray paths contain approximately the same length of ray path within the sedimentary units above the Arbuckle Group, as compared to the more distant earthquakes. This normalization affects the longer-distance ray paths significantly less, as seen in the average δt for the red earthquakes, which contain little variance. These ray paths “sample” a larger area, but the effects of the heterogeneity of the pressure field are reduced because of the averaging over a greater distance. (ii) Local site anisotropy: Shallow sedimentary layers (~1 to 1.5 km thick in this region) may introduce uncertainty to anisotropy estimation. However, long ray paths (~10 to 100 km) in the fractured basement likely dominate short travel path anisotropy effects through the horizontally lying, shallow sedimentary layers (29). Furthermore, local site anisotropy effect is expected to remain constant and not contribute to temporal anisotropy variability because fluid pressure changes are confined to the deepest sedimentary unit, the Arbuckle Group (~300 m thick), and the basement below. (iii) Regional anisotropy change is not taken into account and could exist over the study area. Differences in anisotropy from central Oklahoma through northern Kansas could cause δt of the red earthquakes to be lower than the δt of green and blue events. The northernmost red earthquake near 37° latitude (Figs. 2 and 3) is the earliest time event in the entire data set and has the highest average δt of the red time period, indicating that anisotropy could be greater in the northern part of the region. Identifying anisotropy in these data associated with the region without the influence of changing pore fluid pressure would be nearly impossible, given that the data source is earthquakes that are likely caused by a change in pore fluid pressure.

The S-wave splitting analysis presented here is the first direct evidence of increasing pore pressure in the region detected by seismic observations. This increase in pore fluid pressure is the hypothesized cause of the increase in seismicity in the midcontinent (112). Modeling studies have suggested a pressure plume from wastewater injection advancing through central and northern Oklahoma, and southern Kansas is inducing the observed seismicity (4, 7). Earthquake occurrence (Fig. 1) suggests a northward progression of seismicity over time, which is also supported by more than 1600 earthquake observations near the ZA since 2015 (fig. S1). Most of the observed seismicity near Wellington (fig. S1) has occurred in swarms of earthquakes, as noted in other studies of injection induced seismicity (4, 11). A patchy, or heterogeneous, pressure field as opposed to a uniformly expanding pressure pulse may be a more realistic representation of critical pore pressure distribution in the subsurface and its contribution to induced seismicity. Our results show that analyzing the change in anisotropy of the basement is an effective means of identifying critical changes in pore fluid pressure that are the likely cause of fault reactivation and earthquakes in the region (112). However, this data set is limited in its ability to map local regions of pore fluid pressure change and is only capable of identifying a regional change. A more comprehensive data set with better station coverage and a greater number of earthquakes could be used to identify specific areas of pore fluid pressure increase. This methodology could be applied to other regions of potentially induced seismicity to verify that increasing pore fluid pressure related to deep-well injection is the underlying cause of seismicity increases.

MATERIALS AND METHODS

Earthquake data were obtained from networks hosted through the Incorporated Research Institutions for Seismology (IRIS) including EarthScope Transportable Array (TA), which occupied the region in 2010 and 2011, Nanometrics Research Network (NX), the USGS Networks (GS), and the Wellington, Kansas CO2 Sequestration Monitoring network (ZA network) (Fig. 1). A total of 150 earthquakes of M 2 or greater were used for this study. Most events were in the range of M 2 to 3, the largest event was an M 4.3, and earthquakes occurred primarily in the deeper sedimentary section and the shallow basement (depth, 1 to 11 km; fig. S2). The M 2 criterion was chosen so that each earthquake was clearly visible at all monitoring stations with good signal-to-noise ratio. Earlier events used in the study (2010–2012) came from northern Oklahoma, given the scarcity of earthquakes in southern Kansas at that time period.

Earthquake data for 2015–2016 were obtained from the Wellington, Kansas CO2 Sequestration Monitoring network (ZA). Earthquakes were identified from raw data, and magnitudes were calculated by spectral analysis using SEISAN (30). The earthquakes from the EarthScope array as well as the Nanometrics and USGS arrays (2010–2015) were identified through the USGS earthquake catalog. Waveforms were collected from the IRIS data repository, and magnitudes were obtained through the USGS catalog. Earthquake epicenters and station locations are shown in Fig. 1. Histograms of magnitude and depth of the earthquakes used in this study are in fig. S2.

The S-wave splitting analysis was performed using the methods presented by Silver and Chan (31), using the processing technique of Zinke and Zoback (20). A MATLAB code was modified from SplitLab1.0.5 (32, 33). This method performs a grid search for φ and δt, which best removes the anisotropy by calculating eigenvalues that correspond to the covariance matrix of the two orthogonal components (20, 30, 31). Here, the calculation minimizes the second eigenvalue (fig. S3). The φ and δt are calculated from all stations that have a signal-to-noise ratio high enough to visually identify the first-arrival S wave in both waveform plots (for example, in Fig. 4 and fig. S4) and cross-plots (hodograms) of the horizontal channels (figs. S5 and S6). Figure 4 shows an example of raw-waveform S-wave arrivals exhibiting fast and slow components. The corresponding S-wave split analysis results in φ and δt space can be seen in fig. S3.

Fig. 4 Raw seismograms depicting S-wave arrival.

S-wave arrival from WK07 of the Wellington, Kansas CO2 Sequestration Monitoring Network (ZA) with no processing shown. The earthquake occurred 2 March 2016, 15 km from the recording station at a depth of 5.6 km and at M 2.0. Data are shown in north-south (N-S) and east-west (E-W) channels along with the vertical component. Right column figures are expanded views of the left column waveforms between the solid lines depicting S-wave arrival picks. Dashed lines mark the time separation in S-wave arrivals.

All δt and φ values were determined for all stations that had S-wave arrival picks. The δt values were corrected for the distance of the earthquake hypocenter to each recording station, and the average δt (ms/km) of all stations was computed for each earthquake (Fig. 2B). The φ values were plotted as individual values of all stations for each earthquake and mirrored across 180° in the polar plot (Figs. 2, A1 to A3, and 3).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/12/e1700443/DC1

fig. S1. Plot of the Wellington CO2 Sequestration Monitoring network (ZA) earthquake catalog consisting of 1676 events ranging in M from 0.4 to 4.3 and depth from 1 to 11 km.

fig. S2. Depth and magnitude distributions of the 150 earthquakes used in this study.

fig. S3. Plot of the minimization of the second eigenvalue (λ2) in ϕ and δt space from waveforms shown in fig. S4.

fig. S4. Plot of raw channel data from station WK15 of an M 2.7 earthquake that occurred in July 2015.

fig. S5. Hodogram plots of 0.1-s increments corresponding to the 2-s time window identified in fig. S4.

fig. S6. Hodogram plot of S-wave splitting that aligns with the maximum horizontal stress at approximately 75° (marked with red dashed lines).

Additional Acknowledgments

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: The Wellington, Kansas CO2 Sequestration Monitoring network (ZA) seismic instruments were provided by IRIS through the Portable Array Seismic Studies of the Continental Lithosphere Instrument Center at the New Mexico Institute of Mining and Technology. The facilities of the IRIS Consortium are supported by the NSF under Cooperative Agreement EAR-1261681 and the U.S. Department of Energy (DOE) National Nuclear Security Administration. The facilities of IRIS Data Services, IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope Proposal of the NSF under Cooperative Agreement EAR-1261681. Pressure data were provided by J. Victorine (KGS) and P. Simpson (Trilobite Testing). D. Wreath (Berexco) provided access to the Wellington field site. This report was prepared as an account of work sponsored by an agency of the U.S. government. The U.S. government, any agency thereof, or any of its employees do not make any warranty, express or implied, or assume any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed or represent that its use would not infringe privately owned rights. Reference herein 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 U.S. government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. government or any agency thereof. Funding: This research was supported in part by the DOE–National Energy Technology Laboratory under grant no. DEFE0002056. Additional support was provided by the University of Kansas Department of Geology and the Kansas Interdisciplinary Carbonates Consortium Project INS0074816. Author contributions: K.A.N. carried out the bulk of detailed analysis, with continual supervision and guidance by G.P.T. Pressure data and additional scientific input were provided by T.S.B. and W.L.W. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Project data are available through the IRIS Data Management Center. Additional data related to this paper may be requested from the authors. The earthquake catalog may be accessed at http://earthquake.usgs.gov/earthquakes/search/. The IRIS repository may be accessed at http://ds.iris.edu/ds/.
View Abstract

Navigate This Article