Stability of the Fault Systems That Host‐Induced Earthquakes in the Delaware Basin of West Texas and Southeast New Mexico

The Delaware basin of West Texas and Southeast New Mexico is one of the world’s most productive petroleum basins (Fig. 1). Approximately 16,000 horizontal wells have been drilled and hydraulically fractured in Permian, shale‐dominated reservoirs (IHS Markit, 2021). In 2020, the basin produced >700 million bbls of oil and >3.1 trillion cubic feet of natural gas (IHS Markit, 2021; Energy Information Administration [EIA], 2021). Since 2010, this unconventional development has required the subsurface disposal of 8.6 billion bbls of wastewater (saltwater disposal [SWD]) into strata both above and below petroleum‐producing intervals (Lemons et al., 2019), and there are currently ~1200 active SWD wells in the basin (IHS Markit, 2021; Fig. 1). The basin has also experienced a significant number of earthquakes that are spatiotemporally linked to unconventional oil and gas development. Frohlich et al. (2016) discussed the history of induced seismicity in the region, and Frohlich et al. (2020) demonstrated that induced seismicity in the Delaware basin initiated in 2009 and accelerated significantly in 2016. From January 2017 through June 2021, the TexNet Earthquake Catalog (Savvaidis et al., 2019) and New Mexico Tech Seismological Observatory (NMT) have cataloged >3500 $ML$ 2.0+ earthquakes in the basin. Lomax and Savvaidis (2019), Savvaidis et al. (2020), and Skoumal et al. (2020) provided spatiotemporal causal assessments, each proposing complex linkages of earthquakes to both hydraulic fracturing and SWD. Skoumal et al. (2020), Tung et al. (2020), and Gao et al. (2020) provided analyses linking new earthquake sequences in the northern Delaware basin to SWD in deep strata regionally. Dvory and Zoback (2021a) show that pressure depletion from prior production from the Delaware Mountain Group and Bone Spring Formation limits the likelihood of shallow seismicity in some areas of the basin. Zhai et al. (2021) use geomechanical modeling to propose that shallow SWD is the primary driver of earthquake inducement in the basin. Skoumal and Trugman (2021) expand the earthquake analysis of Frohlich et al. (2020) and conclude that shallow SWD in the southcentral Delaware is the primary driver of earthquake inducement.

Human influence leading to earthquake inducement occurs within a complex geomechanical system. With stratified pore pressure (Pp) regimes, multiple faults systems, regionally heterogeneous tectonic stressing, and highly complex operational history, the Delaware basin is the most challenging region impacted by induced seismicity that has been studied to date. Given the accelerating earthquake rate in the basin, it is vital to develop causal assessments that are grounded within the complex geological architecture of the basin so that mitigation strategies can be implemented that target specific agents of causation. To advance this objective, we present comprehensive new fault interpretations, deterministic analysis of fault‐slip potential (DFSP), and the association to earthquake sequences. Our DFSP maps of multiple extant fault systems extend throughout the region of the Delaware basin and Central basin platform. They provide a resource not only for the understanding of induced seismicity but also for assessments of fault stability in application to landform evolution, for hazard assessment for sequestration, and for hydrocarbon production mechanics.

The Delaware basin is a subbasin of the Permian basin province of West Texas and Southeast New Mexico. The Delaware basin and Central basin platform are defined structurally by a complex network of BR faults (Fig. 1). The basin and its flanks have been affected by several extensional and contractional tectonic events cumulatively leading to the present‐day tectonic fabric (Horne et al., 2021 and references therein). The basin formed in Mississippian through Permian time in the autochthon of the Ouachita–Marathon orogeny and was defined structurally as a foreland basin, as it became downthrown to the west‐vergent, BR, contractile Central basin platform that formed as a prominent constituent of the Ancestral Rocky Mountains.

