## Abstract

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 (*1*) (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 (*2*). 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 (*2*, *3*). 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 *M*_{c} = 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*) = *R*_{0}/(*t*/*t*_{0})^{p} to predict future seismicity rates, where *R*_{0} is the seismicity rate at the start of significant injection rate decrease in 2015, *t*_{0} 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 (*3*). 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 (*4*).

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 *R*_{0} = 80 and *t*_{0} = 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 *M*_{c} within magnitude-dependent space-time windows (*5*). 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 (*R*_{0} = 22, *t*_{0} = 2). In either case, the resulting values of *p* are significantly below the assumed value of 2 by L&Z.

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 form(1)where *P*_{e} is the probability of observing an event with *M* ≥ *M*_{i}, λ 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 *M*_{c}. 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*. (*6*)]. 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).

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 (*7*, *8*), 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 (*9*–*11*). 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 (*7*) 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 (*11*). 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.

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:**We thank the Oklahoma Geological Survey and Oklahoma Cooperation Commission for making the seismicity and fluid disposal data publicly available.

**Funding:**This research was supported by the U.S. Department of Energy–Basic Earth Science program (award DE SC0015539).

**Author contributions:**All authors contributed to researching and writing this technical comment.

**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. Additional data related to this paper may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).