# Seismological Research Letters

- © Seismological Society of America

## ABSTRACT

The 2016 *M*_{w} 5.8 Pawnee earthquake occurred in a region with active wastewater injection into a basal formation group. Prior to the earthquake, fluid injection rates at most wells were relatively steady, but newly collected data show significant increases in injection rate in the years leading up to earthquake. For the same time period, the total volumes of injected wastewater were roughly equivalent between variable‐rate and constant‐rate wells. To understand the possible influence of these changes in injection, we simulate the variable‐rate injection history and its constant‐rate equivalent in a layered poroelastic half‐space to explore the interplay between pore‐pressure effects and poroelastic effects on the fault leading up to the mainshock. In both cases, poroelastic stresses contribute a significant proportion of Coulomb failure stresses on the fault compared to pore‐pressure increases alone, but the resulting changes in seismicity rate, calculated using a rate‐and‐state frictional model, are many times larger when poroelastic effects are included, owing to enhanced stressing rates. In particular, the variable‐rate simulation predicts more than an order of magnitude increase in seismicity rate above background rates compared to the constant‐rate simulation with equivalent volume. The observed cumulative density of earthquakes prior to the mainshock within 10 km of the injection source exhibits remarkable agreement with seismicity predicted by the variable‐rate injection case.

Animations of the evolution of pore pressure, vertical displacement, and cylindrical tensor strains in the injection simulation domain, and a set of input commands for running simulation A using poel in ASCII format.

## INTRODUCTION

Injection‐induced earthquakes are well documented in Oklahoma (Ellsworth, 2013; Keranen *et al.*, 2013; Walsh and Zoback, 2015), where the rate of earthquakes has increased sharply since 2009. Studies of the recent sequences in Oklahoma explain the increased seismicity as a result of increased pore pressures reactivating pre‐existing faults (e.g., Keranen *et al.*, 2014). Hubbert and Rubey (1959) demonstrate this behavior analytically, whereby increases in fluid pressure can cause slip on pre‐existing structures. The influence of fluid pressure on the frictional stability of faults in the subsurface has been argued to explain the earliest known cases of injection‐induced seismicity at Rocky Mountain Arsenal and Rangely (Evans, 1966; Healy *et al.*, 1968; Raleigh *et al.*, 1976; Hsieh and Bredehoeft, 1981).

Although direct pore‐pressure effects are typically believed to be the cause of injection‐induced earthquakes, poroelastic effects may be important to consider. For instance, recent studies consider the possibility that aseismic processes resulting from injection can be responsible for changes in seismicity on faults with frictional instability (Bourouis and Bernard, 2007; McClure and Horne, 2011; Segall and Lu, 2015; Deng *et al.*, 2016; Fan *et al.*, 2016; Norbeck and Horne, 2016a). In general, this has been difficult to establish because of a lack of observations of pore‐fluid pressure in proximity to earthquake sequences, but measured slip histories from controlled fluid‐injection experiments (Guglielmi *et al.*, 2015) demonstrate this effect clearly. Furthermore, numerical simulations (Zhang *et al.*, 2013; Chang and Segall, 2016; Norbeck and Horne, 2016b) demonstrate that fault‐related variations in permeability can influence the location and timing of induced seismicity; careful analyses of seismicity and injection patterns (McNamara *et al.*, 2015; Yeck, Weingarten, *et al.*, 2016) support this conclusion.

In this article, we explore the possibility that poroelastic effects arising from fluid injection influenced the *M*_{w} 5.8 Pawnee earthquake (Fig. 1). Notably, new injection data in nearby Osage County show significant variability in the years preceding the earthquake, whereas other wells show relatively constant rates of injection (Fig. 2). This provides an opportunity to compare the influence of the direct pore‐fluid pressure effect and the poroelastic effect on the generation of seismicity, focusing on the difference between constant‐rate and variable‐rate wastewater disposal. We compare simulations of the time‐varying poroelastic stresses and pore pressure resulting from transient fluid injection with the equivalent constant‐rate injection history based on an equivalent volume of fluid. Our results show clear differences in spatial and temporal patterns in Coulomb failure stress, depending on whether or not poroelastic effects are considered. Our model predicts more than an order of magnitude increase in the maximum seismicity rate prior to the mainshock, about a factor of 3 larger than the constant‐rate result because of differences in the stressing rate history.