For this analysis, we refer to the following earthquake catalogs to cover the entirety of the Delaware basin (Fig. 1). For New Mexico, we use NMT. For Texas, we use TexNet for 2019 through June 2021 (Savvaidis et al., 2019) and Lomax and Savvaidis (2019), who provided advanced relocation of TexNet for 2017–2018 (L&S). The TexNet and L&S catalogs have an $Mc=1.5$ (Lomax and Savvaidis, 2019; Savvaidis et al., 2019). There are no duplicate earthquakes from these combined catalogs. To assist with fault mapping and earthquake‐to‐fault association, we use P. Li, and A. Savvaidis (unpublished manuscript, see Data and Resources) who provide relative relocations for TexNet and L&S (TNr, Fig. 2a).

Dvory and Zoback (2021a) show from both direct measurements and calculations that are based on Mohr–Coulomb criteria that the crust is in critically stress state in the seismically active area of the Delaware basin. Lund Snee and Zoback (2018, 2020) compiled various types of data in the region to define the maximum principal horizontal stress azimuth (⁠$SHmax$ Az) and stress ratio (⁠$Aϕ$⁠). They show that $SHmax$ Az varies systematically across the region, smoothly transitioning from northwest–southeast in the southern Delaware basin to east–west in the center of the basin to almost north–south in its northern areas. The interpretation of $Aϕ$ varies from about 0.6 (normal faulting) in the western boundary of the basin to ~0.85 in the Central basin platform and to ~1.0 in the Midland basin to the east, indicating a transition to combined strike‐slip and normal faulting. Fault‐plane solutions for earthquakes in the Delaware basin indicate normal faulting (TexNet) with planes striking parallel to the local direction of $SHmax$⁠, as expected. Dvory and Zoback (2021b) provide a smoothed and interpolated stress field (Fig. 1) that closely fits the discrete stress observations. We use this interpolation for our DFSP analysis and assume an intermediate $Aϕ$ value between the Delaware basin and the Central basin platform of 0.7.

Many faults deform the Delaware basin and the Central basin platform; the most significant are the BR reverse faults that formed in late Paleozoic time and control its primary structural architecture (Fig. 1). Modern rupture along these faults must, therefore, represent normal slip reactivation. Horne et al. (2021) provided a new and comprehensive 3D interpretation of these faults from well‐based framework mapping, 3D seismic data, and prior publications. Horne et al. (2021) classify their faults in terms of mapping confidence, and we use these faults for our DFSP analysis.

The primary BR fault fabric strikes north‐northwest, is dominantly reverse, and forms the first‐order structural relief in the region. This primary fabric is most intensely developed along the western margin of the Central basin platform. A secondary fault fabric strikes west‐southwest to west‐northwest and is dominantly associated with reverse and strike‐slip offset. This secondary fabric occurs in concentrated zones that are distributed from south to north along the overall strike of the basin (e.g., Grisham fault zone, Culberson‐Mentone earthquake zone). Complex deformation occurs at the points of convergence between the primary and secondary fault fabrics. The BR faults cut up section from basement to levels as shallow as the Wolfcamp Formation (Fig. 2b). The total fault‐trace length (length) of BR faults is ~6500 km, with individual segments ranging in length from 5 to >100 km. The BR faults have throws from 50 to >1000 m and a mean throw to length ratio of ~1:25. In the regional‐scale 3D fault framework in Horne et al. (2021); the mean surface area of the BR faults is $58km2$ with $n=638$⁠. This is an underestimate, because the interpreted faults are not interpreted deeply into basement.

The great majority of recent earthquakes do not occur on BR faults, but on northwest‐trending shallower normal (SN) faults in the central Delaware basin that do not extend to the depth of the BR faults. Our interpretation of these seismogenic faults is newly presented here and based on integrating 3D reflection seismic data (extending the work of Charzynski et al. 2019; Cook et al., 2019), unpublished records from Railroad Commission of Texas, and as directly reported to us by petroleum operators (Fig. 2a). Many of the SN faults follow linear patterns of Interferometric Synthetic Aperture Radar (InSAR) surface deformation (Staniewicz et al., 2020), and the TNr earthquakes closely follow many of the mapped SN faults. Therefore, both the InSAR lineations and TNr earthquakes inform our SN fault interpretation. We are confident in combining these data to map SN faults, because all these indicators come together in areas where we can independently confirm the interpretation using 3D reflection seismic data (e.g., data region 3Dc, Fig. 2a). We employ the same high‐ and moderate‐confidence mapping criteria for the SN faults as Horne et al. (2020).

