Oklahoma's induced seismicity strongly linked to wastewater injection depth

See allHide authors and affiliations

Science  16 Mar 2018:
Vol. 359, Issue 6381, pp. 1251-1255
DOI: 10.1126/science.aap7911

Injection depth matters for induced earthquakes

Wastewater injection has induced earthquakes in Oklahoma, but the relative importance of operational and geologic parameters in triggering such earthquakes is unclear. Hincks et al. developed an advanced Bayesian network to determine the interplay between these parameters in Oklahoma. The injection depth above the crystalline basement was the most important parameter when considering the potential for release of seismic energy. This modeling strategy may provide a way to improve forecasts of the impact of proposed regulatory changes on induced seismicity.

Science, this issue p. 1251


The sharp rise in Oklahoma seismicity since 2009 is due to wastewater injection. The role of injection depth is an open, complex issue, yet critical for hazard assessment and regulation. We developed an advanced Bayesian network to model joint conditional dependencies between spatial, operational, and seismicity parameters. We found that injection depth relative to crystalline basement most strongly correlates with seismic moment release. The joint effects of depth and volume are critical, as injection rate becomes more influential near the basement interface. Restricting injection depths to 200 to 500 meters above basement could reduce annual seismic moment release by a factor of 1.4 to 2.8. Our approach enables identification of subregions where targeted regulation may mitigate effects of induced earthquakes, aiding operators and regulators in wastewater disposal regions.

Oklahoma’s ∼900-fold increase in the annual rate of seismicity since 2009 (1) (fig. S1) makes it the most seismically active region in the contiguous United States. Wastewater disposal via deep injection wells explains the surge in earthquake activity (17). In many cases, pressurized fluids migrate from target formations downward into the crystalline basement (1, 710), where most earthquakes occur (1118). The disposed fluids may locally increase pore fluid pressure and reduce the effective stress along faults (2, 9, 11), resulting in induced seismicity (2, 19). In September 2016, a magnitude (Mw) 5.8 earthquake centered near Pawnee caused injury and damage to buildings, prompting the immediate suspension of 37 disposal wells (20). The number of damaging earthquakes (5, 6, 21) resulted in litigation against well operators (22), increased risk to critical infrastructure (23), and considerable financial loss in the reinsurance industry (24).

The escalating seismicity compelled Oklahoma regulators to introduce emergency directives aimed at managing injection (25), informed by an empirical understanding of how seismicity responds to changes in injection rates or volumes (2, 4, 6). These measures appear to have decreased the earthquake count (Mw ≥ 3) in the period 2015 (21) to the present (fig. S1). Considering the purported strong association between seismicity and oil production in Oklahoma (26), a major driver of this trend could have been abated injection due to relatively low oil prices. However, a case has been made (27) that the probabilities of occurrence of larger-magnitude earthquakes are not decreasing as rapidly as other models suggest (1). The total seismic moment release has declined only modestly despite decreased injection volumes (fig. S1). The ultimate impact of reduced injection may take years to materialize given observed lagging in the response of seismicity levels to injection (11, 27, 28).

The Oklahoma Corporation Commission (OCC) requires operators to prove that their injection does not encounter basement lithologies (29), but debate remains on the precise influence of injection depth. Although a study encompassing the wider Central and Eastern United States (CEUS) found no significant correlation between proximity to basement and earthquake occurrence (4), other studies suggest that injections near, or within, crystalline basement might increase the likelihood of seismicity locally (12, 30, 31). Therefore, we need to analyze the joint effects of spatial, operational, geologic, and seismicity parameters (4, 19, 21).