## DATA AND METHODS

### Injection Data

There are numerous wells within 15 km of the mainshock epicenter, but we restrict our study to the nine wells (Table 1) with injection rates into the Arbuckle Group greater than 0.3 kbbl/mo (Fig. 1), wells that could have had the largest influence on the seismicity. Much of this injection activity has been either fairly small compared with similar wells in Oklahoma or continued at a relatively constant rate (Fig. 2). This includes the two closest wells to the mainshock epicenter, OLDHAM and SCROGGINS, ∼2.4 and ∼5.4 km from the mainshock, respectively; these wells were injecting at relatively constant rates between 146 and 211 kbbl/mo since 2012. Because we are interested in looking at the interplay between the poroelastic effects and hydrologic effects of variable fluid‐injection rates, we examine the two closest wells with the largest changes in injection behavior: SOUTH BEND 2A‐11 (OS6301) and RICE 2A‐9 (OS6379), located ∼6.5 and ∼8.7 km from the epicenter of the mainshock, respectively (SOUTH BEND is ∼3.2 km from RICE). More specifically, injection rates at SOUTH BEND and RICE increased from 0 to 288 and 404 kbbl/mo, respectively, over a period of ∼6 months beginning in early 2013. This rapid change is followed by a steady reduction in rates to level the characteristic of the average in this period for all wells within 15 km of the mainshock, ∼92 kbbl/mo, and following the *M*_{w} 5.8 mainshock these wells were shut‐in. (All wells within a zone of regulatory action were shut‐in after the mainshock.) These data are used in our poroelastic simulations.

### Poroelastic Reservoir Modeling

Our model domain comprised four major lithologic units representing the relatively uniform lithology in north‐central Oklahoma (Luza and Lawson, 1980; Christenson *et al.*, 2011). Injection is assumed to be entirely within the Arbuckle Group, a laterally extensive basal carbonate formation group which we treat as a single layer between 1300 and 1900 m depth (below ground surface [bgs]) based on information from well‐completion diagrams. The injection interval overlies semi‐infinite crystalline basement rock and is overlain by a confining layer of finite thickness. An unconfined layer representing sedimentary units overlies the confining layer. Figure 3 shows the general model setup. We assume that all layers behave as homogeneous linear poroelastic media with varying responses, which we describe in detail next.

We solve the general boundary value problem of fluid‐mass flux into a system governed by the equations of linear poroelasticity (compare with, Wang, 2000); these are based on isothermal equilibrium, Hooke’s law of elasticity, and Darcy’s flow law (conservation of mass). We use the spectral element method of Wang and Kümpel (2003), which calculates time‐varying cylindrically symmetric solutions in a layered half‐space. We assume an unconfined free‐surface boundary condition (*p*=0) and a volumetric source time function for injection based on real data (described above). The injection simulation captures the injection‐rate variability in Osage County from 2013 through late 2016, predominately from the wells RICE and SOUTH BEND (see Fig. 2).

The particular poroelastic parameters specified for each layer in the half‐space include Skempton’s coefficient *B*, the drained and undrained Poisson’s ratios *ν* and *ν*_{u}, the hydraulic diffusivity *D*, and the elastic shear modulus *μ*. Poroelastic solutions are affected by four coupled parameters (1a)(1b)(1c)(1d)(Kümpel, 1991; Wang and Kümpel, 2003; Segall, 2010), in which *α* is Biot’s coefficient of effective stress, *β* is the bulk compressibility, *χ* is the Darcy conductivity, and *λ* is Lamé’s first constant of elasticity. We treat these coupled parameters as homogeneous within each layer except in the basement, for which we use a variable hydraulic diffusivity that decreases logarithmically from the top of the layer to 8 km depth (then is constant at all depths below). We combine this depth‐dependent diffusivity with a homogeneous shear modulus to represent a typical depth‐dependent permeability profile in fractured crystalline rock (Ingebritsen and Manning, 2010). This type of permeability structure is known to control hydraulic conductivity in critically stressed crystalline rock (Barton *et al.*, 1995; Townend and Zoback, 2000).