Throughout the seismically active area, the SN faults strike parallel to $SHmax$ azimuth, varying smoothly from northwest‐striking in the south to west‐northwest‐striking in the north. The faults often occur in pairs forming narrow graben bounded by steeply dipping faults and have a mean throw to length ratio of ~1:100. They appear stratigraphically and mechanically bound to the Delaware Mountain Group, Bone Spring Formation, and uppermost Wolfcamp Formation. The SN faults occur primarily in Reeves County but extend into adjacent Texas counties including Pecos, Ward, and Culberson. Where we have 3D seismic control, the SN faults have a mean surface area of $4.2km2$ with $n=41$⁠. The total fault‐trace length of SN faults is ~1450 km, with individual segments ranging in length from 0.5 to 20 km. Of the ~1450 km of mapped SN faults, ~550 km (38%) of fault length is considered to be high‐confidence interpretations; the remaining 900 km (62%) of segment length is interpreted as moderately confident. The estimated ~1450 km of SN fault‐trace length is a minimum, as there is strong anecdotal evidence of many more SN faults within Delaware basin to the northwest and south of the distribution we show (Fig. 2a) but where data control for mapping is not available to us.

There was no previously recognized surface expression of these faults until the observation of recent InSAR displacements (Staniewicz et al., 2020). Anderson (1981) found that some SN faults in the region resulted from organized flowage and dissolution of evaporite‐dominated strata above the Delaware Mountain Group. Cook et al. (2019) used analysis of 3D seismic anisotropy to map the fault zones from the Delaware Mountain Group into the Wolfcamp Formation. Charzynski et al. (2019) concluded that in some areas, the SN faults extend downward as permeable fracture zones that cause wells targeting the Wolfcamp Formation to have diminished hydrocarbon productivity, and produce anomalous volumes of $H2S$ and water.

We assess the slip potential of the BR and SN faults in our interpretations using a deterministic approach (DFSP), which assumes that the Coulomb failure criterion is applicable and describes the Pp increase (⁠$ΔPp$⁠) needed to achieve criticality. Critical $ΔPp$ is the perturbation in Pp from ambient, therefore low values of $ΔPp$ or DFSP indicate a greater sensitivity (Zoback, 2007). Other similar works assessing the slip potential of induced earthquakes include Walsh and Zoback (2016) and Lund Snee and Zoback (2018). These works employed probabilistic approaches (PFSP) to investigate the slip potential of fault systems assuming hypothetical yet plausible global increases in $ΔPp$⁠. The PFSP analysis of Hennings et al. (2019) compared unperturbed versus hypothetical yet plausible global increases in $ΔPp$⁠. Hennings et al. (2021) assessed the evolution in PFSP using a deterministic regional model of temporally varying $ΔPp$⁠. Morris et al. (2021) also employed a deterministic assessment of slip potential for 3D fault surfaces in the Delaware basin region using the Horne et al. (2021) interpretation but focused exclusively on the BR faults.

The DFSP analysis uses fault strike and dip, the depth of the trace that samples the fault surface, the in situ stress tensor, and scalar geomechanical parameters such as Pp and the coefficient of fault friction. Based on the work of Luo et al. (1994), we employ a hydrostatic condition of 0.01 MPa/m, which is representative of the native pressure state of the Delaware Mountain and Ellenburger Groups prior to widespread petroleum operations. Following Dvory and Zoback (2021a), we assume acritical fault friction of $μ=0.7$⁠. The strike and dip of the BR faults come from the upthrown intersection of the faults and the Ellenburger Group in Horne et al. (2021). Thus, our DFSP calculations represent the slip potential of BR faults that cut from basement through to levels as shallow as the Wolfcamp Formation, the uppermost of which is the current primary target for petroleum production using hydraulic fracturing. The strike and dip of the SN faults come from the intersection of the faults and the Delaware Mountain Group. Fault dip was determined explicitly for the faults constrained by 3D seismic data. We use the mean of the explicit dip (72°) for SN faults that lack 3D seismic control. The traces for both BR and SN fault systems are sampled at 1 km increments laterally for the DFSP analysis.