We developed a Bayesian network (BN) to quantitatively evaluate correlations between well operational parameters, geologic setting, and seismicity, focusing on the “Oklahoma induced seismicity zone” (5) (Fig. 1). Our BN framework capitalizes on mathematical advances (32) and innovative software, UNINET (33), implemented in R (34). UNINET models joint dependency using continuous data distributions, bypassing the limiting assumptions of discrete approximations (32). Monthly underground injection control (UIC) well data for January 2011 to September 2016 were obtained from OCC (35), including injected volume, pressure, duration of well activity, and total well depth. Earthquake location and magnitude data, from the Advanced National Seismic Systems’ (ANSS) earthquake catalog, ComCat (36), were used to calculate total seismic moment release. We considered earthquakes of Mw ≥ 2, with a hypocentral depth <10 km (fig. S3), and a hypocenter falling within 20 km of an injection well (32). We chose this distance on the basis of the documented occurrence of associated earthquake swarms within 10 km (37) to 20 km (3) of high-rate disposal wells, supported by additional calculations that indicate strong correlations between injection and moment up to ∼20 km (32). Basement depths at well sites were estimated by spline interpolation (32, 38) using 1232 depth records from Oklahoma Geological Survey (39) (fig. S5).

Fig. 1 Wastewater injection and seismicity in Oklahoma.

(A). Map showing the surface expression of key geologic provinces of Oklahoma (41) and class II UIC wells (35) (related to oil and gas production). Black polygon defines the “Oklahoma induced seismicity zone” delineated by the USGS (5). (B) Map showing total injected volume (bbl) for 2011 to 2015 (30-km Gaussian kernel density). Symbols denote locations of earthquakes from 2011 to 2016, including the 5 November 2011 Prague (7, 10, 11) and 3 September 2016 Pawnee (18, 21, 28, 44) events. (C) Map showing total seismic moment release from 2011 to 2016 (20-km Gaussian kernel density), with mapped faults in the sedimentary cover including the Nemaha uplift (arrow) also shown (42) [this does not include all large-scale basement structures (17)]. Yearly injection volumes and seismic moment release are shown in fig. S2.

Over the study period, >10,000 active class II UIC wells in Oklahoma (Fig. 1, A and B) collectively injected on average ∼2.3 billion barrels (bbl) of fluid per year via saltwater disposal (SWD) and enhanced oil recovery (EOR) (35) (figs. S1 and S2 and table S1). We considered class II injections into the full range of target geologic formations, to capture all potential contributions to induced seismicity (32). Wastewater disposal is particularly intensive within the Cambrian–Ordovician Arbuckle Group, a sedimentologically and structurally heterogeneous sequence of dolomitized carbonates and breccias (40) ∼0.15 to 2 km thick throughout most of the study area. Sedimentary cover is hydraulically connected to Precambrian crystalline basement (1, 79) (fig. S6) via an extensive fault network (14, 17, 41, 42) (Fig. 1C).

High-rate injection has been linked to increased seismicity (4). However, analysis using BN and regression models indicate that other parameters have an important, but largely unquantified, impact on seismic hazard. We explored correlations between key operational and spatial characteristics and seismicity, comparing a number of network configurations with a simple linear regression (Fig. 2) to determine the most appropriate model (32). We defined nodes for key observables, including well operating parameters (depth, monthly, annual and cumulative injected volume), geologic parameters (distance to basement), and total seismic moment release within 20 km of the well in the year ahead (Fig. 2, A and B) (32)—a standard time scale for operational and seismic hazard forecasts (5, 6, 21). To improve forecasting for deeper and higher-volume wells associated with higher seismicity (Fig. 3 and fig. S7) (32), our models focused on wells ≥1 km deep and cases of monthly injection >10,000 bbl, i.e., those of greater concern to regulators (20). We used UNINET in “data mining” mode to estimate the set of joint normal copulae (multivariate probability distributions) that most closely represented the input data (32). We used the resulting models to estimate moment release for any given set of operating conditions.

Fig. 2 BN structure, unconditional rank correlation coefficients, and cumulative distribution functions (CDFs) for Test 2.

