Aftershock deficiency of induced earthquake sequences during rapid mitigation efforts in Oklahoma

Abstract Induced seismicity provides a rare opportunity to study earthquake triggering and underlying stress perturbations. Triggering can be a direct result of induced stress changes or indirect due to elastic stress transfer from preceding events leading to aftershocks. Both of these processes are observable in areas with larger magnitude induced events, such as Oklahoma. We study aftershock sequences of M2.5 to M5.8 earthquakes and examine the impact of targeted injection rate reductions. In comparing aftershock productivity between California and Oklahoma, we find similar exponential scaling statistics between mainshock magnitude and average number of aftershocks. For events with M≥4.5 Oklahoma exhibits several mainshocks with total number of aftershocks significantly below the average scaling behavior. The sequences with deficient aftershock numbers also experienced rapid, strong mitigation and reduced injection rates, whereas two events with M4.8 and M5.0 with weak mitigation exhibit normal aftershock productivity. The timing of when aftershock activity is reduced correlates with drops in injection rates with a lag time of several days. Large mainshocks with significantly reduced aftershocks may explain decreasing seismicity rates while seismic moment release was still increasing in Oklahoma in 2016. We investigate the expected poroelastic stress perturbations due to injection rate changes within a layered axisymmetric model and find that stresses are lowered by 10s to 100s kPa within the injection-affected zone. For earthquakes induced by poroelastic stress-increase at several kilometers from wells, the rapid shut-in of wells may lead to elastic stress reductions sufficiently high to arrest unfolding aftershock sequences within days after mitigation starts.

Induced seismicity provides a rare opportunity to study earthquake triggering and underlying stress perturbations. Triggering can be a direct result of induced stress changes or indirect due to elastic stress transfer from preceding events leading to aftershocks. Both of these processes are observable in areas with larger magnitude induced events, such as Oklahoma. We study aftershock sequences of M2.5 to M5.8 earthquakes and examine the impact of targeted injection rate reductions. In comparing aftershock productivity between California and Oklahoma, we find similar exponential scaling statistics between mainshock magnitude and average number of aftershocks. For events with M≥4.5 Oklahoma exhibits several mainshocks with total number of aftershocks significantly below the average scaling behavior. The sequences with deficient aftershock numbers also experienced rapid, strong mitigation and reduced injection rates, whereas two events with M4.8 and M5.0 with weak mitigation exhibit normal aftershock productivity. The timing of when aftershock activity is reduced correlates with drops in injection rates with a lag time of several days. Large mainshocks with significantly reduced aftershocks may explain decreasing seismicity rates while seismic moment release was still increasing in Oklahoma in 2016. We investigate the expected poroelastic stress perturbations due to injection rate changes within a layered axisymmetric model and find that stresses are lowered by 10s to 100s kPa within the injection-affected zone. For earthquakes induced by poroelastic stress-increase at several kilometers from wells, the rapid shut-in of wells may lead to elastic stress reductions sufficiently high to arrest unfolding aftershock sequences within days after mitigation starts.