The DFSP results are shown in the map in Figure 3 and as distributions in Figure 4. About 62% of BR fault length have a very high DFSP of >5.0 MPa (implying stability). Only 18% are considerably less stable with DFSP ≤2 MPa. The most sensitive BR faults occur primarily along the southwest flank of the Central basin platform, the southeast axis of the Delaware basin, and along segments of the secondary fault fabric in the central Delaware basin (Fig. 3). The DFSP for the SN faults is very different—71% (~1030 km) have a DFSP of ≤2.5 MPa indicating that they are prone to reactivation in response to modest $ΔPp$⁠. Nearly all of the SN faults are inherently sensitive to slip.

Only a few of the earthquake sequences potentially coincide spatially with BR faults (Fig. 3). The most noble of these is the Culberson–Mentone earthquake zone, where there are numerous highly unstable BR faults associated with recent earthquake sequences. Other areas include in Lea County, New Mexico; along the Texas/New Mexico border; and in the Waha area. The SN faults, which are uniformly unstable, spatially correlate to many earthquake sequences in the central Delaware basin.

Combining the interpretation of both the fault systems—the DFSP and the earthquake epicenters—we conclude that rupture has occurred in distinct geographic groupings. The majority of the earthquakes that have occurred along deeper BR faults are concentrated in the western part of the basin, especially in the Culberson–Mentone earthquake zone. We concur with the previous works of Gao et al. (2020), Savvaidis et al. (2020), Skoumal et al. (2020), Tung et al. (2020), and Zhai et al. (2021) that the most plausible causal agent for earthquake inducement in this area is from $ΔPp$ from SWD into strata above basement—principally units of Devonian. Our DFSP analysis indicates that rupture of many faults in this area can occur with very small $ΔPp$⁠.

We show here that induced earthquakes in the southcentral part of the basin have occurred mainly along the SN faults, which, as a system, have uniformly low DFSP and are highly sensitive, requiring ≤2.5 MPa $ΔPp$ to achieve criticality. The faults closely follow the azimuth of $SHmax$ (Fig. 1). Ge et al. (2020) show that SWD into the Delaware Mountain Group has caused $ΔPp≥2MPa$ in the same areas as these active faults. We, therefore, concur with Skoumal and Trugman (2021) that $ΔPp$ from SWD into the Delaware Mountain Group is a primary causative agent for fault rupture and seismicity in this area. Savvaidis et al. (2020) demonstrate that hydraulic fracturing is an additional causative agent for more isolated earthquake sequences in this area. With SWD into the Delaware Mountain Group and hydraulic fracturing of the underlying shale intervals, it is, therefore, likely that the induced seismogenic ruptures are occurring at depths shallower than basement or the strata immediately above.

Although the velocity structure of the Delaware basin is highly complex, TexNet has relatively sparse station spacing and uses a coarse, regional‐scale, 1D velocity model for locating hypocenters. Therefore, the reported uncertainty (e.g., one standard deviation [st. dev.]) in depth provided by TexNet should be used to identify the comprehensive uncertainty that is much larger (three st. dev. accounts for 99.7% of the depth estimates for each earthquake). The hypothesis of shallow rather than deep earthquakes in the central Delaware basin is supported by the high‐resolution earthquake depth mapping by Sheng et al. (2020). They show that earthquakes in their study area (see Fig. 2a) are principally concentrated in deeper levels of the Delaware Mountain Group and upper levels of the Bone Spring Formation rather than in units under the Wolfcamp Formation and downward into basement. Combining these factors, along with their inherent instability, we conclude that most of the earthquakes in the central Delaware basin occur along the SN faults as we have mapped them. Many of these faults are now neotectonic. We concur with Zhai et al. (2021) that SWD into the Delaware Mountain Group provides the causative stress change required for fault rupture but, rather than stress being transmitted downward 3–4 km causing rupture of exclusively deep faults, it is the SN faults in direct communication with $ΔPp$ that host most of the earthquakes in the central Delaware basin.