(A) The saturated network and (B) the unsaturated network. Arcs denote influence between nodes. (C) Unconditional rank correlation coefficients for all nodes with log10 annual moment release (empirical data and BNs). This correlation is a measure of the strength of influence of each observable on the node of interest (moment release). The BN uses joint normal copulae to represent the empirical data—the closer the values in the BN rank correlation matrix are to the empirical rank correlation matrix, the better the model approximates the data (matrix determinants can also be compared). Full correlation matrices are provided in table S5. (D and E) CDFs of log10 annual moment release for the BN forecast (mean and median estimate), linear regression model, and empirical data for (D) saturated BN and (E) unsaturated BN. This function shows how closely the model output distribution matches the observations. Results are shown for Test 2–learning using 90% of the observations, testing with the remaining 10% (see fig. S11).

Fig. 3 Probability distributions for well depth and distance to basement using BN inputs for January 2011 to December 2015 (wells ≥1 km deep and injection records >10,000 bbl/month).

Probability density for (A) well depth below surface; (B) depth relative to basement and (C) distance to basement; identifying (i) cases with low or no seismic activity (log10 moment in 1 year ≤13; shaded blue) and (ii) cases with relatively high seismic activity (log10 moment in 1 year ≥15; shaded pink). The distributions in (B) and (C) show that wells drilled closer to the basement are associated with higher seismic moment release. Distributions for injection volume are shown in fig. S7.

The saturated BN (Fig. 2A) (32) gave the most complete representation of the observed data under stationary conditions, outperforming the regression model (in-sample testing, table S3 and figs. S9 and S10). However, the unsaturated BN (Fig. 2B) provided better forecasting under the changing injection regime of 2015 to 2016 (fig. S1B). This is because rank correlations between volume variables in the saturated BN preserved relationships (associated with increasing injection) that were no longer valid under post–2015 operating conditions.

We used two test data sets for model verification (fig. S11 and table S3). The first demonstrated operational use of the BN to forecast future moment release, learning with observations from January 2011 to June 2015, and forecasting from July to December 2015. The second test used 90% of the full data set for learning (January 2011 to December 2015, randomly sampled) and the remaining 10% for verification (fig. S11) (32). We compared forecast results using the BN and regression models (Fig. 2, figs. S12 and S13, and table S6) and found that in both cases, the BN gave the best performance.

The spatial distribution of earthquakes (Fig. 1, B and C, and fig. S2) implies important underlying geologic characteristics influencing seismicity. The BN cannot explicitly represent these characteristics because of the variability and paucity of geologic data at relevant depths and resolution across the region. To address this issue, we applied the models iteratively to identify where (spatially) they systematically under or overforecast, because of latent features outside the scope of the model. We used both BN and regression models to estimate seismic moment for a subset of the training data and used the forecast error to compute a “geospatial correction” map—a proxy for unobserved geologic data (32) (fig. S8). We repeated the learning step with the training data using the geospatial correction as an additional input node. We used this revised network (Fig. 2) for final forecast verification (figs. S10, S12, and S13) (32).

The geospatial correction shows the strongest (negative) correlation (–0.55) with annual moment release (Fig. 2C), supporting the observation that lithologic (9) and fault network characteristics (14, 1618) collectively play a critical role in determining susceptibility to induced seismicity. Root mean squared errors, mean absolute errors, and median average deviation for the BN and regression models, both with and without the geospatial correction node, show the improvement in forecast skill resulting from this additional step (table S6, Test 2). This degree of improvement suggests that the geospatial correction effectively characterizes the impact of important, but unobservable (latent), factors. It could additionally prove valuable in mapping spatial variation in induced earthquake hazard (21).

Our analysis with UNINET shows that distance to the basement interface is more strongly correlated (–0.34) with annual seismic moment release than injected volume; and cumulative injected volume has a stronger correlation (+0.26) with annual seismic moment release than annual (+0.20) or monthly (+0.18) volume (Fig. 2C, empirical rank correlations). However, we expect these volume correlations to change with time owing to the effects of regulation, market-driven changes in oil production and associated disposal, and time lags between injection and induced seismicity. For example, cumulative volume will become less well correlated with annual moment if the current trend of reduced seismicity continues (fig. S1A).

