Pore-pressure diffusion, enhanced by poroelastic stresses, controls induced seismicity in Oklahoma.

Induced seismicity linked to geothermal resource exploitation, hydraulic fracturing, and wastewater disposal is evolving into a global issue because of the increasing energy demand. Moderate to large induced earthquakes, causing widespread hazards, are often related to fluid injection into deep permeable formations that are hydraulically connected to the underlying crystalline basement. Using injection data combined with a physics-based linear poroelastic model and rate-and-state friction law, we compute the changes in crustal stress and seismicity rate in Oklahoma. This model can be used to assess earthquake potential on specific fault segments. The regional magnitude-time distribution of the observed magnitude (M) 3+ earthquakes during 2008-2017 is reproducible and is the same for the 2 optimal, conjugate fault orientations suggested for Oklahoma. At the regional scale, the timing of predicted seismicity rate, as opposed to its pattern and amplitude, is insensitive to hydrogeological and nucleation parameters in Oklahoma. Poroelastic stress changes alone have a small effect on the seismic hazard. However, their addition to pore-pressure changes can increase the seismicity rate by 6-fold and 2-fold for central and western Oklahoma, respectively. The injection-rate reduction in 2016 mitigates the exceedance probability of M5.0 by 22% in western Oklahoma, while that of central Oklahoma remains unchanged. A hypothetical injection shut-in in April 2017 causes the earthquake probability to approach its background level by ∼2025. We conclude that stress perturbation on prestressed faults due to pore-pressure diffusion, enhanced by poroelastic effects, is the primary driver of the induced earthquakes in Oklahoma.

The recent increase in the number of earthquakes in the central and eastern United States since 2008 is attributed to the massive deep subsurface injection of saltwater (e.g., refs. 1–4). The spatial proximity of the seismicity to the injection wells and that many of these events are preceded by month to years of high-volume fluid injection suggest a link between the observed seismicity and injection (e.g., refs. 5–10). The elevated rate of seismicity and chance of damaging earthquakes raise broad societal concerns among the public and regulators (1). Despite significant efforts in improving the monitoring capability for detecting induced seismicity (5) and understanding the underlying mechanism (4, 10–13), much more needs to be done to successfully quantify and forecast the associated time-varying seismic hazard (14–18).

Since 2008, central and northern Oklahoma has experienced a 900-fold increase in seismicity (Fig. 1 and ), including 4 major events with magnitude greater than or equal to 5: 2011-11-09, moment magnitude (Mw) 5.7, Prague; 2016-02-13, Mw 5.1, Fairview; 2016-09-03, Mw 5.8, Pawnee; and 2016-11-07, Mw 5.0, Cushing (). In response to the seismicity surge during 2015 and to mitigate hazards, the regulator reduced the total volume of disposed brine within areas of elevated seismicity in 2016 to less than 40% of the 2014 total amount (). Despite injection reduction and decline of seismicity rate, seismic moment release soared within the injection-regulation zones, culminating in several major events in late 2016, such as the Pawnee and Cushing earthquakes (19). This observation indicates that the relation between fluid injection and the associated induced seismic hazard is complex and that the underlying physics associated with the process of induced seismicity is not captured in the current induced earthquake-forecasting models (14–18).

Fig. 1.

Fluid injection and seismicity in Oklahoma from 1995 to 2017. (A) Map showing the locations of M3+ earthquakes after declustering and Arbuckle wastewater disposal wells. Blue circles are recent earthquakes during 2008–2017, and black circles indicate the historical earthquakes before 2008. Triangles represent the injection wells, colored and scaled according to average monthly injection rates. The 4 M5+ earthquakes are shown as black stars and are summarized in . Black lines are mapped faults (56). The black dashed line is the Nemaha Fault Zone and Uplift dividing CO (red box) and WO (blue box). A, Inset shows the location of the study region (black dashed box). (B) Time series of total monthly injection volume (red) and a histogram of all recorded seismicity (black) after declustering in CO. (C) Time series of total monthly injection volume (blue) and a histogram of all recorded seismicity (black) after declustering in WO. Ave., average; CO, central Oklahoma; EQ, earthquake; inject., injection; WO, western Oklahoma.