In Figure 5, we show the magnitude histories of earthquakes from the Culberson–Mentone earthquake zone, which has an increasing maximum magnitude with time, and the southcentral Delaware basin earthquake zone, which has a steady maximum magnitude with time. Since 1 January 2018, the Culberson–Mentone earthquake zone, where the earthquakes are occurring along BR faults, has had >1300 earthquakes of $ML≥2$ and 19 with $ML≥3.5$⁠, including the 26 March 2020 $ML$ 4.6. Conversely, in the entire footprint of the SN faults since 1 January 2017, there have been >1000 earthquakes with $ML≥2$ and only one $ML≥3.5$⁠, the 22 December 2018 $ML$ 3.8. As described previously, only a few earthquake sequences have occurred on the BR faults. This may be beneficial, given that BR faults can have large surface areas and cut relatively stiff country rock (e.g., basement and platform carbonates), which implies the potential for larger seismic moments during slip (Zoback and Gorelick, 2012). In contrast, the SN faults are strata‐bound, have smaller surface areas, and cut rocks with a relatively lower stiffness (e.g., sandstones and shales), implying a reduced potential for larger seismic moments during slip. The earthquake history in Figure 5 reinforces this hypothesis. With this hypothesis as a prompt, we encourage additional quantitative study comparing the earthquake data from these two regions.

The DFSP map we provide should be used to assess the fault‐slip hazard related to petroleum operations in the Delaware basin region and as a general guide for mitigation. It can also be used for hazard assessment associated with sequestration. However, our characterization is based on the strength of the available data, and our fault interpretations should be considered as inherently incomplete. Hazard assessment and mitigation at the local scale should be performed using the best data available for fault and geomechanical parameterization.

In their natural state, the faults in both the systems are relatively stable, which explains why the earthquake rate was low prior to the initiation of unconventional petroleum operations.

BR faults that are potentially unstable exist throughout the Delaware basin region, but only a few earthquake sequences are associated with these faults, most prominently the Culberson–Mentone earthquake zone. These faults typically have large surface areas and cut relatively stiff strata; therefore, their maximum earthquake moment can be large.

SN faults are numerous in the central Delaware basin; they are highly sensitive to reactivation, they spatially correlate with many recent earthquake sequences, and they have recently become neotectonically active. About 71% of their length becomes unstable at Pp changes of ≥2.5 MPa—a level of $ΔPp$ that has been shown to have been created in the region due to SWD into strata cut by these faults. Perturbance from hydraulic fracturing operations is also implicated in causing earthquakes along these faults. Given the mechanically stratified architecture of these faults and the stiffness of the rocks they cut, maximum earthquake moment may have an upper limit, perhaps to what has been observed thus far.

All existing data we use can be accessed at the respective journal cited. The data on faults and fault‐slip potential we discuss are available for download at the Texas Data Repository at doi: 10.18738/T8/TBTRXM. FSP software can be accessed at https://scits.stanford.edu/software (last accessed June 2021). The unpublished manuscript by P. Li, and A. Savvaidis (in revision), “Cross‐correlation relocation to identify active faults in the Permian Basin,” submitted to Seismol. Res. Lett.

The authors acknowledge that none have conflict of interest

The authors thank TGS for access to 3D seismic data and ExxonMobil for data that assisted with fault mapping. Funding support for Hennings, Horne, Li, and Savvaidis was provided by State of Texas through the TexNet Earthquake Monitoring Program and the sponsors of the Center for Integrated Seismicity Research (CISR) at The University of Texas. The sponsors of the Stanford Center for Induced and Triggered Seismicity (SCITS) provided funding support for Dvory and Zoback. The authors thank two anonymous reviewers for reviews that greatly improved this article.