Comment on "How will induced seismicity in Oklahoma respond to decreased saltwater injection rates?" by C. Langenbruch and M. D. Zoback.

The state of Oklahoma has experienced an unprecedented increase in earthquake activity since 2009, likely driven by large-scale wastewater injection operations. Statewide injection rates peaked in early 2015 and steadily decreased thereafter, approximately coinciding with collapsing oil prices and regulatory action. If seismic activity is primarily driven by fluid injection, a noticeable seismogenic response to the decrease in injection rates is expected. Langenbruch and Zoback suggest that "the probability of potentially damaging larger events, should significantly decrease by the end of 2016 and approach historic levels within a few years." We agree that the rate of small earthquakes has decreased toward the second half of 2016. However, their specific predictions about seismic hazard require reexamination. We test the influence of the model parameters of Langenbruch and Zoback based on fits to observed seismicity distributions. The results suggest that a range of realistic aftershock decay rates and b values can lead to an increase in moderate earthquake probabilities from 37 to 80% in 2017 without any further alteration to the model. In addition, the observation that all four M ≥ 5 earthquakes to date occurred when injection rates were below the triggering threshold of Langenbruch and Zoback challenges the applicability of the model for the most societally significant events.

In their recent work, Langenbruch and Zoback () (hereafter referred to as L&Z) quantify the seismogenic response to a reduction in disposal rates based on an adapted empirical model, originally developed for isolated injection, for example, in geothermal or hydrocarbon reservoirs (). The adapted model assumes delayed earthquake triggering, a specific triggering threshold in addition to Gutenberg-Richter magnitude scaling, Omori-like rate decay, and a seismogenic index related to fault density and stress state. Whereas spatial variations in earthquake response to injection are governed by the seismogenic index, temporal seismicity rate changes are thought to be controlled by injection volume and the p value of the Omori-like rate decay after peak injection (, ). All the modeling parameters in L&Z are determined empirically by fitting the observed data and are then applied to predict seismicity rates and magnitude distributions.

There are two primary factors that directly affect the predicted probability of moderate earthquakes in the statistical model, that is, the temporal seismicity decay after peak injection and the b value used to describe magnitude distributions. In this comment, we examine the observed range of p and b values and report implications for exceedance probabilities of moderate earthquakes in 2017 and 2025. This exercise shows that the probability of having a significant earthquake could be underestimated by a factor of 2, even within the confines of the model assumptions. We additionally show that all four M ≥ 5 earthquakes occurred at times when injection rates were below the triggering threshold in their model, suggesting that the model structure may also require reexamination.

Here, we use the same data sources as L&Z and include the slightly more recent seismic data now available. These data include monthly injection rates in all Arbuckle wells from the Oklahoma Corporation Commission and seismicity records from the Advanced National Seismic System earthquake catalog between 2009 and 2017. The study region is constrained to the two areas in central and northwestern Oklahoma that have been subject to injection rate reductions directed by the Oklahoma Corporation Commission. The earthquake record, on average, is complete to magnitude 2.5, with some fluctuations above this value, so that a magnitude of completeness of Mc = 3 is selected in agreement with L&Z. Our analysis, similar to observations by L&Z, found that the monthly rate of M ≥ 3 events is lower in 2016 than in 2015, approaching values close to 2014.

We first address the Omori-like seismicity rate decrease after peak injection. L&Z use a modified Omori relationship of the form R(t) = R0/(t/t0) to predict future seismicity rates, where R0 is the seismicity rate at the start of significant injection rate decrease in 2015, t0 is the injection time window above the triggering threshold, and p is the rate decay exponent. Their value of p = 2 is taken from an earlier study of isolated injection operations in geothermal reservoirs (). In the supplement, L&Z even suggest that p = 2 is a conservative estimate, whereas the vast majority of studies find tectonic p values close to 1 ().