We perform two simulations of the same variable injection rates with the goal of isolating the sensitivity to elastic parameters in the four layers (see Table A1). The first, referenced as simulation A, uses elasticity constants representative of the general rock lithology, adapted from Roeloffs (1996), Wang (2000), and Jaeger *et al.* (2007). The second, referenced as simulation B, fixes the parameters in each layer to those for the basement layer in simulation A. In both simulations, the Skempton’s coefficient is set to a single value across all layers, *B*=0.75; this may be unrealistically high for some units in the sedimentary section (Lockner and Stanchits, 2002), but is generally valid for granitic basement (Wang, 2000; Jaeger *et al.*, 2007), and direct pore‐pressure observations (Kroll *et al.*, 2017) suggest a comparable value for the Arbuckle Group.

We set the drained Poisson’s ratio to *ν*=0.25, and the undrained Poisson’s ratio to *ν*_{u}=0.38. With the Skempton’s coefficient stated earlier, this gives an effective stress coefficient of *α*=0.75. This value of *α* is higher than the value generally assumed for bulk granite, *α*=0.47, because the undrained Poisson’s ratio is higher than expected for the same rock type (see Wang, 2000, their table C.1); however, published values of *α* are generally calculated from laboratory measurements (e.g., Rice and Cleary, 1976) and are thus subject to the influence of measurement variability. We use statistical inference to quantify the conditional probability distributions of each parameter because uncertainties are not generally reported, to our knowledge. Based on 10^{4} realizations in a Monte Carlo simulation of *α* (equation 1a), with the values used in simulation A for the basement layer (Table A1) and two assumed levels of relative uncertainty (5% and 10%), we find broad distributions of *α*: standard deviations are 0.11 and 0.23 for 5% and 10% relative uncertainty, respectively, and lower bounds are 0.53 and 0.27 at 95% confidence. The details of these distributions demonstrate that, although the inferred value of *α* is perhaps unlikely for bulk material, it is not implausible statistically. Furthermore, the effects of pervasive hydraulically conductive (critically stressed) fractures in crystalline rock include elevated compliance and specific storage (Barbour, 2015; Xue *et al.*, 2016); thus, although our simulations may represent a limiting case of bulk material properties, they likely capture the enhanced sensitivity to external stresses (e.g., injection of wastewater) expected for the Precambrian basement.

### Seismicity Rate Modeling

The coupled geomechanical and hydrologic simulations provide spatial and temporal distributions of stress and fluid pressure throughout the model domain. We assume that the orientation and location of the fault plane that hosted the mainshock follows the distribution of relocated aftershocks (Yeck, Hayes, *et al.*, 2016) and is vertical. The top of the fault plane begins at the intersection of the injection layer and basement. Then, we calculate the Coulomb stress changes along the fault (Pollard and Fletcher, 2005) caused by fluid injection (2)in which Δ*τ* is the change in shear stress, Δ*σ* is the change in normal stress on the fault, Δ*p* is the change in fluid pressure within the fault zone, and *f* is the static friction coefficient. Tensile normal stresses are taken as positive in this sign convention; hence, positive Δ*σ*+Δ*p* indicates unclamping of the fault. Furthermore, positive Δ*τ* represents left‐lateral shear; hence, positive Δ*τ*_{f} represents failure‐promoting shear stress changes, assuming the maximum principal stress is compressional and oriented at N85°E (Walsh and Zoback, 2016). At each point along the fault, the changes in shear and normal stresses are obtained by rotating the principal stress tensor from the poroelastic simulation into the direction of the fault pole. The fault surface is assumed to have a uniform static coefficient of friction *f*=0.65, as laboratory studies of basement rock in Oklahoma indicate (Carpenter *et al.*, 2016).