UNINET uses joint normal copulae (32,33). As several input distributions are highly skewed and seismic moment is censored (because of catalog incompleteness for Mw < 2.5; fig. S4), the normal copula cannot perfectly represent the data. We see this in the difference between the empirical and BN unconditional rank correlations (table S5) and cumulative distribution functions (CDFs; Fig. 2 and figs. S12 and S13). However, calculations suggest that using optimal copulae would yield only a small improvement, negated by greatly increased computational complexity (32). The asymmetry of the data also means that a linear model will not perform as well as the BN (e.g., Fig. 2, D and E).

We demonstrate how the BN could be used to quantify the expected impact of UIC regulations with three applied examples. We can (i) reduce injection depth by a fixed amount [to simulate well plug back (29) or optimal depth of a new well], (ii) cap injection volume according to imposed regulatory limits (25, 43), and (iii) restrict injection depth in the Oklahoma induced seismicity zone, by enforcing a minimum permissible distance above the basement interface. We intend these examples as demonstrations of how the BN could be used in practice. Before using this approach to inform operational or regulatory practices, further detailed simulations would be required to account for local-scale susceptibility to induced seismicity.

For a subset (10%) of well records (fig. S11), we reduced well depth by 1000 m, and the resulting “test data” were used to produce a set of revised forecasts for annual moment release. We calculated the change in mean moment for each individual case due to the change in depth (Fig. 4). For above-basement wells, raising the well by 1000 m shifted injection further from the basement (shallower) and led to a reduction in seismic moment release. We found the greatest reduction in wells that inject higher volumes (Fig. 4A). Conversely, increasing depth to inject closer to (or into) the basement led to an increase in seismic moment release (fig. S14).

Fig. 4 Simulated impact of raising the injection well level or capping monthly injection volume.

(A) For a subset of records, well depth was reduced by 1 km and used to generate revised estimates for annual moment release using the BN (colored points) and linear regression models (black points). The plot shows change in annual moment versus original well depth and predicts the greatest reduction in moment release for wells initially closest to the basement interface. The joint effect of cumulative volume is also seen—wells with higher cumulative injections are expected to experience proportionally greater moment reduction. (B) Predicted change in mean annual moment resulting from a regulatory limit of 15,000 bbl/day (450,000 bbl/month) versus monthly injection with no cap. This plot shows cases of high injection only (where the limit is applied). Points are colored by total reduction in annual volume to illustrate the importance of maintaining injection limits over time. Inset: CDFs for moment release (BN mean estimates) show the overall effect on all high-volume wells.

To model the potential impact of regulating volume, we applied a daily limit of 15,000 bbl per well, as imposed by Oklahoma regulators (25), similar to the Kansas 16,000 bbl/day limit for areas of “seismic concern” (43). We applied this limit to the full data set and used the revised monthly injections to recalculate annual and cumulative volumes for each well (32). We then used the BN to simulate the revised annual moment release for all individual monthly well records with reduced injection (Fig. 4B).

As expected, our simulations showed that wells subjected to the largest drop in monthly and annual injected volume yielded the greatest reduction in moment release (Fig. 4B). We found a small number of exceptions in which large changes in mean annual moment were predicted for wells with relatively low annual injection rates (Fig. 4B). These cases largely corresponded to wells injecting close to the basement (fig. S15), where more modest changes in injection volume can yield proportionally larger decreases in moment (Fig. 4). This demonstrates the importance of quantifying the joint effect of operational parameters (depth, cumulative volume and injection rates) on future seismicity levels.

We performed a final set of simulations to quantify the statewide impact of depth restriction, using depth limits of (i) 200 m and (ii) 500 m above basement, wherever this could be practically applied. Using well records for December 2015 to predict annual moment release in the year ahead, we found an indicative (mean) reduction in total annual moment release for the Oklahoma induced seismicity zone of a factor of ∼1.4 for the 200-m limit and a factor of ∼2.8 for the 500-m limit (32).