Overall, earthquake hazard is proportional to the seismicity rate (14), which is determined by changes in crustal stress (20). As the fluid is injected into the target formations and pore pressure diffuses, the stress field is perturbed, reducing the shear strength of prestressed faults and promoting their slip (4, 13, 21). Moreover, the maximum magnitude of injection-induced events is controlled by several factors, including total injected volume, regional tectonics, and local hydrogeology (22–25). The total seismic moment is correlated with both injection volume (26) and basement depth (27). Further, the occurrence of moderate- to large-magnitude induced earthquakes is determined by background tectonic stress and basement fault structures (25). This body of evidence highlights that a successful forecasting model requires full integration of the physics governing the processes of fluid diffusion in a poroelastic medium and induced earthquake nucleation (28).

Results and Discussion

The monthly injection volumes at 867 high-volume wells spanning January 1995 to April 2017 were obtained from the Oklahoma Corporation Commission (9). We removed wells with significant data gaps, leaving us 715 wells within north-central Oklahoma (Fig. 1 and ). The Oklahoma Geological Survey compiled the earthquake catalog used here. The seismic moment released in both central Oklahoma (CO) and western Oklahoma (WO) was dominated by the magnitude (M) 3+ earthquakes that mostly nucleated within the upper few kilometers of the crystalline basement, with far fewer in the Arbuckle group (2, 15). In CO, although the fluid injection commenced in 1995 and increased over time, the seismicity only began in 2008 and peaked in 2015 (Fig. 1). In WO, injection began in 2005, but the sharp increase in seismicity occurred in 2013, coinciding with a rapid rise of fluid-injection rate (Fig. 1). Focusing on the mainshocks, we followed the declustering approach of Reasenberg (29) recently applied to seismicity catalog in Oklahoma (16, 30, 31) to remove the dependent earthquakes and identify the events directly linked to deep fluid injection ().

Availability of subsurface stratigraphy, local hydrogeology, seismic tomography, seismicity migration analysis, and Earth tide strain analysis allowed us to characterize parameters of a linear, layered poroelastic Earth model ( and ). Poroelasticity describes the mechanical coupling between pore pressure and poroelastic stresses, and the fluid-diffusion rate is controlled by the hydraulic diffusivity. The permeable Arbuckle sandstone–carbonate formation in Oklahoma possibly has different hydraulic diffusivities within the geologically distinct CO and WO regions (32). Studies have indicated a diffusivity of 1–2 m2/s for CO (2, 9, 33, 34) and 2–7 m2/s for WO (35–37). We chose diffusivities of 1.5 and 4.0 m2/s for CO and WO, respectively. As discussed later, these values enable the best representation of observed seismicity in the study area (see Fig. 3). Moreover, records from wells drilled into the crystalline basement (38) provided estimates of spatially variable basement distance to the injection well bottom (). This distance was suggested to correlate with seismic moment release (27) and was accounted for to determine the injection depth for the layered poroelastic model ().

Fig. 3.

Observed and predicted M3+ earthquakes in CO and WO through the physics-based approach. (A–C) Observed (Left) and simulated (Center, only PP; Right, PP and PES) earthquake magnitude–time distribution (blue dot) and the associated earthquake count (black) and density (red) for M3+ earthquakes in CO. (D–F) Observed (Left) and simulated (Center, only PP; Right, PP and PES) earthquake magnitude–time distribution (blue dots) and the associated earthquake count (black) and density (red) for M3+ earthquakes in WO. CO, central Oklahoma; EQ, earthquake; PES, poroelastic stresses; PP, pore pressure; WO, western Oklahoma.