Introduction
Waste fluid disposal has been recognized as a major cause of induced earthquakes with several events above magnitude 5 in Colorado and Oklahoma (e.g. Ellsworth, 2013;Keranen et al., 2013;Rubinstein et al., 2014;Weingarten et al., 2015;Bao and Eaton, 2016;Yeck et al., 2017). Observations of larger magnitude induced events also extend to geothermal reservoir stimulation which likely led to a Mw5.4 event in South Korea in 2017 (Kim et al., 2018). Commonly, injection induced events are thought to be a result of effective stress reduction due to pressure diffusion along prestressed faults (e.g. Shapiro et al., 1997;Raleigh et al., 1976;Hsieh and Bredehoeft, 1981). However recent observations highlight the importance of aseismic processes and poroelastic coupling within the injection zone, which can intensify the seismogenic response * Corresponding author. at large distances and depths (e.g. Segall and Lu, 2015;Chang and Segall, 2016a;Barbour et al., 2017;Goebel et al., 2017a;Goebel and Brodsky, 2018;Guglielmi et al., 2015;Bourouis and Bernard, 2007). Differences in mechanistic views on induced events also impact our understanding of seismicity migration patterns, spatial decay and the choice of hazard mitigation strategies.
While the overall fraction of injection wells associated with induced earthquakes is small (about 1 to 10%), some regions are especially susceptible to induced stress changes (Weingarten et al., 2015;Göbel, 2015;Skoumal et al., 2018). The seismic hazard in regions with high induced event rates seems to be controlled by geological factors, particularly the relative distance between injection and the crystalline basement (e.g. Hincks et al., 2018;Skoumal et al., 2018). The placement of injection wells directly above the crystalline basement may be especially problematic because of the large spatial footprint of induced poroelastic stresses or fluid pressure diffusion in the presence of high-permeability faults Goebel et al., 2017a;Goebel and Brodsky, 2018  Both direct effects due to increased fluid pressures and poroelastic stresses as well as indirect effects due to event-event interactions and aftershock triggering can promote productive induced earthquake sequences. Fluid injection has been observed to raise aftershock productivity and background rates (Llenos and Michael, 2013;Zaliapin and Ben-Zion, 2016;Brodsky and Lajoie, 2013) such that both power-law aftershock decay and swarm-like sequences of events with similar magnitudes are observed (e.g. Sumy et al., 2014;. The complex interaction of the involved processes complicates induced seismicity mitigation. Previous mitigation efforts focused on restricting the minimum distance between wells and active faults (e.g. Rubinstein et al., 2018;Brown et al., 2017), as well as injection rate reductions and temporary well shut-in within specific distances from larger magnitude earthquakes in Oklahoma (Baker, 2016). In addition, injection directly into basement formations was suspected to be especially problematic for induced earthquake activity. In early 2015, the Oklahoma Corporation Commission issued a directive that injection wells would not be allowed to be completed into basement so operators were required to prove a well's true depth or be shut-in. 1 By July 2015, more than 120 wells were mandated to inject into the shallower Arbuckle formation instead of basement 2 ; however the effectiveness of shallower injection remains questionable considering the burst of large induced events in 2016 (Yeck et al., 2017). The diffuse nature of seismicity and disposal wells in Oklahoma motivated another mitigation strategy in the form of state-wide restrictions to monthly injection rates in 2016, which were set to 40% of 2014 rates (Baker, 2016).
To date, little is known about how ongoing induced seismicity sequences react to rapid injection rate reductions at broad regional distances (e.g. 10s of km), especially in areas with complex eventevent interactions and aftershock triggering. Here, we analyze aftershock productivity and decay rates after moderate earthquakes throughout the state of Oklahoma between 2009 and 2018. We specifically focus on seismicity rate changes that coincide with targeted injection rate reductions, providing a natural laboratory for the study of earthquake triggering.
We start by separating the seismicity into fore-main-aftershock sequences and background events. We then investigate average scaling statistics of aftershocks in Oklahoma and California and identify deviations from the average trend in Oklahoma. The detected anomalous aftershock sequences are then correlated with injection activity in the corresponding region. Lastly, we compute expected poroelastic stress changes using the recorded injection data and average well-earthquake distances.

Data and method
We analyze and compare seismicity records in Oklahoma and California based on data from the Oklahoma Geological Survey 1 https://www.occeweb .com /News /2015 /01 -30 -15EQ %20ADVISORY.pdf. 2 http://www.occeweb .com /News /DIRECTIVE -2 .pdf. and Southern California Earthquake Data Center. The primary focus of this study are seismicity rate changes in Oklahoma between 2009 and 2018, which are investigated through statistical analysis of standard and relocated earthquake catalogs (see Table 1 for details). The magnitude of completeness, M c , for each catalog is determined by fitting magnitude distributions with the Gutenberg-Richter relationship and determining the magnitude cut-off that minimizes the misfit between the exponential fit and observed distribution (Clauset et al., 2009). In addition to the earthquake data, we utilize well locations and injection rates archived by the Oklahoma Corporation Commission (OCC, 2018).

Aftershock detection and productivity relation
We separate the record of earthquake locations, origin times and magnitudes into clustered and background events. This separation is based on nearest-neighbor event-pairs which are determined from spatial-temporal distances of event pairs scaled by parent event magnitude (Zaliapin et al., 2008;Zaliapin and Ben-Zion, 2013). The observed nearest-neighbor distance distributions are compared with randomized Poissonian catalogs that have the same number of events and magnitude distributions as the original catalogs. Events are categorized as clustered at distances and times below the 99th percentile of the randomized catalogs (Fig. S1). The described analysis has the advantage that clustering characteristics are determined based on the observed data distributions without ad-hoc choices of space-time windows. Thus, there are few adjustable free parameters (here Gutenberg-Richter b-value, completeness magnitude and fractal dimension of seismicity). Previous studies highlighted that clustering characteristics are largely insensitive to parameter variations within plausible ranges (Zaliapin and Ben-Zion, 2013), and the method has successfully been applied to tectonic and induced events (Schoenball et al., 2015;Zaliapin and Ben-Zion, 2016).
The nearest-neighbor distances of clustered events can be used to assemble chains of triggered event families, which can be further classified into fore-main-and aftershocks based on relative timing and magnitudes of events within each family (Zaliapin and Ben-Zion, 2013). We find that the identified aftershocks within event families show Omori-type power-law decay of seismicity rates from the mainshocks (Fig. S3). This strong temporal clustering demonstrates that triggering mechanisms beyond purely injectioninduced forcing have to be considered to explain seismicity in Oklahoma (Llenos and Michael, 2013). In the following, we use the classified main-and aftershocks to determine productivity relations, i.e. the total number of aftershocks for each mainshock with a specific magnitude.
The average number of aftershocks, N AS , within specific mainshock magnitude bins (here we use 0.1 magnitude intervals) can commonly be described by an exponential distribution of the form (e.g. Reasenberg and Jones, 1989;Helmstetter et al., 2005): where M is the mainshock magnitude and α is the scaling exponent describing the relative number of aftershocks for small and large mainshock magnitudes. The exponent α is usually close to 1 but may vary as a function of tectonic regime and in the presence of geothermal activity (e.g. Helmstetter et al., 2005;Trugman et al., 2016).
In this study, background events are defined as events that occur at distances and times statistically indistinguishable from randomized earthquake catalogs. Since the randomization of the earthquake catalog is performed for the whole record, we implicitly assume a constant average background rate. Background rates in Oklahoma are thought to be driven by both tectonic and induced forcing (e.g. Llenos and Michael, 2013), giving rise to strong seismicity rate variations associated with injection rate increase and decrease. To address temporal background rate variations, we follow the strategy to first analyze the entire earthquake record and then compare the results with the analyses of records divided in 2 and 3 yr intervals for which the assumption of constant background rates is more realistic (see supplement Section S1).

Aftershock statistics
We compare productivity statistics for mainshocks with magnitude between 2.5 and 6 in California and Oklahoma. The average number of aftershocks within mainshock magnitude-bins can approximately be described by Eq. (1) within a magnitude range of 3.0 to 4.9 in California and 3.0 to 4.4 in Oklahoma. Corresponding values for α are 1.0 and 1.1 (Fig. 1). In addition to the slightly larger scaling exponent in Oklahoma, we observed overall higher aftershock numbers for a given mainshock magnitude. Several mainshocks above M4.4 are significantly deficient in aftershocks so that aftershock numbers are 2 to 10 times lower than expected from the average trend, and fall outside of the 95% confidence interval of the exponential fit. Mainshocks with deficient aftershock numbers include the M5.8 Pawnee, the M5.7 Prague, the M5.0 Cushing, and the M4.7 Cherokee event. Two clear exceptions from these observations are the M4.8 and 5.1 Fairview events with ∼1.5 times higher than average aftershock numbers. We tested the robustness of these observations by varying the magnitude of completeness between M c = 1.8 and 3.0 and using a relocated catalog from Schoenball and Ellsworth (2017). Results are consistent across the explored values of M c except for a change in α to 1.2 for M c = 3.0 in Oklahoma.
A potential issue when comparing aftershock productivity in Oklahoma, e.g. between the M5.7 Prague event in 2011 and the M5.1 Fairview event in 2016, are strongly varying background rates associated with substantial fluctuations in fluid injection rates. Comparing mainshocks in periods with low and high background rates may introduce a bias towards lower aftershock productivity when background rates are comparably low such as in 2011. We test the robustness of the determined productivity variations in sub-catalogs divided in 2 and 3 yr time intervals. The results confirm that the identified mainshocks with deficient aftershock activity are a robust trend in the data and are unlikely to be solely driven by background rate variations (Fig. S4).
One important difference between California and Oklahoma is the additional contribution of large-volume waste-water disposal. In Oklahoma, the increase in fluid injection activity has been associated with higher seismicity rates and event clustering (e.g. Ellsworth, 2013;Llenos and Michael, 2013) but mitigation and reduced injection rates may also be consequential for earthquake activity. We examine changes in fluid injection within 15 km from each mainshock and find that mainshocks with depleted aftershock sequences were also subject to rapid, targeted mitigation in the form of reduced injection rates or complete well shut-in (Fig. S6). Mitigation efforts can easily be identified within the injection rate time series by rapid, step-like rate reductions e.g. within days after the Cushing (Fig. 2) and Pawnee events (Fig. S6A). These observations are roughly in agreement with documented mitigation efforts by the Oklahoma Corporation Commission (Baker, 2016). An exception is the 2011 M5.7 Prague event which shows no indication of monthly injection rate reductions within 15 km (Fig. S6F). The lack of discernible mitigation effects may in part be due to the monthly time sampling of injection in 2011 (daily rates became available in September 2014). At the time, the Oklahoma Corporation Commission undertook localized mitigation by requiring several wells, near the Wilzetta Fault, to be shut-in for several days during which formation pressures were tested (Keller and Holland, 2013).
While the wider area around the Prague event hosted more than 10 high-volume disposal wells within 10 km from the main-shock, very little injection occurred near Fairview within 15 km of the M4.8 and M5.1 events in 2016. These are also the two mainshocks with productivity above the average trend. The induced events in this area were previously identified to be caused by farfield poroelastic stresses from field-wide injection activity at more than 30 km distance (Goebel et al., 2017a). Prior to the events, the regulators' mandate to reduce injection rates mainly focused on halting or reducing injection within 10 to 15 km from larger magnitude events. The regulator formulated a response plan for the Fairview area that included an injection rate reduction of 25% and injection depth reduction for two wells in November 2015 as well as a 18% injection rate reduction for 27 wells in the larger Fairview area in February 2016 (Baker, 2016). This response plan did not have a visible effect on cumulative, monthly injection rates of wells within 30 km of the M5.1 Fairview event until 100 days after the mainshock (Fig. 2). Even after 100 days, the decrease in injection rates was significantly less abrupt and of relatively smaller volume compared to other regions.
We further investigate potential effects of injection rate reductions on seismicity rates by analyzing aftershock rates after the M4.7 Cherokee event in November 2015. The aftershock rates are substantially lower than rates of the comparable M4.8 Fairview event but can approximately be described by a modified Omori relationship until 7 days after the mainshock (Fig. 3). The arrest of the aftershock sequence and deviation from Omori-type decay occurred about 5 days after the start of injection rate reduction and within 1 day of complete shut-in of two high-rate injectors at 4 km distance from the earthquakes (see Fig. S6C for inj. rates). The close association of well and earthquake locations and rate changes suggest a potential causal relationship.

Expected pore pressure and elastic stress changes due to mitigation
In the following, we numerically explore expected poroelastic stress changes as a result of long-term high-rate injection and rapid injection rate reductions in the area of notable earthquakes in Oklahoma. The effect of active mitigation is evaluated at 4 km distance and 7 days after the start of mitigation based on observations for the M4.7 Cherokee event (Fig. 3).
We determine poroelastic stress changes in a fully-coupled, layered, axisymmetric model comprised of the 500 m thick Arbuckle  injection formation, a 500 m thick impermeable cap-rock above the Arbuckle and the crystalline basement beneath (Fig. 4). The thickness of the crystalline basement is varied between 100 m to 10 km to explore the effect of boundary conditions. The model is built using the Comsol Multiphysics finite-element package (COM-SOL, 2017), and geometry and boundary conditions are shown in Fig. 4. Displacements and fluid pressures are continuous throughout the model domain.
The hydraulic and elastic parameters are taken from published values for the Arbuckle injection reservoir (Franseen and Byrnes, 2012;Morgan and Murray, 2015;Kroll et al., 2017). The cap rock in our model has strongly reduced permeability inline with a lack of observed upward fluid mobility. The crystalline basement has a permeability value three orders of magnitude lower than the Arbuckle formation (Table 2). We use a simplified, linear injection rate function based on observed rates near the Cushing events to examine the effect of a sharp decrease within days after the mainshocks. The model parameters for each layer are summarized in Table 2.
The main focus of our model are the poroelastic stress changes after mitigation efforts. Nevertheless, we also examine poroelastic stress increase before mainshocks, for which the model predicts an elastic stress increase between 10 −3 to 10 −1 MPa within the crystalline basement (Fig. 5). Elastic stress changes are highest within the upper 500 m to 1 km of the basement, as a result of larger stresses within the injection zone and continuity of displacements across the reservoir-basement interface. The amplitudes of resolved stress changes are comparable to previous modeling results for the Pawnee earthquakes (Barbour et al., 2017).
We evaluate the numerical modeling results in a benchmark test, using analytical solutions for pressure and elastic stress changes due to injection in a vertically confined injection layer from Helm (1994) and Rudnicki (1986) (Fig. S7). The numerical model accurately captures both the rapid near-well pressure decay Fig. 5. Poroelastic stress change in the radial direction (σ rr ) during high-rate injection before the Cushing event. Elastic stresses may have been sufficiently high to cause induced earthquakes at several kilometers distance and depth from the injection wells. out to about 2 km and the far-field response which is dominated by power-law elastic stress decay ahead of the diffusing pressure front (Fig. S8). A shut-in of the injection operation results in a reduction of pressures and elastic stresses within less than a kilometer from the well in both numerical and analytical solutions. At larger distances, on the other hand, fluid pressures continue to increase as a result of delayed diffusion effects, whereas elastic stresses show a notable decrease at distances beyond the largest pore-pressure gradients.
We further test this observation in the layered, numerical poroelastic model, focusing on poroelastic stress changes within the upper basement. We find a similar pattern of decreasing elastic stresses both close to the well and out to several kilometers distance beyond the pressure diffusion front (Fig. 6A). Interestingly, while the pressures within the injection formation behave analogously to the simple benchmark solutions, pressures within the basement continue to increase across all distances. This increase is a result of the low basement permeability which delays the diffusive pressure response beyond the evaluated time window of 7 days.
We find a maximum zone of influence of poroelastic stress decay of about 5 km within 7 days after injection shut-in (Fig. 6B).
At these distances, peak stress-decrease is expected to be ∼0.1 MPa within the injection zone and one order of magnitude smaller within the underlying basement. The relatively far reach of elastic stress decrease is partly driven by the rapid change in poroelastic loading within the overlying reservoir. Poroelastic effects may be further amplified within the basement if the contrast in elastic moduli across the sediment-basement interface exceeds the expected reduction in poroelastic coupling coefficient (Fig. S11).
We performed several sensitivity tests of the observed changes in fluid pressures and elastic stress. We test the influence of boundary conditions by varying the distance to the closest boundary (lower basement boundary) between 100 m to 10 km and find consistent near and far-field reductions in pressure and stress. We test different durations of injection intervals before shut-in between 7 to 300 days and basement permeability between 10 −17 and 10 −15 . The maximum spatial extent of elastic stress decrease depends partially on the duration of injection and permeability    6. A) Pore fluid pressure (blue) and radial elastic stress (red) as a function of horizontal distance from the well within the basement. Pressures and stresses are shown before (dashed curves) and 7 days after (solid curves) well shut-in. Poroelastic stress changes are dilational close to the well due to aquifer movements away from the well and change to compressive stresses at larger distances (see labels and black arrows). Note that at larger distances within the basement, pressures are still increasing while elastic stresses already started to drop. B) Spatial extent of absolute changes in pressure and stress after well shut-in. Pressure reduction is shown within the injection reservoir since basement-pressures are still increasing within 7 days after shut-in.
structure. Longer injection durations and high permeability favor a larger zone of influence of compressive stress reduction.
A key result from the numerical and analytical models is the effect of mitigation on fluid pressures and elastic stresses at kilometer scales within days from well shut-in. For the tested range of parameters, we find that fluid pressures may be reduced within the injection reservoir up to 1.5 km within days from shut-in. The response in the underlying basement is dominated by poroelastic effects. The corresponding elastic stress may decrease at distances of 5 km or more.

