We utilize quantitative risk assessment (QRA) to calculate the conditional probability of slip on mapped faults in response to injection-related increases in pore pressure in north-central Oklahoma (USA) where widespread injection of produced saltwater has triggered thousands of small to medium-sized earthquakes in the past 7 yr. The conditional probability incorporates the uncertainty in each Mohr-Coulomb parameter (stress tensor, pore pressure, coefficient of friction, and fault orientation) through QRA. The result is a cumulative distribution function of the pore pressure required to cause slip on each fault segment. The results can be used to assess the probability of induced slip on a known fault from a given injection-related pore pressure increase. After dividing north-central Oklahoma into six study areas, we invert earthquake focal plane mechanisms in each area to constrain the orientation and relative magnitude of the principal stresses. The QRA identifies the potential for slip on the fault that produced the M 5.6 Prague earthquake in 2011 and the northeastern extension of a mapped fault associated with the M 5.1 Fairview earthquake sequence that occurred in early 2016, and, had the 289°-striking fault of the September 2016 M 5.8 Pawnee event been mapped, it would have been identified as potentially active.
One mechanism by which injection-related pore pressure increases can trigger seismic slip is well understood in the context of Mohr-Coulomb failure criteria (e.g., National Research Council, 2013). Wastewater disposal has been associated with induced earthquakes at a number of sites in the central and eastern United States (Horton, 2012; Keranen et al., 2013). In north-central Oklahoma, the increase in seismicity that began in 2009 (Fig. 1) is caused by large rates of produced-saltwater disposal in the area where the earthquakes are occurring (Walsh and Zoback, 2015). We present a methodology to quantify the conditional probability of inducing slip on known faults in response to injection-related pore pressure changes given reasonable assumptions using a Mohr-Coulomb slip criteria. While other physical processes can affect triggering, we identify the faults of most concern.
We first review existing information on faults in Oklahoma. We then discuss the available earthquake focal mechanisms within each of six study areas to constrain the stress field in each area. The study areas are defined to encompass most of the recent seismicity and the best seismic coverage (McNamara et al., 2015). The focal mechanism inversions yield information about both stress orientation and relative magnitudes. We next demonstrate how quantitative risk assessment (QRA) can be used to assess the conditional probability that a known fault might slip in response to injected-related pore pressure increases through a Mohr-Coulomb mechanism. A key component of this analysis is constraining the uncertainties associated with each of the geomechanical parameters. It should be noted that within the seismic hazard community, the term “risk” has a specific definition that incorporates seismic hazard, exposure, and the vulnerability of the community and structures affected by shaking. In our use of the term “QRA”, we are actually assessing the conditional probability of fault slip, which is, strictly speaking, related only to seismic hazard.
The map in Figure 1 shows fault segments compiled from published sources and information donated by petroleum operators (Darold and Holland, 2015). The database contains 26,313 fault segments, each defined by two connected coordinate points. No information about fault dip or depth is included. The fault map is complete to varying degrees, as areas of tectonic uplift expose faults at the surface and sediments conceal them, unless known from wells or seismic data. This results in fewer mapped faults in the study areas considered here (as 2–3 km of sedimentary rock overlies crystalline basement), but the earthquakes in these areas are evidence that the lack of mapped faults is a reflection of undersampling, not the absence of faulting. The majority of earthquakes in north-central Oklahoma are at 5–6 km depth (McNamara et al., 2015), ∼2–4 km into crystalline basement where faults are difficult to image. As can be inferred from the distribution of M ≥2.9 focal plane mechanisms shown in the Figure 1 inset, the majority of the earthquakes are not associated with mapped faults. To calculate the probability of fault slip in a given stress field, we need information about both fault strike and dip, ideally its actual three-dimensional geometry. In the absence of dip information, we consider a probabilistic distribution of dips.
INVERTING MOMENT TENSORS FOR THE TECTONIC STRESS FIELD
We characterize the stress field in the study areas utilizing 335 focal plane mechanisms from the St. Louis University Moment Tensor program (Herrmann, 2016), determined by waveform modeling, and from the U.S. Geological Survey National Earthquake Information Center (NEIC) catalog (http://earthquake.usgs.gov/earthquakes/). Each focal mechanism provides the strike, dip, and rake of two nodal planes, one of which represents the slip plane and the other an auxiliary nodal plane. We constrain the tectonic stress orientations and relative magnitudes using a linearized inversion of focal plane mechanisms (after Michael, 1984). Ideally, one would invert at least 25–30 well-constrained focal mechanisms per area (Townend and Zoback, 2001). Thus, there is a trade-off between the size of the area selected and the number of focal mechanisms in each area that are available to constrain inversions.
Following Angelier (1990), we represent the relative magnitudes of the principal stresses with the parameter ϕ which is defined by the magnitude of the intermediate principal stress (S2) with respect to the maximum (S1) and minimum (S3) principal stresses:
When combined with the orientation of the principal stresses, the ϕ value describes the style of faulting in a given area. Table DR1 in the GSA Data Repository1 presents a summary of the stress inversion results. We also show in the Data Repository the convergence and decrease of uncertainty of the stress orientation and ϕ values for each of the six study areas as focal mechanisms are added to the inversion. Consistent values of stress orientation and relative magnitude are obtained in each area, and there is typically <5° between the stress orientation obtained by the inversion of the earthquake moment tensors and independent wellbore stress information in each study area (Alt and Zoback, 2014; Heidbach et al., 2010). The stress field in area 2 appears to be more complex than those of the other areas, and therefore area 2 was split into a northern (area 2N) and southern area (area 2S), each with a comparable number of moment tensors. In area 2S, the focal mechanism inversions do not match the orientation of maximum horizontal stress from wellbores (Fig. DR4 in the Data Repository) nor the focal mechanism inversions in the surrounding areas. This discrepancy may reflect that there are not a sufficient number of focal mechanisms to satisfactorily constrain the stress in area 2S or that there are localized variations of the stress field. Whatever the reason, we do not have confidence that the stress field is sufficiently constrained in area 2S to analyze it.
While the orientation of maximum horizontal stress, SHmax, throughout central Oklahoma is fairly uniform (Fig. 1), relative stress magnitudes vary. Study areas 3, 4, 5, and 6 are characterized by strike-slip faulting (SHmax > SVertical > Shmin; hmin is minimum horizontal), areas 1 and 2S are characterized by both strike-slip and normal faulting (SHmax ≈ SVertical > Shmin), and area 2N, much like southern Kansas (Rubinstein et al., 2015), is characterized predominantly by normal faulting (SVertical > SHmax >Shmin). We bootstrap the moment tensor inversion to constrain the uncertainty in stress orientation and ϕ value that we input into the QRA.
THE PROBABILITY OF TRIGGERED FAULT SLIP
In the context of Coulomb faulting theory, determining whether a fault will slip in response to fluid injection depends on the magnitude and orientation of the stress field, the orientation of the fault, pore pressure, and material parameters such as the coefficient of friction. Because there is uncertainty in each of these parameters, we evaluate the possibility that the mapped faults might slip in response to pore pressure increases using QRA, a Monte Carlo method used to evaluate the probability of an uncertain outcome incorporating uncertainty in the input parameters. In applied rock mechanics, QRA has been used to evaluate wellbore stability (Moos et al., 2003) and fault seal (Jones and Hillis, 2003). Leakage of CO2 from carbon sequestration reservoirs through wells, faults, and fractured caprock was investigated using QRA by Chae and Lee (2015). Chiaramonte et al. (2007) used QRA to evaluate if increases in pore pressure from a CO2 injection pilot project might induce fault slip. The Chiaramonte et al. (2007) analysis assumed either purely strike-slip or normal faulting and only considered uncertainty in fault orientation. In this study, we generalize the QRA to include potential slip in any direction on mapped faults from uncertain stresses. Use of QRA is essential as there is uncertainty in all of the model parameters and no information available on mapped fault dip.
Figure 2 shows the distributions of ϕ, friction, pore pressure, and stress used in the QRA of study area 6. Similar figures for the other study areas are included in the Data Repository, as well as how the uncertainty in each parameter is established and a table of values. Figures 2C and 2D also show response surfaces in red, which use the most likely value in each distribution in Figure 2 (indicated by the vertical dotted black lines in the other distributions) to show the required pore pressure to induce fault slip based on a fault’s dip (Fig. 2C) or strike (Fig. 2D). The black horizontal line in Figures 2C and 2D shows the 2 MPa expected pore pressure perturbation. This magnitude is based on Nelson et al. (2015) and the fact that wellhead pressures remain subhydrostatic after injection stops. Note that few fault strikes are present in study area 6 that would be expected to slip in response to a 2 MPa pore pressure change.
RESULTS OF THE QRA
In Figure 3, we apply QRA to the mapped faults in study area 6 utilizing the uncertainty distributions shown in Figure 2. As with the bootstrap of the moment tensor inversions, we evaluate 10,000 random combinations of parameters for each mapped fault segment to evaluate the conditional probability of slip as a function of pore pressure perturbation given the model assumptions described above in the context of Mohr-Coulomb faulting theory. We present the result as an empirical cumulative distribution function (CDF) showing the probability of slip on a known fault as a function of pore pressure increase. In Figure 3A, we map 484 fault segments in area 6 colored by their corresponding curves in Figure 3B. Figure 3B shows the probability of fault slip as a function of pore pressure change. We use a traffic-light color scale where fault segments with <1% conditional probability of slip in response to a 2 MPa pressure perturbation are colored green. Those with a >33% conditional probability are colored red. The large north-northeast–striking Wilzetta fault is colored green in Figure 3 as it strikes at too high an angle to the direction of maximum horizontal compression to be activated by a relatively small pore pressure perturbation. The Wilzetta splay fault just west of Prague (source of the 2011 M 5.6 earthquake) is colored red and would have had nearly a 50% probability of slip in response to a 2 MPa pressure perturbation. Note that there are other mapped potentially active faults. Maps and CDF curves similar to those in Figure 3 for each of the other study areas are presented in the Data Repository.
Figure 4 maps the probability of fault slip in response to a 2 MPa pore pressure change in all six study areas and an indication of the style of faulting (strike-slip, strike-slip–normal, or normal faulting) based on the focal mechanism inversions as summarized in Table DR1. We observe that the majority of mapped faults are not likely to be activated by modest pore pressure changes. Figure 4 also maps M 3+ earthquakes and saltwater disposal wells that injected more than 300,000 barrels in any month from 2009 through 2014. The 13 February 2016 M 5.1 earthquake near Fairview, Oklahoma, is circled in the southwestern part of area 1. The focal mechanism of the earthquake indicates a steeply dipping, northeast-striking fault plane that aligns with similarly striking faults mapped to both the southwest and northeast of the epicenter, colored varying shades of yellow.
The results presented in Figures 3 and 4 show the conditional probability of fault slip, given the uncertainties in the model parameters. The conditional probability of fault slip differs from the probability of an earthquake on a given fault because we don’t know where the fault is in the earthquake cycle. That is, we assume that if a fault is optimally oriented to slip in response to a 2 MPa pore pressure perturbation, it is likely to. However, had there been an earthquake on the fault in the past few thousand years, the stress drop in that earthquake might not have recovered, reducing the probability of fault slip. Furthermore, our methodology considers the probability of fault slip, not the size of any possible earthquake. In other words, while large faults can produce larger earthquakes, the potential for slip over part (or all) of a fault is unknowable. The San Andreas fault in California (USA) produces approximately one million M ∼2 events for every M ∼8 event. The fact that the majority of earthquakes are occurring on unmapped faults limits the applicability of our analysis to known faults and requires constraining other parameters such as the stress and pore pressure. We ignore poroelastic effects as these are likely to be minor in response to such small pore pressure changes. Similarly, while the Prague main shock may have been caused by a combination of static stress transfer from the M 5 foreshock and hydraulic triggering (Keranen et al., 2013), our analysis identifies that this fault could be triggered by a small perturbation.
We have taken a conservative approach to evaluating the potential for mapped fault slip in a Mohr-Coulomb framework. We also approach the magnitude of pore pressure changes in a conservative manner. In other words, assuming that an ∼2 MPa pore pressure change is being transmitted throughout the region of seismicity and to the ∼5–6 km depth of the earthquakes, we are considering a worst-case hydrologic scenario. It would be preferable to use constrained hydrologic modeling to assess the pore pressure change resulting from injection at the fault. As discussed by Keranen et al. (2014), this is difficult because many of the key hydrologic parameters are unknown. For example, if one were to assume injection is into a relatively isotropic, permeable medium, pore pressure would spread out uniformly from an injection well, and the pressure change at some distance could be estimated in a straightforward manner. However, permeability of faulted basement is likely to be complex. For example, some faults could be barriers to cross-fault flow but could have highly permeable damage zones that allow for rapid fluid pressure propagation parallel to the fault plane (e.g., Hennings et al., 2012). This complexity could result in relatively large pressure changes in some places at significant distance from injection wells (compared to what would be predicted from an isotropic permeability model) but smaller pressure changes in others.
We present a framework for calculating conditional probability of fault slip from a pore pressure perturbation by modeling Mohr-Coulomb slip, incorporating the uncertainties in relevant parameters. The application to induced seismicity in north-central Oklahoma confirms the potential for slip on the northeast-striking splay fault that produced the 2011 Prague earthquake but calculates a low probability of induced slip on several other large faults in the region (such as the Nemaha fault) in response to the relatively small pore pressure changes. The 2016 Fairview earthquake occurred on the previously unmapped extension of a fault mapped as having moderate but varying probability of slip. Had the 289°-striking nodal plane of the focal mechanism of the M 5.8 Pawnee event been mapped, it would have been identified as a potentially active (red) fault.
The stress orientations obtained by inverting the focal mechanisms generally compare well with stress determinations made in wells in the same areas, building confidence in the inversions. The variation in stress has a significant effect on the orientation of potentially active faults within each area. This approach requires mapped faults, but the majority of the earthquakes in Oklahoma are on unmapped faults. It is also necessary to constrain the geomechanical parameters to assess the probability of injection-induced fault slip, and additional mechanisms could be considered. Despite these obstacles, this approach provides a rigorous quantification of mapped fault slip potential incorporating uncertainties in relevant parameters.
The Stanford Center for Induced and Triggered Seismicity funded this work. We thank Stanford Professors Greg Beroza, Bill Ellsworth, and Jack Baker for helpful conversations, and Hiroki Sone for coding the Michael (1984) inversion and the resolving stress on fault codes. We also thank Peter Hennings, Cliff Frohlich, and one anonymous reviewer.
- Received 21 June 2016.
- Revision received 14 September 2016.
- Accepted 16 September 2016.
- © 2016 Geological Society of America