Using the Earth model and time series of injected fluid volume, we solved for the spatiotemporal evolution of the pore pressure and poroelastic stresses in the crust (4, 28). Accounting for both poroelastic stresses and changes in pore pressure, we calculated the total Coulomb Failure Stress (CFS) (Fig. 2 and ) with a northeast-oriented receiver fault geometry at the Arbuckle-basement boundary, where it is hydraulically connected to the underlying crystalline basement through high-permeability fault zones ( and ). We found that changes in pore pressure dominate the spatial and temporal patterns of CFS change. The simulated pore pressure, poroelastic stresses, and CFS change all increased before the hypothetical injection shut-in in April 2017 and decayed afterward. The corresponding temporal evolution of averaged Coulomb Stressing Rate (CSR) is dominated by pore-pressure diffusion and has a nonlinear pattern (Fig. 2). Two peaks characterize the time series of CSR in the CO region at the end of 2008 and 2014; in the WO area, there was an extended period of elevated CSR during 2013–2015 (Fig. 2). Both patterns are consistent with the temporally variable injected volumes (Fig. 1 ). Our sensitivity tests using 2 optimal, conjugate fault orientations with northeast and southeast strikes, suggested by Alt and Zoback (39) and Holland (22), showed that the CSR time series remains unchanged given either of the conjugate faults, suggesting that pore-pressure diffusion dominates the CFS.

Fig. 2.

Simulated time series of average Coulomb Stressing Rate (CSR) and seismicity rate in CO (A) and WO (B). Solid blue and red curves are the time series of average CSR, accounting for both pore pressure and poroelastic stresses along northeast fault geometry and the resulted relative seismicity rate, respectively. Three (two) snapshots of spatial distributions of CSR and logarithmic relative seismicity rate are also shown for CO (WO). Dashed blue and red lines indicate the individual effect of pore pressure and poroelastic stresses on the CSR and relative seismicity rate. Black solid circles are the incremental M3+ seismicity either from the beginning of the observation period in 1995 or the timing of the previous snapshot. Model parameters are provided in . Note that the curves for southeast trending fault, which are not shown, are the same as that for northeast trending fault. CO, central Oklahoma; WO, western Oklahoma; Ch., change.

Assuming that the background stressing rate and characteristic relaxation time are the same for both CO and WO, we applied a rate- and state-dependent seismicity rate model (13, 20, 28) to simulate the change in earthquake count relative to the background seismicity, as a result of imparted CFS changes. This model is, however, applicable only if fault systems are critically stressed before injection. The Arbuckle group, where most injection occurred, is naturally underpressured throughout most of the midcontinent (40, 41). Thus, during the early stage of injection, the fluid was used to compensate for the pressure deficit, and only when the pressure was high enough to propagate into the basement was it able to trigger seismicity. This hypothesis is consistent with the observations that elevated seismicity in CO and WO regions began ∼13 and ∼8 y after injection began, respectively (Fig. 1 ). Accounting for these delays and informed by the time series of CSR, we solved for the temporal evolution of relative seismicity rates in the CO and WO regions (Fig. 2 and ), as well as at the epicenters of the 4 M5+ earthquakes using the fault plane solutions obtained from the focal mechanisms () through an adapted version of the seismicity rate model. This model (Eq. in ) considers a regional critical time (2008 for CO and 2013 for WO) when the seismicity rate starts to deviate from its background value ( and ). In 2015, the time series of seismicity rate for both regions were characterized by a significant peak, while that of CO showed an additional, smaller peak in 2010 (Fig. 2). The snapshots of the spatial distribution of the modeled seismicity rate showed outward propagating seismicity front with a decreasing amplitude after injection shut-in at some high-volume wells (). This observation is consistent with the notion that after injection, shut-in fluid continues to diffuse away from injectors until pressure equilibrium is achieved. Before that, the state of effective stress in the medium changes continuously as the pore pressure evolves in space and time. Noteworthy, we found a causal relationship between the locations of observed seismicity and that predicted by our physics-based model, despite considering only homogeneous values for the model parameters in each region. Moreover, the model shows that the Prague Mw 5.7 and Fairview Mw 5.1 earthquakes happened during increasing seismicity rate, while the Pawnee Mw 5.8 and Cushing Mw 5.0 earthquakes occurred when the seismicity rate was at its maximum (); however, the application of such physics-based models for forecasting large events is yet to be tested.