Following the empirical approach introduced originally by Dieterich (1994) and expanded upon by Segall and Lu (2015) and Chang and Segall (2016), we use a rate‐and‐state friction framework to model the evolution of seismicity rate along the fault in response to changes in the pre‐existing state of stress. In this model, seismicity rate evolves according to (3)in which is the Coulomb stressing rate, is the background stressing rate, is a characteristic timescale over which the seismicity rate returns to background levels, *A* is the rate‐and‐state direct‐effect parameter, is the effective normal stress, and *R* is defined as the ratio between the current seismicity rate and a reference seismicity rate. When *R*=1, the predicted seismicity rate is equal to the background rate. We solve equation (3) numerically using an explicit Runge–Kutta method (Griffiths and Smith, 2006). The rate‐and‐state frictional properties used in our model are representative of granite at hydrothermal conditions (Blanpied *et al.*, 1995) and are listed in Table A2.

An advantage of this frictional model is that it is possible to consider arbitrary loading conditions. In this study, we examine the effects of a relatively rapid change of injection rates over a period of roughly 3.7 yrs in two wells near the fault, in Osage County; seismicity rate calculations are based on the simulated Coulomb stress changes (from the poroelasticity model) associated with the idealized disposal well operations. We repeat these calculations with the direct pore‐pressure effect isolated, wherein only pore‐pressure changes are taken from the solution; we verified this approach with a traditional hydrologic model using an equivalent layered structure (P. Hsieh, personal comm., 2016). It is important to recognize that these idealizations are not a representation of the full history of injection leading up to the mainshock; but they can help explain the effects of variations in injection rate. It is also important to recognize that the seismicity rate model used in this study does not simulate the nucleation, rupture, or arrest processes of individual earthquake events in a deterministic manner: changes in seismicity rate can only be interpreted in terms of their effect on earthquake statistics.

## RESULTS

### Stress Changes from Injection

The poroelastic injection simulations show features consistent with general hydrologic models. That is, the strongest changes are seen in pore‐fluid pressure changes, which are strongly dependent on the hydraulic diffusivities of each rock layers. The presence of a confining layer promotes fluid diffusion radially in the Arbuckle‐Simpson layer and, to a lesser extent, in the vertical direction into the basement, resulting in a pore‐pressure distribution that decays both in depth and radial distance from the injection source. Based on supplementary simulations, the realistic depth‐decreasing diffusivity in the semi‐infinite basement layer enhances the depth‐decay in pore‐pressure diffusion rates.

Despite the general consistency between our poroelastic simulations and hydrologic models, the poroelastic model includes significant variations in the stress field at all locations that cannot be replicated in a decoupled or partially coupled hydrologic model. Figure 4 shows snapshots of the evolution of stresses and pore pressure for simulation A (Table A1). There is a nonlinear relationship between the change in excess pore pressure Δ*p* and the change in mean stress Δ*σ*_{m}, defined by changes in principal stress (Δ*σ*_{1}≥Δ*σ*_{2}≥Δ*σ*_{3}) (4)indicating the influence of laminar fluid‐volume flux (*χ*∇*p*) through the model domain. Figure 4 also shows elevated differential stresses (5)that are initially concentrated near the injection source during the injection transient and eventually focus around the mainshock centroid depth. Although elastic stresses initially dominate pore‐pressure changes, increases in Δ*σ*_{d} are on the order of the size of Δ*σ*_{m} and Δ*p*, indicating that poroelastic shear stresses exhibit considerable influence on the stress state of the rock, even at distances much greater than the injection source depth. The strongest changes in pore pressure and stresses are close to the injection source, and the sharpest contrasts in pore‐pressure and stress changes are related to contrasts in material properties. In the sedimentary layer for example, significant stress changes are generated but very little pore‐pressure diffusion occurs because of the relatively low diffusivity of the confining layer.