Although there are legitimate questions about the applicability of a single decay law for statewide seismicity, we proceed here taking the L&Z model at face value. We determine p from a simple nonlinear least-squares fit using the start time of aftershock-like decay in the study by L&Z in June 2015 until the end of 2016, with R0 = 80 and t0 = 2, and find that p = 1.4 ± 0.16 (Fig. 1A). Although L&Z do not specifically fit any catalog to determine p, the comparison in their Fig. 3 implies that the Omori-like decay is justified by observed rate changes in a catalog that was partially declustered for aftershocks associated with M ≥ 4.7 mainshocks. We follow a similar strategy but consistently remove aftershocks and secondary aftershocks for all events above Mc within magnitude-dependent space-time windows (). The resulting rate of seismicity decay is substantially reduced, that is, seismicity background rates decay gradually over time, supporting a best-fitting value of p = 1.2 ± 0.08 (R0 = 22, t0 = 2). In either case, the resulting values of p are significantly below the assumed value of 2 by L&Z.

Fig. 1

Annual exceedance probability as a function of p and b values in 2017 and 2025.

For 2017, the probability of a M ≥ 5 event varies between 37 to 80% and remains significantly above the tectonic probabilities. Probabilities remain far above historic (that is, tectonic) values even in 2025. (A) Modified Omori aftershock decay following L&Z for the complete catalog (black) and a catalog with aftershocks removed (gray). The best-fitting p values are 1.4 (yellow) and 1.2 (red), respectively, which is significantly smaller than p = 2.0 (green) in the study of L&Z. (B and C) Magnitude distributions (colored circles) for different time periods in the study of L&Z (see legend) and maximum likelihood b value fits (colored lines) with Mc = 3.0 (triangles) for western and central Oklahoma. The best-fitting values of b vary between 1.14 and 1.46 (see legend). (D) Annual exceedance probability for tectonic events before 2009 (gray) and predicted behavior based on different values of p and b in 2017 (orange) and 2025 (blue). The values used in the study of L&Z are shown as solid blue and orange curves.

Next, the observed magnitude distributions within the time windows analyzed by L&Z are examined to determine b values from a maximum likelihood fit (Fig. 1B). We find that b values vary significantly between ~1.14 and ~1.46. In the following, we consider these values as a plausible range within the statistical model, whereas L&Z restricted their attention to b values of 1.33 and 1.41. The consequences of both b value variations and a more subtle rate decay can be quantified by recalculating the exceedance probability of events with moderate magnitudes. These probabilities were derived by L&Z from a cumulative Poissonian distribution of the formwhere Pe is the probability of observing an event with M ≥ Mi, λ is the rate parameter (that is, observed rate, modeled rate based on a seismogenic index, or predicted rate based on the modified Omori relationship), and b is the exponent that describes the magnitude distribution above the magnitude of completeness Mc. For estimates of future earthquake rates, λ depends mainly on the p value of the modeled seismicity decay rate. Lower p values lead to prolonged periods of high seismicity rates, and lower b values increase the relative number of large magnitude events (Fig. 1C). As a consequence, the probability of exceeding a magnitude 5 event in 2017 may be as high as 80%. This is more than twice the value of 37% inferred from the parameters in the study of L&Z. It should be clear from this discussion that the exact values of p and b are debatable, and therefore, the entire range of probabilities in Fig. 1C should be viewed as plausible outcomes of the model.

Although we have successfully reproduced the quantitative results of the statistical model in Fig. 1C, we are perplexed by the qualitative statement by L&Z that “the probability of potentially damaging larger events, should significantly decrease by the end of 2016 and approach historic levels within a few years.” Both the results by L&Z and the reproduced Fig. 1C do not seem to support this assessment because none of the corresponding parameter choices lead to 2025 probabilities as low as historic values (that is, before 2009).