Several hydrogeological and nucleation parameters control the time-dependent behavior of the simulated seismicity rate. We performed sensitivity tests using the CO datasets to investigate the effects of diffusivity D, characteristic relaxation time t (by varying ), background stressing rate , and critical time t () on the seismicity rate model. We altered only one parameter at a time, fixing others to a reference value (). At the local scale, varying all 4 parameters can change the timing, amplitude, and shape of the predicted seismicity rate curve (), as suggested by previous studies (13, 28). This is due to the complex evolution of fluid diffusion and earthquake-nucleation processes in time and space. At the regional scale (), however, the timing of the average seismicity rate curve is insensitive to the variation of these 4 parameters, because the seismic hazard is proportional to the aggregated seismicity rate over the whole region. In particular, in Oklahoma, the injectors are widely distributed, and the significance of individual diffusion time in the observed average seismicity rate is largely reduced.

Moreover, given that the characteristic relaxation time is far longer than the duration of our model simulation, the observed lag between the seismicity rate and monthly injection rate (Figs. 1 and 2) is likely due to the gradual increase of the average cumulative CFS (Fig. 2). This means that during the initial stage of injection-rate reduction (e.g., early 2015), the cumulative CFS continued rising despite the countering effects of continuous fluid diffusion. Also, the long characteristic relaxation time had a small impact on the seismicity rate increase that primarily corresponded with the increasing cumulative CFS.

The variation of the hydrogeological and nucleation parameters can influence the amplitude of the average seismicity rate, with t and t having the greatest and lowest sensitivity within reasonable ranges, respectively (). The bimodal pattern of average seismicity rate curve in terms of peak-to-peak ratio, which is defined by the ratio of the peak amplitudes of seismicity rate, is sensitive to D and t, but insensitive to t and . Conversely, the amplitude of the seismicity rate showed a larger sensitivity to D than to t. Noteworthy, the half-year uncertainty chosen for the t sensitivity test is reasonably well-constrained by observation (31). Table 1 summarizes the overall conclusions reached following the sensitivity tests on model parameters.

Table 1.

Summary of model parameter sensitivity tests

Parameter	Regional seismicity rate curve
	Timing	Maximum amplitude	Peak-to-peak ratio
D	Low	Intermediate	High
Aσ¯	Low	High	Low
τ˙_0	Low	Intermediate	Low
tcrit	Low	Low	High

The regional seismicity rate curve is characterized by 3 properties of timing, maximum amplitude, and peak-to-peak ratio. The degree to which they depend on each parameter is qualified with 3 levels of low, intermediate, and high sensitivity based on the performed tests. D, hydraulic diffusivity; , background stressing rate; t, critical time. Characteristic relaxation time t is varied by altering .

These sensitivity tests indicated that the parameters used for the physical modeling of the fluid diffusion and earthquake nucleation processes were reasonable and well-constrained by observations. Moreover, the trade-off between diffusivity and characteristic relaxation time is negligible at the regional scale. An additional test without considering a critical time resulted in an entirely different seismicity rate curve (), implying that the critical time is required for an accurate assessment of the seismicity rate. Considering that the background stressing rate and characteristic relaxation time in WO are same as that of CO, for WO, we only performed additional sensitivity tests investigating the effects of hydraulic diffusivity and critical time () and obtained similar results to that of CO.