Discussion
In this work, we presented observations of targeted injection rate reduction and coinciding earthquake sequences with low aftershock productivity in Oklahoma. The systematic temporal and spatial correlations between mitigation and deficient aftershock sequences suggests a likely physical connection. Mitigation efforts started as early as one day after the mainshocks, and affected wells within a distance range of 4 to 15 km. Our observations suggest that rapid mitigation can lower earthquake activity over large areas and potentially reduce induced seismic hazard from fluid-injection operations. Previous work suggests a potential connection between injection rates or volumes and the expected number of induced earthquakes or seismic moment release (van der Elst et al., 2016;McGarr, 2014;Weingarten et al., 2015). This effect may also explain some local differences in seismic activity between different areas in Oklahoma.
Our analysis of earthquake rate changes after rapid mitigation provide some insight into the distance range at which injection rate changes can influence seismicity. We found several examples, such as the Pawnee and Cushing events, for which active mitigation within 15 km from the mainshocks was sufficient to reduce the rate of aftershocks. The examples from the Fairview area showed that continued high-rate injection at 30 to 40 km distance likely kept the number and rate of earthquakes high throughout the aftershock sequence. More detailed observational constraints on maximum expected distances at which mitigation is effective is difficult in Oklahoma because many wells at varying distances likely contributed to the observed earthquakes (e.g. Animation S1).
Our observations suggest that not all mitigation efforts in Oklahoma were successful or included wells at sufficiently large distances. Examples of documented mitigation efforts include a 25% injection rate reduction within two wells in the Fairview area in November 2015, which failed to reduce seismicity rates or pre-vent the M5.1 Fairview event in February 2016 (Baker, 2016). The lack of effectiveness of mitigation in the Fairview area may be explained by the presence of far-field poroelastic stresses caused by densely-spaced high-rate injectors at more than 15 km from the earthquakes (Goebel et al., 2017a). Similarly, the Cushing area initially experienced only localized regulatory actions in the form of 25% injection rate reduction, which decreased seismicity rates at the end of 2015 and beginning of 2016. However, seismic activity again surged following injection rate increase in the broader Cushing area which eventually led to the M5.0 Cushing event in November 2016 (Fig. S12, Animation S1). Injection activity in the area was subsequently reduced at a larger scale leading to the described reduction in aftershock rates. These observations highlight the overall complexity of the coupling between injection and seismicity rates, which complicates an a priori prediction of the effectiveness of mitigation efforts.
Poroelastic effects can induce deep and distant earthquakes by fluid injection and may also impact earthquake rates during active mitigation efforts. The here described observations of reduced seismic activity after targeted mitigation may be explained by poroelastic coupling between fluid pressures and solid stresses, which allows rapid stress reductions within days after well shut-in. The difference in spatial-temporal scales of solid stresses vs. fluid pressures in a fully coupled model have been explored previously to explain post-shut-in seismicity rate jumps (Segall and Lu, 2015). In contrast to the previous work, observations from Oklahoma suggest optimally-oriented faults that enhance induced earthquake activity during elastic stress increase. Rapid mitigation, on the other hand, is associated with elastic stress reductions which potentially remove previously induced shear stresses, leading to the observed aftershock deficiency.
The present work further highlights specific conditions that promote poroelastic stress contributions to induced earthquake triggering. Such conditions include: (i) densely-spaced injection wells that collectively act as a finite source (Goebel et al., 2017a), (ii) effective elastic coupling at large distances and depths (Chang and Segall, 2016b;Barbour et al., 2017;Goebel and Brodsky, 2018), and (iii) elastic stress variations at short time scales after injection rate changes as shown here and in Segall and Lu (2015).