At depths of the mainshock centroid, the effect of material contrasts controls regions of high Δ*σ*_{d} that are localized. In particular, the depth‐dependent hydraulic diffusivity structure has the strongest influence on localization, as the following comparison between the two simulations demonstrates. In simulation B, the elastic parameters are set equivalently in all layers, based on the values for the basement fault (Table A1). Even though this fixed‐elasticity equivalent of simulation A is unrealistic, especially in shallow sedimentary layers, it is an instructive exercise for understanding the sensitivity of the simulated stress changes to variations in elastic parameters of the rock. Figure 5 compares the evolution of stresses and pore pressure at the centroid of the mainshock for both simulations; the simulations are nearly identical at each point in time, with the exception of the minimum principal stress at early times and the pore pressure at later times. We find that elastic stresses dominate pore‐pressure effects at the mainshock centroid when the principal stress aspect ratio (6)(e.g., Bott, 1959) is less than 1/2. The transition from *r*≤1/2 to *r*>1/2 occurs ∼1.7 yrs after the variable‐rate injection begins, around when peak differential stresses occur, with peak pore‐pressure changes occurring roughly 0.5 yrs later. Elevated levels of differential stress at centroid depths are apparent in both the homogeneous elasticity and layered elasticity cases, but the maximum difference between the two solutions is less than 1 kPa; hence, even with homogeneous elastic structure, stress localization is controlled by the hydraulic diffusivities in each layer. The pronounced dependence of both pore‐pressure and stress changes on hydraulic diffusivity is a direct consequence of interrelated hydromechanical parameters (Kümpel, 1991; Barbour and Wyatt, 2014).

In addition to the effects of elasticity, there is another difference between the poroelastic and hydrologic models; the timing of peak pore‐pressure changes at centroid depths is different. Although a hydrologic model predicts a pressure increase everywhere in the model domain assuming boundary continuity, the poroelastic simulation predicts a decrease in pore pressure in the basement caused by instantaneous elastic strain in the rock at the initiation of injection that is roughly 0.5 yrs long. This is a phenomenon sometimes observed in shallow water wells in crystalline aquifers (e.g., Wolff, 1970). Despite the small size of the pressure reduction relative to pressure increases at later times, the observed lag represents a temporary fault‐stabilizing pressure reduction (Δ*p*<0) prior to imminent pressure rise (Δ*p*>0) from pressure diffusion. The duration of this delay depends on distance to the source and exhibits sensitivities to the difference between the drained and undrained Poisson’s ratios (Rice and Cleary, 1976) *ν* and *ν*_{u} and the initial conditions of the simulation. The magnitude of the initial change in injection rate is particularly important, but the opposite effect occurs in the transition to zero injection; thus, in certain faulting regimes it is theoretically possible to mitigate damaging effects of rapid shut‐in by carefully tapering injection rates (Segall and Lu, 2015).

### Seismicity Rate Evolution

In the seismicity rate change model, the extent of the Sooner Lake fault is represented as a 4 km (along depth) by 8 km (along strike) vertical plane with a strike of N65°E and an orientation of ∼45° to the injection source. The fault tip is located at ∼4.2 km relative to the injection source (RICE) and 1.9 km below the surface. The Coulomb stress change acting on the fault is used to inform the seismicity rate model, and the stressing rate is assumed to be constant between the weeklong timesteps from the injection model, and is not sensitive to small variations in dip angle. In Figure 6, we show the distribution of Coulomb stresses resolved on the fault, and the associated seismicity rate changes, at several snapshots in time during simulation A. There is a general correspondence between the spatial pattern of stress and seismicity rate, with some aftershocks following isostress contours. The maximum Coulomb stress change is on the order of Δ*τ*_{f}≈40 kPa at the mainshock centroid, and the maximum seismicity rate change is *R*≈16, suggesting more than an order of magnitude increase above the background rate.