The change in poroelastic stresses contributes to the process of induced seismicity in addition to pore pressure (13). However, its role in evaluating induced seismic hazard ranges from being ignored (17) to being considered as the primary driver of large-footprint seismicity (10) when injecting into permeable layers. Our model analysis allowed us to evaluate the contribution of poroelastic stresses in a complex setting such as Oklahoma. We separately calculated the CSR due to pore pressure and poroelastic stresses and investigated the associated seismicity rate changes. We found that pore-pressure diffusion was the dominant factor determining the seismicity rate in contrast to the effect of poroelastic stresses (Fig. 2). But, adding the contribution of poroelastic stresses to the pore pressure amplified the estimated seismicity rate by 6-fold and 2-fold for CO and WO, respectively, compared with the case that only considered the contribution of pore-pressure diffusion. This is because seismicity rate is sensitive to stress fluctuations and is exponentially proportional to CFS change (13, 33), meaning that a small shift in CFS may result in a substantial increase in earthquake number due to the nonlinear seismic response to the stressing rate. The more prominent poroelastic effect on seismicity rate in CO is due to the slower fluid-diffusion rate.

We next evaluated the Gutenberg–Richter frequency–magnitude relationship (42) for the recorded seismicity before the seismicity increase in the CO and WO regions (). We found that the k values (annual rate of 10) are ∼102.71 and ∼101.7 y−1 for the 2 regions, respectively, with the same b value of 1.09 (a and b are the parameters determining Gutenberg–Richter law). Availability of Gutenberg–Richter parameters that characterize the background seismicity rate before seismicity increase and the relative seismicity rate change obtained through modeling allowed us to estimate the probability density function of absolute seismicity rate as a function of earthquake magnitude and time (). Estimating the probability density function of given earthquake magnitudes (e.g., M3+) at any time interval, the magnitude–time space can be discretized and randomly sampled through iterations (Fig. 3). In CO, 2 peaks of increased seismicity rate and earthquake magnitude were accurately recovered in 2010 and 2015, considering the joint effect of pore pressure and poroelastic stresses. In WO, both predicted and observed earthquake magnitude–time distributions were characterized by a single peak in 2015. In both zones, the predicted numbers of M3+ earthquakes were similar to that of observed seismicity. Moreover, we tested the effect of poroelastic stresses on earthquake magnitude–time distribution and found that ignoring poroelastic stresses could lead to substantial decreases of forecasted earthquake number and magnitude and thus underestimation of seismic hazard (Fig. 3).

We further estimated the annual earthquake magnitude exceedance probabilities in CO and WO (Fig. 4 and ). The annual probability for exceeding M5 increased with time until 2015, from <1% in 2008 to 43% in 2015 for CO and from <1% in 2012 to 45% in 2015 for WO. Due to the mandated injection-volume reduction in 2016, our model predicts a significant decrease in the probability of exceeding M5 (down to 23%) in WO. However, the reduction of probability is negligible in 2016 for CO. The reason for different responses to injection-volume reduction is that the earthquake exceedance probability is highly time-dependent because of the temporally variable nature of fluid injection, fluid diffusion, and earthquake nucleation. The nucleation parameters of characteristic relaxation time and background stressing rate were the same for CO and WO in our model; however, the injection volume reduction in WO was more significant than that of CO. Also, the Earth model in WO was characterized by a larger hydraulic diffusivity, which resulted in the more rapid diffusion of fluid and decrease in seismicity rate and thus exceedance probability.

Fig. 4.

Annual earthquake magnitude exceedance probabilities in central Oklahoma (CO) (A) and western Oklahoma (WO) (B). Injection data are updated to April 2017 when a hypothetical injection shut-in is imposed for simulation. In both regions, seismic hazard is expected to be reduced to 2008–2012 levels by 2025 after the hypothetical injection shut-in based on the poroelasticity and rate-and-state friction models.