Up to this point, the discussion has been focused on assessing parameter sensitivity within the confines of the L&Z statistical model. We now explore observations that raise issues about the model framework. Equation 1 assumes that the occurrence of large-magnitude events can be determined by extrapolating the rate of small-magnitude events, even if observational and predictive time windows are short. We test this assumption by comparing trends in cumulative seismic moments and cumulative earthquake number between 2009 and 2017 (Fig. 2A). The latter is dominated by small-magnitude events, whereas cumulative moment release is mainly controlled by large-magnitude events. Seismicity rates seem to decrease after the Fairview earthquakes in February 2016, whereas seismic moment release continued to accelerate until the end of 2016 partially due to the M = 5.8 Pawnee and M = 5.0 Cushing earthquakes. This increase in seismic moment is in apparent violation of the prediction in L&Z [see also Yeck et al. ()]. Moreover, the Prague M = 5.7 event in 2011 occurred at a time when the L&Z model predicted a very low probability of an M ≥ 5 event, and all four M ≥ 5 events (Prague, Fairview, Pawnee, and Cushing) occurred when injection rates were below the triggering threshold in the study of L&Z even when accounting for the 2- to 5-month time shifts (see green lines in Fig. 2, B and C).

Fig. 2

Trends in earthquake numbers are largely uncorrelated from trends in cumulative seismic moment, suggesting a different behavior of seismicity rates of small-magnitude versus large-magnitude events.

(A) Cumulative seismic moment (red) and cumulative earthquake number above M = 3 (black) with earthquakes M ≥ 5 highlighted by black arrows and vertical dashed lines. The inset shows a map of the study region and earthquake locations. W.O., western Oklahoma; C.O., central Oklahoma; (B) Injection (Inj.) (blue) and seismicity rates (red) in central Oklahoma as well as the triggering threshold (blue dashed line) for induced seismicity in the study of L&Z. None of the M ≥ 5 events occurred at times when injection rates exceeded the triggering threshold even when accounting for the proposed time delay (green lines) by L&Z. (C) Same as (B) for the western Oklahoma region.

A potential solution to this discrepancy would be to require longer time lags between injection rate changes and seismogenic response. We find that injection and seismicity rates are correlated with time lags of 9 to 14 months for both the entire study area and the central and northwestern regions similar to long time lags in previous studies (, ), but in contrast to the 2- to 5-month time lag of L&Z. This time lag approximately matches the observed onset of injection rate increase in 2012 and seismicity rate increase in early 2013 (Fig. 3 by L&Z). For the longer time lags, all the M ≥ 5 earthquakes other than Prague would occur within time windows when injection rates exceed the triggering threshold. However, this solution would also imply that the predicted seismic hazard in Fig. 1C in 2017 and 2025 remains elevated for an additional 4 to 12 months.

Although extending the time lag is a satisfactory solution based on the above statistics, there are physical and observational reasons to resolve the quandary by examining the model structure. Recent studies have shown that injection-induced earthquakes can involve a variety of physical mechanisms beyond the direct effective stress reduction in the seismogenic index model. These mechanisms include poroelastic stress transfer, aseismic creep, and static stress change (–). The different mechanisms have specific implications for the statistical approach to hazard assessment. For instance, the elastic stress changes external to a pressurized unit have been shown to induce seismicity at distances >30 km for some of the cases in Oklahoma () and would therefore modify the inferred stresses in the basement rock. In addition, some theoretical studies suggest that the magnitude distribution changes after injection and later earthquakes should be preferentially large (). The associated seismogenic response is likely to vary spatially, depending on local triggering mechanism and injection operations, so that aggregate, large-scale injection rate changes may not capture the underlying seismic hazard. Some of these mechanisms are acknowledged by L&Z as specifically missing from their model, and adding in more sophisticated physics may provide alternative solutions to the moderate magnitude issue discussed above.

In summary, we find that realistic parameter choices in the statistical model result in an increase in the probability of moderate earthquakes in 2017 by a factor of 2 over the reported probabilities by L&Z. In addition, the model is not able to predict several moderate earthquakes that have already occurred, and thus, the probability of future damaging earthquakes may also be underestimated.