Potential mechanisms for rapid aftershock reduction
The effect of targeted injection rate reductions on aftershock rates can be considered a natural laboratory for the study of earthquake triggering. Lower aftershock productivity after mitigation Fig. 7. Schematic image of two scenarios that may explain the observed aftershock deficiency after targeted mitigation efforts in Oklahoma. Aftershock rates are shown in red and fluid injection rates in blue. Insets show maps of hypothetical earthquake and well distributions. Left: If injection rate reductions directly affect an unfolding aftershock sequence, deviations from Omori-type behavior are expected. Right: If injection rate reductions mainly lower the rate of distributed background events, aftershock decay is expected to remain Omori-like but at lower productivity level. The new productivity level with induced forcing removed is similar to purely tectonic aftershock decay. may be explained by two processes: (1) If fluid injection activity leads to a large-scale increase in background stressing rates, this may also lead to higher triggering susceptibility and more productive aftershock sequences (Dieterich, 1994). Rapid mitigation removes the additional induced forcing so that aftershock productivity returns to tectonic values (Fig. 7 right). High background rates complicate a clear statistical separation of background and aftershocks within the aftershock zone. A potential contribution of background rate variations is supported by the observed higher overall productivity in Oklahoma compared to California, assuming tectonic and induced stressing rates in Oklahoma vs. only tectonic rates in California. This is in agreement with previous work which found higher aftershock productivity when injection activity is increasing (Llenos and Michael, 2013). Moreover, after mitigation some of the aftershock sequences in Oklahoma exhibit productivity values comparable to purely tectonic sequences in California (Fig. 1b).
(2) Another potential explanation is connected to the direct effect of mitigation on fault stresses. Rapid mitigation may significantly lower the stress level on faults around the mainshock, thereby inhibiting or arresting a developing aftershock sequence (Fig. 7 left). Such arrests are observed during the Cushing and Cherokee aftershock sequences. The presented modeling results indeed suggest that induced stress variations affect large-enough areas to stop an unfolding aftershock sequence. Thus, induced stress reductions of 0.1 MPa in the reservoir and 0.01 MPa in the basement are potential upper thresholds for aftershock triggering stresses.