Our model predicts a probability of ∼10% to exceed M5.8 during 2016 in CO where the M5.8 Pawnee earthquake occurred in September 2016. Injection-induced seismicity follows a Gutenberg–Richter relationship, which generally holds for the distributed events. One limitation of applying this law to estimate the maximum earthquake that can be hosted by the local faults is that the linear relationship associated with the Gutenberg–Richter law may not capture the background magnitude statistics of the rare large-size earthquakes. The distribution of b-value-associated with fluid-induced seismicity can vary significantly in space (43, 44). Moreover, some induced sequences exhibit deviations from the Gutenberg–Richter distribution (45–47). For these reasons, the calculated probability should be viewed as an upper bound estimate to exceed M5.8 at the regional scale, which can be complemented by local seismic hazard evaluation efforts accounting for locations of preexisting faults (). Other factors, such as the length scale of the seismogenic faults and fault connectivity, should also be considered for physics-based forecasting.

Assuming a hypothetical injection shut-in in April 2017, the probabilities within both CO and WO continue to decrease, as expected from postinjection decaying patterns of both pore pressure and seismicity rate, and approach the background tectonic level by approximately 2025. This indicates that the probability of large earthquakes may not decrease immediately following injection shut-in, primarily due to the time-dependent nature of fluid diffusion, which continues after injection stops (4, 13, 33). Factors affecting the probability decay include injection rate, background stress, and local hydrogeology.

Our model incorporates coupled porous media flow modeling and a physics-based seismicity rate model to statistically forecast the magnitude–time distribution of induced earthquakes as well as to calculate the earthquake magnitude exceedance probability, which is physically different from other approaches that predict the induced seismicity frequency (15–18). Seismogenic index theory, which considers injection rates to reconstruct a Gutenberg–Richter law, has been applied to M3+ induced earthquakes in Oklahoma (15). However, this theory is suggested to only apply for magnitude ≤2.0 earthquakes and nondecreasing injection rate (48). The follow-on modification of seismogenic index theory (17) replaces the injection rate with pore pressure obtained from an uncoupled groundwater flow model. This modification has no direct physical interpretation, and the estimated seismogenic index for the same area is significantly different from one study to another (15, 17). As indicated by Langenbruch et al. (17), the distribution of seismogenic index is controlled by the pore pressure obtained from a flow model, and thus it resembles a free parameter that one can tune to fit observed seismicity, raising the issue of model uniqueness and sensitivity. Norbeck and Rubinstein (16) use a reservoir-compression model (e.g., no flow condition) to calculate pore pressure, which then is combined with seismicity rate theory. This model overestimates the pore pressure by up to 2-fold, and the pore pressure cannot decrease because fluid-diffusion processes are neglected. Thus, the seismicity rate is likely overestimated. Ignoring fluid diffusion also prevents the model from forecasting future seismicity rate. A recent model (18) combines fluid-diffusion modeling and seismicity rate theory in WO, but the distributed injectors are superimposed and simplified to a broad zone of a single injection. Since the diffusion length and time are related to hydraulic diffusivity, and the length scale is greatly increased, this simplified model requires a diffusivity of 44–277 m2/s, which is 1–2 orders of magnitude larger than that in the Arbuckle layer. Moreover, these cited studies ignore the effects of poroelastic stresses on the forecasted induced seismicity at the regional scale, whereas previous studies (13, 33) and our results (Fig. 2) show that poroelastic stresses are essential in the process of induced seismicity. Our approach relies on choices in rock mechanical properties and coupling a fluid-diffusion simulation and seismicity rate model that each exhibit strong nonlinearities as the associated parameters vary (). The involved parameters can be estimated and constrained independently from geomechanical, seismic, and geodetic data (e.g., ref. 49), rather than through a calibration procedure, which is required for empirical approaches (e.g., refs. 15 and 17).