In our simulations, the maximum predicted seismicity rate occurs roughly one year following the period of peak injection rate. During 2015, the predicted seismicity rate declines steadily until early 2016, at which point the rate stabilizes at an average rate of *R*≈5. The predicted seismicity rate tracks the stress changes closely and resembles a time‐shifted and smoothed version of the injection rate, which appears to be a consequence of the sensitivity of equation (3) to the stressing rate for these loading conditions. We do not observe an immediate correlation between the predicted peak seismicity rates and the timing of the seismic sequence. Comparing the seismicity rate calculation using the full poroelastic solution with the pore‐pressure effect isolated—only fluid pressure changes affect the state of stress along the fault—reduces the maximum rate to *R*≈3. Furthermore, the value of *R* appears to be relatively insensitive to elastic structure, which again implies that hydraulic diffusivities exhibit the strongest influence on the fully coupled solution.

## DISCUSSION

The frictional stability of the Sooner Lake fault prior to the 2016 *M*_{w} 5.8 sequence may have been reduced by active wastewater injection in close proximity to the fault. Although pressure diffusion is indeed the dominant mechanism for reducing the effective stress on the fault, the magnitudes of shear and normal stresses induced by coupling between elastic deformation of the solid matrix and pressure diffusion are comparable to pore‐pressure changes, to within a factor of 2. The distribution of Coulomb failure stresses predicted by our simulations between the centroid depth and the top of the fault at the same radial distance (shaded region in Fig. 7) implies that even though high‐permeability pathways may be important to consider, they are not necessary to explain the relative timing between this injection transient and the seismic sequence; a stronger control on this timing is due to gross permeability structure, which depends on both the hydraulic and elastic properties.

We have not yet considered the effects of steady injection rates at wells closer to the rupture (e.g., wells in Pawnee County). Appealing to the same physical mechanism used to inform the seismicity rate changes, we assume that injection at these wells must have influenced stressing rates on the fault; however, without strong changes in injection rate their influence may have been subtle, albeit important. For example, the seismicity rate model (equation 3) is highly sensitive to stressing rates as we found in the variable‐rate case; but, after nearly a decade of steady injection, as annual records of injection volume extending back to 2006 indicate (e.g., OLDHAM), stress and pore‐pressure accumulation would be relatively slow. To test this, we simulated the effect of long‐term constant injection beginning in 2011, two years prior to the previous simulations. To make the results comparable, we use a constant‐rate injection with the same total volume of fluid as in the variable‐rate case. We find that constant‐rate injection loads the fault slower than in the transient cases (Fig. 7). Additionally, this constant‐rate simulation produces a much different temporal change in *R*, which roughly approaches steady‐state conditions with the square root of time, and decreases rapidly after shut‐in (Fig. 8). This comparison strongly suggests that long‐term injection may have been responsible for a gradual loading of the fault to the point where it primed the fault for failure triggered by the short‐term high‐rate injection transient. In the absence of transient injection activities, however, frictional failure may still have occurred at a much later time.

In addition to computing the rate increase *R*, we examine the integral of the rate increase to test our results against observed seismicity (Fig. 8). This quantity, referred to as Σ*R*, represents the forecast increase in earthquake rate at an individual point in time. In some sense, it is analogous to an earthquake probability density function, whereby the probability of an earthquake is the highest at times for which *R* is highest. Accordingly, we could consider Σ*R* to be an empirical cumulative density function. Examining this function for simulation A (Fig. 8) shows that although the instantaneous probability of an earthquake may be lower at later times in the simulation, the overall probability for a given earthquake or a collection of earthquakes to have occurred is still increasing.