Increase in seismic moment during periods of seismicity rate decrease
The observed earthquake sequences with deficient numbers of aftershocks due to mitigation may explain a discrepancy between earthquake rate reduction and increased moment release in Oklahoma in 2016 (Fig. 8) (Yeck et al., 2017;Goebel et al., 2017b). While state-wide earthquake rates started to decrease in early 2016, seismic moment release kept increasing until the end of 2016. The missing aftershocks due to mitigation after large earthquakes in 2016 may explain the overall rapid decrease in statewide seismicity rates. It also suggests that state-wide injection rates provide limited information about the expected induced hazard which may be dominated by substantial, local variations in injection volumes and pressures. The inset shows the regions most affected by injection-induced earthquakes in Oklahoma. B: Evolution in cumulative number of earthquakes (black curve) and seismic moment (red curve) above M3. Note that while seismicity rates started decreasing in March 2016, seismic moment release kept increasing until the end of 2016 so that about half of the total moment within the 9 yr time span was released in 2016.

Conclusion
We presented detailed statistical and numerical analyses of the connection between changes in injection rates and aftershock statistics in Oklahoma. Aftershock sequences in Oklahoma are generally more productive than comparable sequences in California, potentially caused by elevated ambient stress level due to fluid injection operations. Several larger mainshocks with M≥4.5 in Oklahoma exhibit anomalously low productivity. The corresponding regions experienced mitigation and injection rate reductions within 4 to 15 km from the mainshocks. The systematic temporal and spatial correlations between mitigation and deficient aftershock sequences suggests a likely physical connection.
Numerical and analytical models show that mitigation effects are dominated by fluid pressure reductions close to injection wells within the Arbuckle and elastic stress decreases at larger distances within the basement. Both induced stress increase before the mainshocks and stress reduction during the aftershock sequences can be viewed as natural laboratory for the study of earthquake triggering. Assuming that mitigation efforts directly affected unfolding aftershock sequences, we estimate upper thresholds for aftershock triggering stresses of 0.1 MPa in the Arbuckle and 0.01 MPa in the basement where most earthquakes occur. These estimates are maximum values and actual stress changes may be even lower. Our observations suggest that rapid mitigation can affect earthquake rates over large areas and potentially reduce induced seismic hazard within days from fluid-injection operations.