A distinct advantage of the BN is the ability to update the model and produce revised forecasts as data become available. For example, using newly released (September 2017) OCC injection data (35), we tested BN forecasting for 2016, first using the model as developed above (fig. S16) and second, using an updated geospatial correction and moving window for learning, which (by partially accounting for temporal changes in the system) improves performance (figs. S17 and S18) (32). The empirical and unconditional BN rank correlations (table S7) computed with the complete updated data set are not considerably different from those for 2011 to 2015 (table S5), lending credence to the basic model.

Our results indicate that seismic moment release in Oklahoma is strongly correlated with injection proximity to crystalline basement (Figs. 3 and 4A and fig. S7), corroborating previous hypotheses (30, 31) but contrary to observations encompassing the CEUS (4). We attribute this disparity to scale dependence (19). In Oklahoma, the permeability structure of the Arbuckle Group permits downward fluid migration into crystalline basement (1, 9), causing reactivation of optimally oriented strike-slip faults (14, 1618). On a continental scale, lithologic and structural heterogeneity and concomitant stress-field variations could obscure any association between depth and seismicity, as some lithologies are more susceptible to earthquakes (9) (fig. S6); e.g., where fracture flow causes localized pore fluid pressure increases, allowing rupture along critically stressed faults (9). Regions of relatively high-permeability basement in other parts of the CEUS will conceivably permit greater pore fluid diffusion and pressure release, inhibiting earthquakes (9). Induced seismicity is also controlled by local hydrogeologic conditions (30) and the location, scale, and orientation of critically stressed faults—all unique aspects of the geotectonic setting. This highlights the importance of local-scale assessments.

Our Oklahoma assessment exploited a 6-year record of fluid injection to develop a comprehensive understanding of the controls on induced seismicity. The two-stage BN quantifies the joint effects of operational parameters and latent spatial features on seismic moment release, facilitates regular model updating, and offers improved forecast performance compared to traditional linear models. This approach has important implications for statewide assessments and regulatory decision-making under geologic uncertainty [note that the Pawnee (16, 18, 44), Cushing (44), Prague (11, 14), and Fairview (45) earthquakes occurred on previously unrecognized faults or splays]. Our framework allows regulatory actions to be evaluated on a rational, quantitative basis in terms of seismic effects. The model could be used to identify evolving vulnerability, complementing zoned seismic hazard forecasts of the U.S. Geological Survey (USGS) (21), and provide a basis for targeted management of induced seismicity in Oklahoma and other wastewater disposal regions. We welcome the support of researchers, operators, regulators, and policy-makers in updating and extending our model.

Supplementary Materials

Materials and Methods

Figs. S1 to S18

Tables S1 to S7

Data File S1

References (4654)

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: T.H. and W.A. were supported in part by the CREDIBLE consortium (UK Natural Environment Research Council grant NE/J 017299/1). W.A. received support in part from Risk Management Solutions (RMS) for a conceptual study associated with induced seismicity attribution. D. Ababei at LightTwist Software provided invaluable assistance with UNINET software and scripting (available at The study used UIC injection well data obtained from the Oklahoma Corporation Commission ( and earthquake data from the USGS ANSS Comprehensive Catalog ( Input data are provided in the supplementary materials. We acknowledge helpful discussions with R. Muir-Wood (RMS, London) and R. J. Briggs (formerly Praedicat Inc., CA). W.A. and T.H. conceived the research. T.H. performed the modeling and analyzed and refined the BN. W.A. provided seismological interpretation. T.G. performed geologic interpretation and coordinated and wrote the paper. R.C. provided support with statistical analysis and performed additional UNINET tests. All authors contributed to the interpretation, discussion, and writing.

Stay Connected to Science

Navigate This Article