Foreshock activity and temporal variations in magnitude distributions (Walter *et al.*, 2017) lend support to this prediction. We also note that seismicity was not observed within a 10 km radius of the modeled injection source until late‐2013, when a clear increase in rate occurred. Even though the majority of these events are not considered foreshocks, their cumulative density function resembles the predicted *R* densities for the variable injection‐rate case. In general, detectable seismicity occurs approximately one year after variable‐rate injection begins, and is nearly coincident with peak injection rates. If the poroelastic stresses modeled for the Sooner Lake fault are equivalent for similarly oriented fault planes around the injection source, assuming cylindrical symmetry is justified (Segall and Lu, 2015), then these events may be related to changes in injection rate. And if earthquake magnitudes are random in time (e.g., Kagan and Knopoff, 1987), the overall probability of a large earthquake owing to increases in injection is expected to increase even though it may be instantaneously lower at later times.

The apparent connection between injection rate, predicted seismicity rate, and observed cumulative density is consistent with the results of Weingarten *et al.* (2015) that demonstrate that, although large total disposal volumes are not a necessary condition for elevated seismicity, high injection rates do imply a higher probability of inducing seismicity. One possible explanation for the lack of seismicity in other parts of Osage County, where moderately high rates of fluid injection occur (Murray, 2014), may be from regional variations pore‐pressure changes that occur because of compartmentalization of the sedimentary section, which includes the Arbuckle Group. This is supported by myriad faults with appreciable throw that show up in comparable hydrostratigraphy in south‐central Oklahoma (Christenson *et al.*, 2011; Mashburn *et al.*, 2013) and a general deepening of basal formation depths to the southwest (Luza and Lawson, 1980).

Fluid pressures in a compartmentalized reservoir will remain high in response to injection mass if flow is inhibited. In a region with highly variable injection histories (Weingarten *et al.*, 2015) and a set of mapped faults with both diverse orientations and preferential regional directions (Holland, 2013; Alt and Zoback, 2016; Walsh and Zoback, 2016), the effect of compartmentalization and its influence on heterogeneous pressure distributions would be especially pronounced. Furthermore, injection in Osage County is predominantly on the west side of a frictionally stable north–south‐trending fault that intersects the Sooner Lake fault (Fig. 1c), which suggests that pore‐pressure changes may have localized until transient postseismic flow (Sibson, 1994) caused pressure equalization (Manga *et al.*, 2016). In any case, the effects of reservoir compartmentalization in the Arbuckle can be investigated directly with data from a network of downhole pore‐pressure monitoring systems.

One complicating factor in these interpretations is the possibility of anisotropic or temporal variations in permeability. Active crustal shear strain maintains shear dilatancy over hydrothermal sealing processes, diminishing the response to quasi‐static deformation in favor of enhanced flow in high strain‐rate plate‐boundary regions (Barbour, 2015). Even though strain rates in Oklahoma inferred from intraplate seismicity (Anderson, 1986) and geodetic measurements (Calais *et al.*, 2006) are low compared to those at plate boundaries, for example, they may still outpace hydrothermal sealing processes in the region (Nathenson and Guffanti, 1988). Spatial variations in permeability exist, especially near faults, but these may be compensated by variations in specific storage, which implies that hydraulic diffusivity may be relatively uniform over large regions (Xue *et al.*, 2016).

A shortcoming of the seismicity rate model used in this study is that it calculates changes relative to the background seismicity rate. The question then becomes: What is the background rate? In the Pawnee study area, it is difficult to infer a reliable estimate of the background seismicity rate because of the extremely small sample size of documented earthquakes prior to 2009 and because all major seismic sequences in Oklahoma have been on unmapped faults (Yeck, Hayes, *et al.*, 2016). The model results we show must be interpreted within the context of a probabilistic analysis; however, it is not possible to quantify the effect on the seismic hazard without establishing a background seismicity rate, which is subject to current debate. The background seismicity rate in Oklahoma has been estimated (Ellsworth, 2013; Langenbruch and Zoback, 2016) to be roughly one *M*_{w}≥3 earthquake per year inside the areas of interest defined by the Oklahoma Corporation Commission (∼26×10^{3} km^{2}). If the seismicity rate for our model domain can be inferred from this regional estimate by accounting for the relative area covered by the model (∼314 km^{2}), the local background seismicity rate is roughly one *M*_{w}≥3 earthquake per century. Our modeling results suggest an increase in seismicity from a background rate of roughly 0.01 earthquakes per year to roughly 0.15 earthquakes per year in mid‐2014. As we mentioned above, seismicity rates did increase in the years prior to the Pawnee earthquake, around the time of peak‐predicted seismicity rate (Fig. 8), but these observations are limited and possibly tied to changes in detection thresholds (i.e., network changes).