Assessing the time-varying seismic hazard due to fluid injection is critically important. Successful efforts to forecast fault activation require accurate quantification of the physics governing the evolution of crustal stresses and seismicity rate (1, 28, 50). Despite the improvements in seismic monitoring capacity and the resulting decrease in the magnitude detection threshold (7, 51), estimates of induced earthquake probability still require models that account for stress and pressure changes. Our results highlight the critical role of fluid diffusion in a poroelastic medium to understand the temporal evolution of the induced seismic hazard. Also, continuously updated information about the probability of a future earthquake is essential for successful operational earthquake forecasting. Thus, the ability to link the evolution of pore fluid pressure and stress change to seismicity rate change presents a proactive approach to quantifying the seismic hazard associated with fluid injection (13) and developing frameworks for operational induced earthquake forecasting (52).

Methods

Poroelastic Modeling.

We employed a coupled poroelastic model to calculate the spatial and temporal evolution of poroelastic stresses and pore pressure due to fluid injection. The theory of poroelasticity accounts for the coupling between the deformation of the porous medium and evolution of the pore fluid pressure. This means a change of pore pressure can deform rocks and vice versa. The full governing equations of linear poroelasticity combine the Navier–Cauchy equation, which is derived by substituting poroelastic constitutive equations into the equilibrium equation, and the diffusion equation, which is obtained by combining Darcy’s law and the continuity equation (53). The governing equations relating the deformation field and pore pressure p, both of which are a function of position and time t, are given (53, 54)

where ∇ is the gradient operator and ∇⋅ is the divergence operator, G is the shear modulus, υ is the drained Poisson ratio, α is the Biot effective stress coefficient (the change in fluid volume per unit volumetric change in medium under drained condition), Q is the Biot modulus, χ is the mobility coefficient defined by the ratio of intrinsic permeability and dynamic fluid viscosity, is the body force per unit bulk volume acting on solid medium, and q is the volumetric fluid injection rate per unit bulk volume. To characterize a linear poroelastic medium, 5 independent parameters are needed, including G, υ, undrained Poisson ratio υ, hydraulic diffusivity D (the ratio of hydraulic conductivity and specific storage), and Skempton coefficient B (the change in pore pressure per unit change in confining pressure under undrained conditions). Parameters α, χ, and Q can be uniquely determined by using these 5 parameters (54). The system of governing Eqs. and is solved by imposing a mechanical boundary condition of zero traction and a flow boundary condition of zero excess pore pressure at the half-space surface (55). Additional details about model setup and implementation are given in .

Seismicity Rate Modeling.

Both laboratory experiments and the rate-and-state friction law predict that a small change in shear or normal stresses may cause a significant change in fault slip rate. Dieterich (20) developed a framework describing the evolution of seismicity rate as a function of background seismicity rate and CFS change. A simplified version was given by Segall and Lu (13) relating the relative seismicity rate R(, t) (rate of seismicity relative to the background seismicity rate) to the CSR . However, the inherent assumption associated with the Dieterich (20) model is that the background stress is sufficiently high relative to the shear resistance that the background stressing rate leads to a nonzero seismicity rate. Given that the Arbuckle formation is naturally underpressured (40, 41), an amount of fluid is initially needed to compensate for the pressure deficit before the seismicity rate increases. Thus, to solve for the seismicity rate associated with the imparted CSR, we considered a regional critical time t (when seismicity rate starts to deviate from its background value):

where is the background stressing rate, A is a constitutive parameter in the rate-and-state friction law, is the background effective normal stress, and is the characteristic relaxation time. Note that this seismicity rate theory relates the seismicity rate history to Coulomb stressing history, meaning the seismicity rate is determined by the entire history of stressing rate. Also, R is a value 1. When the stressing rate is equal to the background stressing rate, R = 1. The regional critical time t is fixed for both CO and WO constrained from observed seismicity for both zones (Fig. 1 ). The piecewise ordinary differential equations (Eq. ) are solved numerically. Additional details about CFS calculation and model setup are given in .