The present set of poroelastic injection‐ and seismicity‐rate models do not incorporate any detailed information about the initial state of stress acting on the fault, the fault’s initial proximity to shear failure, stress transfer effects during earthquake rupture, or variations in stress owing to either geometric variation, frictional asperities, or aseismic slip. In future studies, it would be worthwhile to use deterministic earthquake mechanical and hydrological models to identify the conditions that would have influenced earthquake nucleation and rupture during the Pawnee sequence (e.g., McClure and Horne, 2011; Norbeck and Horne; 2016b).

## CONCLUSIONS

Simulations of time‐varying Coulomb failure stresses on the fault induced by injection activities indicate that the combined effect of stress changes associated with a high‐rate fluid‐injection transient and long‐term injection may have influenced the timing and location of the Pawnee earthquake. Even though diffusion is the principal mechanism for transmitting pore‐fluid pressure changes in basement rock, the effects of strain coupling between the rock and fluid are significant and should not be ignored. Seismicity rate calculations using a rate‐and‐state friction model are sensitive to the rate of change of Coulomb failure stresses and predict increases of more than an order of magnitude above background rates. Predicted seismicity rate increases generally resemble a lagged smoothed version of the injection‐rate history, and the observed seismicity leading up to the mainshock follows the cumulative density function of rates predicted by the variable‐rate injection simulation. If seismicity rates are affected by poroelastic stress changes, our models suggest that they may remain elevated for at least a year following shut‐in of injection wells after the earthquake, which emphasizes the importance of incorporating accurate geomechanical properties of the subsurface in understanding the hazards associated with injection‐induced seismicity, especially in low strain‐rate environments.

## DATA AND RESOURCES

Injection data for Osage County are from the Environmental Protection Agency (EPA) and are available in Ⓔ the electronic supplement to this article; data for other counties are from the Oklahoma Corporation Commission (OCC) Oil and Gas Database (http://www.occeweb.com/og/ogdatafiles2.htm). Mapped faults are from the Oklahoma Geological Survey (http://www.ou.edu/content/ogs/data/fault.html). We used poel (i.e., Wang and Kümpel, 2003), version 2012, for numerical simulations of injection. Regional earthquakes are from the Advanced National Seismic System (ANSS) earthquake catalog (http://www.ncedc.org/anss/catalog-search.html). All websites were last accessed March 2017.

## ACKNOWLEDGMENTS

We thank Nancy Dorsey at the Environmental Protection Agency (EPA) for help accessing injection records from Osage County and thank Xiaowei Chen and Paul Hsieh for helpful discussions. Elizabeth Cochran and Paul Hsieh graciously provided thorough preliminary reviews on short notice. We thank two anonymous reviewers for comments that led to improvements in the clarity and presentation of this work. Jack Norbeck was supported by the Stanford Center for Induced and Triggered Seismicity and a Mendenhall Postdoctoral Fellowship.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

## Appendix. SIMULATION PARAMETERS

In Table A1, we give the layered poroelastic parameters used in simulations of injection. Simulation A represents the case of variable elasticity parameters within each layer based on representative bulk properties (Wang, 2000; Jaeger *et al.*, 2007), whereas simulation B represents the limiting case in which elasticity parameters are fixed for the entire domain based on parameters for granitic basement. In both simulations, the hydraulic diffusivities are variable, but equivalent; thus, the effective permeabilities are different between the two simulations, except in the basement rock.

In Table A2, we give the parameters used in the seismicity rate calculations, including the state parameter, background stressing rate, and effective normal stress.