Hydraulic Fracturing, Fault System Architecture, and the Details of Anthropogenic Earthquakes in the Post-Pennsylvanian Delaware Basin of West Texas Open Access

In this article, we address the question: what, if any, patterns emerge when induced earthquake magnitudes are treated as proxies for displacement along a fault system and over time? We focus on two fault systems within a larger area of interest (AOI) of approximately 8630 km² in the southern Delaware Basin around and south of the town of Pecos (Figures 1(a) and 1(b)). The rapid increase in earthquake frequency in the west Texas Delaware Basin was coeval with the increase in oil and gas production activities in the region [1] (Figure 1(c)). Most of the seismicity in the AOI is confined to the post-Pennsylvanian sedimentary fill of the basin [2], on what we term here shallow normal faults (SNF; Figure 1(a)). Although the current stress regime in both the basement and sedimentary fill of the AOI is primarily extensional [3], SNF is disconnected from the deeper, basement-rooted faulting, which dates from pre-Pennsylvanian, predominantly contractional deformation in the region.

Hydrocarbon production activities in the form of hydraulically fractured (HF) horizontal wells are concentrated at the base of the Bone Spring and top of the Wolfcamp formations, saltwater disposal (SWD) via injection wells is confined to the Delaware Mountain Group (DMG) (Figure 1(a)). Fault systems rooted in the Proterozoic, sub-Ellenburger crystalline basement, and that penetrate upward into the Pennsylvanian strata of the basin fill, are known to exist within the AOI [4-6]. These predominantly contractional and strike-slip systems have been reactivated since about 2009; however, the contribution to the cumulative seismicity from these basement-rooted faults is minor as compared with the predominantly extensional faults (SNF) in the overlying sedimentary basin fill (Figure 1(d)).

Detailed hydrogeological modeling of wastewater injection (SWD) into the DMG indicates an AOI-wide, monotonic increase in pore pressures over the period 2010 to 2021 [7], with 54% of the model experiencing pore pressure increases >50 psi (0.69 MPa) and 12% experiencing >200 psi (1.38 MPa). By the end of 2021, 90% of the earthquakes that had occurred in the AOI were located within the model area experiencing pore pressures of >3000 psi (20.68 MPa; hydrological-unit-thickness-weighted average; equivalent to a pore pressure gradient of 0.51 psi/ft or 0.011 MPa/m; Figure 2(a)) to depths of ~2500 m, within the DMG. Thus, regional scale pore pressure increases generated by SWD can contribute to the destabilization of preexisting fault systems in the ambient stress state of the Delaware Basin [5].

HF can also trigger seismicity, depending on the geological conditions [8, 9], although these effects are likely more localized. It is difficult to pinpoint specific HF completions as likely causes of specific earthquakes; however, there is a very strong AOI-wide similarity in the cumulative record of earthquake occurrences and HF completions (Figure 1(c)).

The pattern of seismicity in the AOI is characterized by curvilinear zones of high earthquake frequency that strongly resemble fault systems (Figures 2(a) and 2(b)). These zones also correlate closely with lineaments expressed in the vertical displacement of the Interferometric Synthetic Aperture Radar (InSAR)-defined topographic surface over the period of January 2015–December 2021 [2] (Figure 2(c)). The lineaments are attributed to anthropogenic reactivation of SNF mapped by Horne et al. [5] and examined by Hennings et al. [4]. Analyses of slip on shallow faults in the region of Pecos, TX, were published, and the inversion of InSAR-derived surface displacement to constrain the depth and magnitude of fault slip is robust [10-12].

We chose two well-defined fault systems to be the subjects of the work described here (Figures 1 and 2). Fault systems 1 and 2 are similar in size (12, 14 km in length) and have similar InSAR signatures (Figure 2(c)). However, they have experienced different pore pressure perturbation histories resulting from SWD and HF activities in their immediate vicinities. Using local magnitudes (M_L) for events from two earthquake catalogs and assuming that these are a reasonable proxy for fault displacement, we calculate along-strike displacement histories for the two fault systems. The strong resemblance between these displacement profiles and those measured in natural and experimentally formed normal fault systems [13] leads us to suggest that the spatial distribution (especially along-strike) of slip events triggered by pore pressure perturbations is determined by the preexisting fault system architecture. This could provide a predictive tool for induced seismicity if the fault system architecture is known prior to the initiation of subsurface fluid injection activities.

The work described here was supported by the Center for Injection and Seismicity Research (CISR) at the University of Texas Bureau of Economic Geology. CISR collects, organizes, and analyzes data (specifically production and injection well data, InSAR data, and commercially collected seismicity catalogs) from a wide range of contributors.

Production well locations and completion dates were provided by a well database compiled by CISR using the commercial S&P Marketplace well drilling database; for the areas of interest, this dataset spans the period of April 22, 2006, to October 29, 2022. Completion dates are likely accurate within approximately 7 days, subject to reporting protocols of the numerous operators in the Basin. Locations of completions are given as the latitude, longitude (or easting and northing), and the total vertical depth of the midpoint of the well’s lateral leg, corrected for surface elevation. Locations and injection histories of SWD wells were obtained from the B3 insight injection well database for the time period January 2010 to December 2023.

Pore pressure perturbations generated by SWD are derived from the hydrogeological model of the DMG developed by CISR [7]. Horizontal spatial resolution of this model is ~1.6 km, and vertical resolution is ~75 m. Here, we use annually extracted pore pressure values from December 31, 2010, to December 31, 2021, to establish a “background” variation in pore pressure.

TexNet [14] has collected and analyzed earthquake data from the Delaware Basin since January 1, 2017, and proprietary earthquake catalogs with higher resolution as provided for our use by CISR sponsors have been available since May 2019. Here, we use the TexNet catalog from January 2017 to May 2019 and a third-party commercial catalog thereafter. Spatial resolution of earthquake hypocentral locations is: vertical = 0.2 – 1.6 km and horizontal = 0.75 km. In the areas of interest analyzed in this work, it is recognized that TexNet hypocentral locations are typically calculated to be deeper than their likely origin [2]. For this reason, plots involving the vertical locations of earthquakes are constructed using only the third-party commercial catalog, although where only horizontal locations are required, both TexNet and third-party catalogs are used. Earthquake data collected after May 2019, and used here, only includes events with magnitudes >1.5. Prior to May 2019, events used in this work include a small proportion (≤12%) of events with magnitudes <1.5. Focal mechanism interpretations for ~100 third-party catalog events are also available within the AOI, and these are used to support interpretations of fault system architecture.

InSAR data collected over the period 2015–2022 [4] provide a measure of ground deformation with high vertical resolution (centimeters), 160 m horizontal resolution, and ~1 year temporal resolution.

There is no strict definition of the term “fault system,” so for clarity, we define our use of the term here. A fault system is a series of faults formed by or responding to a stress state and acting in concert with each other. In our definition, we place no constraints on size—many characteristics of fault systems are independent of scale. The two fault systems analyzed in detail in this article are of similar size and themselves consist of smaller systems and exist within larger systems.

In this work, we refer to “shallow” earthquakes and faults, and we define the term shallow for faults that are typically shallower than Pennsylvanian strata in the Delaware Basin, are often strata bound within Guadalupian or Leonardian (Permian) strata, and appear to be kinematically disconnected from basement-rooted faults in the area of study (Figure 1(a)).

Horne et al. [6] provide a succinct summary of the geologic setting of the Delaware Basin. The present-day structure of the Permian Basin, of which the Delaware Basin is the western sub-basin, is dominated by the results of interactions between the Ancestral Rocky Mountain Orogeny and Ouachita Marathon Orogeny during Late Mississippian to Early Permian [15, 16]. The Ancestral Rocky Mountain Orogeny was responsible for the uplift of significant basement-cored contractional features including the Central Basin Platform, which separates the Delaware Basin from the Midland Basin to the east [15-21], and the development of localized oblique and strike-slip structures [16, 17, 22-27]. The Ouachita Marathon Orogeny resulted from the collision of the Laurentian and Gondwanan plates [18, 28-36] causing the inversion of the Laurentian margin and the development of several flexural foreland basins. The Paleocene age Laramide Orogeny imposed primarily contractional and strike-slip deformation on the region [35, 37, 38], followed by Miocene and younger Basin-Range extension and Rio Grande Rift normal faulting [39-41]. These post-Permian extensional events control much of the present-day stress orientations in the southern Delaware Basin [42-45] and provide the context in which SNF have formed and have been reactivated.

Earthquakes occur on existing faults when conditions such as pore pressure perturbation, fault orientation with respect to ambient stress states, and fault cohesion and friction permit. Local magnitudes (M_L) are available for both the TexNet and third-party catalogs used here. Comparison of these two catalogs for events occurring in the AOI over the period May 2019 to June 2022 yields 715 events that can confidently be identified in both catalogs. A cross-plot for these events indicates that the M_L values reported for these catalogs are functionally the same (Figure 3(a)). Huang et al. [5] establish a correlation between the TexNet M_L and the moment magnitude (M_w) for an event, which they then use to calculate the seismic moment (M_o) for that event (after Hanks and Kanamori [46]). We follow a similar procedure using the reported M_L values from the TexNet and third-party catalogs to calculate both M_w and M_o for events in the fault systems analyzed here in detail. Both earthquake magnitude and seismic moment scale with displacement and slipped area during the event (e.g., Figure 3(b)) [47].

There can be considerable variation in interpreted displacement and rupture areas for a given magnitude (approximately one order of magnitude); however, for a set of geologic conditions (rock types, faulting stress regime, depth, and the concomitant stress-drop during an earthquake), these uncertainties will be decreased to a near linear (log–log) relationship [48, 49]. Given that geologic conditions are likely consistent for seismogenic slip along the SNF in the southern Delaware Basin, we treat earthquake local magnitude (M_L) as a semiquantitative proxy for fault displacement. However, for comparison, we also provide M_o values, calculated as described above and aggregated for time and distance along fault system strike. Both M_L and M_o are useful proxies for fault displacement, although neither provides a direct measure. M_o can be converted to an estimate of displacement if we know both the slipped area and strength characteristics of the fault and its host rock [48], in the absence of this information (as here) and used without conversion it amplifies the contributions of large events and subdues the contributions of small events. M_L is empirically linked to displacement, but by its nature as an exponent, it subdues the contributions of large events. We provide plots of both M_L and M_o as proxies for displacement, although the former is perhaps the more practical for the analysis pursued in this article because it provides a clearer picture of some of the subtleties of displacement distribution. Magnitude of completeness (M_c) for the induced seismicity in the Delaware Basin has been estimated to be M_c = 1.5 [8], and we mainly use events with magnitudes ≥1.5. We aggregate event magnitudes in 500 m sections along inferred fault systems, and these data are then further subdivided into 4-month time increments. The subdivided data can then be used to construct a graphical representation of the evolution of slip on a fault system through time.

Because the reactivated fault systems are shallow and predominantly normal displacement (fault system 1: 3–4 km depth; fault system 2: 2–3 km depth), there is a strong topographic surface response to subsurface faulting. Prior analysis of the comparison of the depth of fault slip and InSAR data collected over the time interval 2015 to the end of 2021 provides subcentimeter resolution of this surface deformation [4]. However, because fault slip is not the only contributor to ground surface deformation—others include local SWD, HF, oil and gas production activities, groundwater withdrawal, and possibly evaporite dissolution—this dataset is a secondary component of the analysis presented here.

Earthquakes generally occur on existing fault surfaces. Focal mechanism analysis provides two mutually consistent candidate surface orientations (the nodal planes) and their respective slip directions for each event. One of these nodal planes will be the reactivated, slipped surface, and the other (the auxiliary plane) will be perpendicular to both this plane and its slip direction. A range of fault surface orientations is likely to exist, and those with high slip tendency (T_s) orientations with respect to the ambient stress state are the most likely to be reactivated [50]. Slip direction on a reactivated fault surface is a function of the orientation of the surface with respect to the ambient stress state and the stress magnitude ratio, A_φ [51]. Consider two fault surfaces (strike/dip values, 180°/60° and 150°/60°) within an ambient normal faulting stress regime where S_Hmax (σ₂) is NS (S_Hmax azimuth = 180°) and S_hmin (σ₃) is EW and limiting friction is 0.6 (Figure 4, column A, rows 1 through 5). Both faults experience T_s > 0.5, but the slip rake on fault 150°/60° varies from −139° (A_φ = 1) to −90° (A_φ = 0); fault 180°/60° has a slip rake of −90° regardless of the value of A_φ. Expanding this to a fault population with high T_s orientations (0.5 < T_s < 0.6) and at intervals of 10° strike and dip, we can see a relationship between the range of slip rakes and A_φ (Figure 4 column B, rows 1 through 5). An efficient way of presenting this relationship is by plotting fault strike against rake for different values of A_φ (Figure 4 column C, rows 1 through 5). In these plots, we have provided more context by including surfaces with lower values of T_s (0.3-0.6). Assuming that those surfaces with the highest values of T_s are most likely to be represented in focal mechanism data, we further analyze the 0.5 < T_s < 0.6 data. The strike value at which the data distribution crosses the −90° rake line is the S_Hmax azimuth. The range of rakes varies with A_φ, from 0° to −180° (A_φ = 1) to 0° (A_φ = 0) (Figure 4 column C). Thus, these plots provide a means for estimating S_Hmax azimuth and A_φ for the ambient stress state. This analysis can also be used to assess the confidence in the stress estimation: The range of nodal plane orientations may be narrow or may not be symmetrically distributed about the S_Hmax azimuth [52]. This will be apparent from the strike versus rake cross-plot. Knowledge of friction and vertical stress then permits complete discovery of the stress tensor field.

It should be noted that the analysis presented in Figure 4 applies to actual slipped surfaces, whereas focal mechanism solutions provide an ambiguous result for each seismic event—two potential slipped surfaces but no indication as to which slipped and which is the auxiliary plane. Plotting strike versus rake for all nodal planes (half of which will be auxiliary planes) produces plots that are more complex than those in Figure 4 column C. However, because the strike and rake ranges of the auxiliary plane population are slightly smaller than those of slipped planes, the estimated values of S_Hmax azimuth and A_φ are essentially unaffected.

Two measures of fault stability are slip tendency (T_s) and critical pore pressure (ΔP_fc) [50, 53]. Using representative fault geometries and stress states [5, 54] for the two fault systems considered here, we calculate the likely values of slip tendency and critical pore pressure in the extant stress state.

T_s is the ratio of resolved effective shear to normal stress on a plane, and ΔP_fc is the pore pressure increase that will cause a fault facet to slip, i.e., reach limiting friction (maximum T_s), which we here assume to be 0.6. Larger values of T_s imply greater instability as do smaller values of ΔP_fc. T_s varies across a fault’s surface controlled in large part by variations in strike and dip of the fault’s surface. ΔP_fc varies with fault surface orientation, but also with depth (Figure 5). The relative stability of a fault in a stress state can be assessed by plotting cumulative area sorted by either T_s or ΔP_fc, such a plot provides information on how much of a fault’s area approaches or exceeds conditions for slip in a given stress state.

We use the location and timing of shallow (depth from ground surface, fault system 1: 2.25–4.25 km; fault system: 2–2.6 km) earthquake activity along two inferred fault systems south and east of Pecos, Texas (Figure 1(b)), to elucidate the history and spatial distribution of displacement accumulation along these systems for the period 2017–2023. Fault systems 1 and 2 are inferred from earthquake locations and InSAR data [4, 5]. Events in the vicinity of these fault systems that are confidently placed in the basement are rare and do not represent slip on the SNF that are the focus of this study. Oil and gas production activities are examined for possible links to the development of seismicity. Patterns of fault displacement are compared with ground-surface deformation, and we propose interpretations of fault system architecture that are consistent with published models of normal fault development. For consistency and context, we examine the inferred fault behavior in the light of ambient stress states as derived from focal mechanism analyses available from the third-party earthquake catalog and coeval with the seismogenic period 2017–2023.

All 17 time steps for fault system 1 are illustrated in online supplementary Appendix 1.

The first recorded and inferred displacement in the timeframe considered here (September 2017 to April 2023) occurred at the NW and SE ends of the fault system, with a less active section between 5 and 11 km measured along strike from the NW end of the system (Figure 6). This pattern is well established by late 2018, with three displacement maxima in the NW-most 4 km of the fault system, and two smaller maxima in the SE-most 3 km of the fault system (steeper sections of the 2018.83 highlighted curves in Figure 6(a); peaks in the lowest yellow highlighted curves in Figures 6(b) and 6(c)). The intervening 6 km contains four much smaller maxima.

The overall pattern remains stable as the whole fault system accumulates displacement (Figure 6). However, the displacement maximum near the SE end of the system (at ~12 km from the NW end) accumulates displacement faster than all other maxima but without significant lateral growth.

The overall rate of displacement accumulation decreases (see, e.g., the more closely spaced curves between 2020.5 and 2023.17; Figure 6(a)). However, the three, initially separate maxima between 2 and 5 km from the NW end of the system merge laterally to form a 3.5 km long displacement maximum. This contrasts with the other displacement maxima, which are on the order of 1–2 km long.

At locations 0.75 and 12.25 km along the fault system, displacement accumulation is nearly continuous, with short periods of stasis, from late 2017 until late 2021. After 2021, the rates of displacement accumulation decrease to near zero (Figure 7(a)). Location 6.75 km experiences low rates of displacement accumulation throughout the period investigated, with two long periods of stasis (Figure 7(a)). The overall distribution of slip in fault system 1 is represented by contours of aggregated M_L (third-party catalog only) as projected onto a vertical along-strike section (Figure 6(c), lower half). The contours delineate the areas of larger cumulative displacement and illustrate the segmented nature of the reactivated fault surfaces. It is clear from the contoured plot that the seismogenic, and presumably the slipped portion(s) of the fault system, is not only horizontally segmented but also limited with respect to depth, i.e., strata bound (Figure 6(c)).

Within fault system 1, there is ongoing seismicity from late 2017 to mid-2023 in the magnitude range 1.5–2.0 (Figure 7(b)). However, three large events define an upper magnitude bound that decreases over time from early 2019 to mid-2023 (Figure 7(b)). All three events occurred within one of the largest displacement maxima along the fault system (at 12.25 km; Figure 6(c)), and two correspond with a high frequency of HF completions together with high rates of SWD injections (Figure 7(c) and next section). More than 400 events occurred on fault system 1 in the period analyzed; total seismic moment (M_o) for the period was 9 × 10²⁰ ergs.

In the period October 2017 to March 2020, there were approximately 50 HF well completions in the vicinity of fault system 1. No well completions are reported after March 2020 for this fault system (Figure 7(c)). Forty-three (60%) of the completions occurred within 5 km of the NW end of the fault system, coincident with three of the largest displacement maxima (Figure 6(b) and (c)). The steepest increase in earthquake frequency and seismic energy release are both coincident in time with the most active period of hydraulic fracturing completions (end of 2017 to early 2020; Figure 7(c)). There is also a broad correspondence between earthquakes, HF, and the cumulative volume of SWD in the vicinity of fault system 1 (Figure 7(c)). Prior to December 2017, SWD injection rates increased from 0 to ~5.1 × 10³ bbls/d, this rose to 7.5 × 10⁴ bbls/d by July 2018 and 1.2 × 10⁵ bbls/d by October 2019 (Table 1, Figure 7(c)); this period of increasing injection rates was approximately coeval with high earthquake frequencies (>0.3 events/d (Figure 7(c)). SWD injection rates decreased to 8.3 × 10⁴ bbls/d from October 2019 to June 2023 (Table 1) corresponding to a decline in earthquake frequency to ≤0.2 events/d (Figure 7(c)).

The predominant vertical motion of the ground surface in fault system 1 over the period 2015 to December 2021 is continuous subsidence, although one area toward the SE end of the fault system experienced subsidence followed by uplift (Figure 8). Subsidence is typically on the order of 3 cm, but locally exceeds 6 cm. Fifteen cross-sections were constructed to elucidate the topography of the InSAR-derived surface; five representative sections are shown in Figure 8.

Only one focal mechanism solution has been calculated for events adjacent to fault system 1. This magnitude 2.7 event occurred in June 2020 at a depth of 5.3 km. Nodal planes have the orientations (right-hand rule, strike/dip/rake) 307°/36°/−117° and 159°/59°/−72°. The negative rakes of this solution accord with the most likely normal faulting stress state at this location [54]. For consistency with the inferred stress state, we choose the steeper of the two nodal planes (159°/59°/−72°) as the likely slipped plane.

All seventeen time steps for fault system 2 are illustrated in online supplementary Appendix 2.

Initial inferred displacements over the period September 2017 to April 2019 are in two distinct 3.5 km long stretches separated by ~1.5 km (Figure 9). Maximum cumulative magnitudes are concentrated at four locations at the NW and SE extremities of each stretch: 0.8, 3.3, 5.8, and 8 km from the NW end of the fault system, respectively. These relatively small displacements are best represented in Figure 9(c), although the M_L 3.82 event (TexNet catalog) at ~7.8 km from the NW end of the fault system creates a very strong maximum in M_o parameter space (Figure 9(a) and (b)).

Displacement accumulates in the two previously developed areas of displacement, and events begin to occur between the two (Figure 9). Initial displacement occurs toward the SE at 10.2 km from the NW end of the fault system.

Continued displacement results in the formation of an 8.5 km long stretch of fault system (Figure 9(a) and (c)). This is largely separated from a smaller displacement maximum ~2 km long at the SE end of the fault system (Figure 9(c)).

Fault system 2 is characterized by irregular displacement accumulation, with long periods of stasis punctuated by short periods of rapid displacement (Figure 10(a)). Overall distribution of displacement shows a concentration in an ~8 km long compound displacement maximum from ~0.2 to 8.2 km along strike. This characteristic is well represented by the cumulative M_L data (Figure 9(a)), the aggregated M_L data (Figure 9(c)), and the contoured aggregated M_L data (Figure 9(c), lower half). The six upper magnitude bounding events all occurred within this compound displacement maximum (Figure 9(a) and (c)). As in the case of fault system 1, fault system 2 is not only simply horizontally segmented but also limited with respect to depth, i.e., strata bound (Figure 9(c)).

Fault system 2 also exhibits ongoing seismicity from late 2017 to mid-2023 but over a wider magnitude range 1.5–2.5 (Figure 10(b)) than fault system 1. Again, the largest events define an upper magnitude bound that decreases over time from early 2019 to mid-2023 (Figure 10(b)). These events all occurred within the large, coalesced, displacement maxima toward the NW end of the fault system (Figure 9(a) and (c)). Fault system 2 experienced 200 events during the period analyzed; however, total seismic moment (M_o) over the same period was 4 × 10²¹ ergs. The largest events in fault system 2 are approximately 0.5 of a magnitude unit larger over the period from late 2017 to early 2023 than those in fault system 1 (compare Figures 7(b) and 10(b)).

In the period August 2017 to February 2022, there were ~15 HF well completions near fault system 2. Eleven of the completions occurred within or adjacent to the large displacement-maximum area of the fault system (Figure 9(c), lower half). Well completions in fault system 2 ended in mid-2022 (Figure 10(c)). SWD injection rates from September 2014 to January 2018 averaged 1.3 × 10⁴ bbls/d, thereafter injection rates rose to 2.1 × 10⁴ bbls/d (until March 2020), then fell to 3.1 × 10³ bbls/d (until September 2021), and rose again to 1.6 × 10⁴ bbls/d (until November 2023) (Table 2, Figure 10(c)). Despite the sporadic nature of earthquake occurrence in fault system 2 (Figure 10(a) and (c)), there is a correspondence between earthquake frequency and SWD injection rates in its vicinity, for example, the lowest injection rate period (March 2020 to September 2021) corresponds to an earthquake frequency of only 0.04 events/d (Figure 10(c)).

Fault system 2 is also characterized by continuous subsidence from 2015 to December 2021, most values being in the 3 cm range (Figure 11), but in places exceeding 5 cm. There are also small areas of reinflation both south and north of the SE-most end of the fault system (Figure 11). The most visible feature is a shallow NW–SE trending trough 1 to 2 km wide and parallel to the regional curvilinear fabric evident throughout the AOI (Figure 2(c)). Seven cross-sections were constructed to illustrate the nature of the InSAR-measured displacement of the topographic surface, and four of these are shown in Figure 11.

We use twenty-six focal mechanism solutions from the third-party earthquake catalog in the period June 2019 to May 2023 to assist with the estimation of stress state in fault system 2. All but one event occurred at depths less than 2.7 km. Strike versus rake analysis of the candidate nodal planes for the twenty-four events with negative rakes (predominantly normal displacement) indicates that S_Hmax azimuth is 133° (Figure 12).

We calculate the rake direction for all candidate nodal planes and for a range of A_φ values (Figure 12). As expected, larger A_φ values yield wider rake ranges than smaller ones, to near 0° to −180° where A_φ = 1, indicating that there is a sufficient range of potential slip surfaces to estimate A_φ. The range of reported focal mechanism rakes is −30° to −150°, with one exception; this range is consistent with an A_φ value of 0.5–0.6 (Figure 12), with our interpretation of normal slip on the faults studied, and with prior analyses of local stress state [3, 54, 55].

Previously mapped faults restricted to shallow strata above crystalline basement in the vicinity of fault systems 1 and 2 [5], when illuminated by the most likely ambient stress states for these locations [54], have critical pore pressure and slip tendency versus cumulative fault area distributions as shown in Figure 13. In the case of fault system 1, if SWD-generated pore pressure change (ΔP) is on the order of 50 psi (0.35 MPa), none of the fault area is likely to have been reactivated. However, in the case where ΔP is on the order of 200 psi (1.38 MPa), approximately 1.3 km² (7%) of the total 20 km² fault area has potential for reactivation (fault system 1; Figure 13(a)).

If SWD-generated ΔP in the vicinity of fault system 2 is on the order of 50 psi (0.35 MPa), none of the fault area is likely to have been reactivated. In the case where ΔP is on the order of 200 psi (1.38 MPa), approximately 3.7 km² (9%) of the total 38 km² fault area has potential for reactivation (fault system 2; Figure 13(a)).

Slip tendency values for previously mapped shallow faults in fault systems 1 and 2 when illuminated by the most likely preindustrial stress state [53] indicate that both systems are stable but are capable of being reactivated (T_s values ≥ 0.47; Figure 13(b)).

Although the two fault systems are broadly similar in that they consist of multiple displacement maxima, there are substantial differences in detail and in their histories of displacement accumulation over their seismogenic episodes. Fault system 1 experienced many more seismic events than fault system 2. These were distributed smoothly over time, with initially high rates of seismicity gradually decreasing over time to lower rates, and periods of stasis were rare and short-lived (Figures 6(c) and 7(a); Table 3). In contrast, seismicity on fault system 2 was irregular, dominated by long periods of stasis punctuated by short bursts of high earthquake frequency (Figures 9(c) and 10(a); Table 3). Three of the largest events on fault system 2 exceed M_L 3, whereas none of the events on fault system 1 exceed M_L 2.81 (Figures 7(b) and 10(b)). Thus, although fault system 1 hosted 418 events to fault system 2’s 200, the total seismic moment experienced by fault system 1 was much less than fault system 2 (Table 3).

The HF well completions recorded in fault system 1 form a near-continuous distribution (Figure 7(c)). Over the same period, the earthquake frequency was 0.31 events/d. Earthquake frequency then dropped gradually to 0.24 events/d by March 2021, and 0.07 events/d by November 2022, with concomitant decreases in cumulative event magnitude (Figure 7(c)). SWD injection rates recorded for the vicinity of fault system 1 varied smoothly over the period analyzed and show trends that correspond with the earthquake frequencies experienced in the fault system (Table 1, Figure 7(c)). By contrast, fault system 2 experienced only fifteen HF completions distributed sporadically over the period October 2017 to July 2022 (Figure 10(c); Table 3). In fault system 2, earthquake frequencies were on the order of 0.05 events/d, punctuated by five short-duration episodes with >0.3 and as high as 0.5 events/d (Figure 10(c)). HF well completions and earthquake frequencies in fault system 2 were both sporadic, and SWD injection rates were lower than in fault system 1 and varied more widely (Tables 2 and 3). There is, however, some correspondence between SWD injection volume rates and earthquake rates (Figure 10(c)).

InSAR measurements of the vertical displacement of the topographic surface reveal a prominent fabric manifest as troughs and ridges with amplitudes of 3 to 6 cm and wavelengths of 1 to 2 km oriented NW–SE.

Morris et al. [54] estimated the most likely stress state for the AOI defined in this work as having an A_φ = 0.794 and an S_Hmax azimuth of 130°. There are currently no reliable focal mechanism solutions in the vicinity of fault system 1; however, the general trend of the system is consistent with this S_Hmax azimuth. Strike versus rake analysis of the focal mechanism solutions available for fault system 2 is also consistent with these stress state parameters (S_Hmax azimuth = 133°), although a slightly lower value of A_φ (~0.6) is indicated.

Slip tendency and critical pore pressure measures on previously interpreted faults in the vicinity of the two fault systems indicate that there is some potential for fault reactivation in these areas. However, the SWD-generated pore pressure increases alone, as modeled by Ge et al. [7], are not sufficient to explain the observed patterns of seismicity. In addition, the analysis summarized in Figure 13, using fault interpretations from Horne et al. [5], indicates that fault system 2 was inherently less stable than fault system 1 prior to the initiation of HF and SWD activities.

The distribution of displacement inferred from the aggregated earthquake magnitude data on the two fault systems considered here strongly resembles the patterns of displacement on normal faults that have developed by initiation, growth, and linking of small faults to form larger faults both layer parallel [13, 56-58] and layer perpendicular [59-61]. Here, we focus on the horizontal manifestation of this process. Faults forming by segment growth and linkage retain a memory of the earlier segmentation in the form of displacement maxima and minima. Not all smaller faults will have linked to create larger ones and will exist as isolated features consisting of a displacement maximum between two zero displacement nodes, i.e., fault tips. We interpret the patterns of displacement to be manifestations of fault segmentation both hard-linked and soft-linked through relay structures. The histories of displacement accumulation on each of the fault systems share characteristics with the histories of HF well completions and SWD injection conducted in and adjacent to the fault systems. Fault system 1 is characterized by initial rapid displacement accumulation which transitions smoothly into decelerating displacements. The initial rapid displacement accumulation is coeval with 72 HF well completions approximately evenly distributed in time and with accelerating SWD injection rates over the same time span (Figure 7(a) and (c)). Displacement history and HF completions on fault system 2 share the characteristic of irregularity (Figure 10(a) and (c)), manifest by its wide range of HF completion rates and short duration but intense seismic episodes (Figure 10(c); Table 3). SWD activity near fault system 2 is approximately an order of magnitude lower and proportionally more variable than fault system 1 (Figure 10(c); Tables 1-3). Both fault systems exhibit a decrease in the magnitude of the largest event for a given period, over time (Figures 7(b) and 10(b)). These large events occur in fault segments of greatest displacement accumulation (Figures 6 and 9) and dominate the overall seismic energy release of each fault system (Figures 6(a), (b), 9(a) and b; Table 3).

We suggest that both fault systems had been primed for reactivation by the regional increase in pore pressure generated by SWD into the DMG and because they are both moderately well oriented for reactivation in the ambient stress state (e.g., Figure 13). However, the details of the resultant seismicity depend on the architecture of the fault system and on more local pore pressure perturbations generated by HF and SWD activities close to the faults. Preexisting fault segment size and connectivity determine the pattern and location of aggregate seismicity experienced by the fault system. The displacement maxima and minima, inferred from M_L and M_o distributions over time, are faithful representations of the distribution, linkage, and size of reactivated fault segments in the subsurface. A priori knowledge of fault system architecture would provide a powerful insight into the likely distribution of induced seismicity. The character of the pore pressure perturbation is also significant. Maximum earthquake magnitudes together with total seismic moment are lower for fault system 1 (despite higher earthquake counts, more HF completions, and higher SWD injection rates) than for fault system 2, which experienced lower intensity but highly irregular and episodic HF and SWD perturbations. This may be interpreted to mean that continuous and sustained pore pressure perturbations maintain fault slip by small increments, whereas episodic pore pressure perturbations allow stress concentrations between events leading to larger magnitudes for some individual earthquakes. An alternative explanation is that the larger events that occurred on fault system 2 were located along a major, well-connected, 8.5 km long segment of the system, while the largest connected segment on fault system 1 is only 3.5 km long. We further suggest that the background of lower magnitude events (Figures 7(b) and 10(b)) represents quasi-continuous destruction of fault system cohesion resulting in linkage of smaller slipped surfaces into a sufficiently large and connected segment capable of a substantially larger magnitude event. The gradual decrease in magnitude of these larger events over time is the result of the generally decreasing intensity of HF and SWD activities driving the triggering pore pressure perturbations.

NW–SE features of the deformed topographic surface represent a faithful map of the deformation resulting from reactivated SNF as related to subsurface oil and gas activities. Combining the data and observations from seismicity, well completions, and InSAR measurements, we develop fault reinterpretations for the two fault systems (Figure 14(a) and (b)). These are nonunique, but they support the conclusions of Horne et al. [5] that this portion of the Delaware Basin is characterized by narrow horst and graben structures oriented NW–SE, and faults that are highly segmented at a 1 to 2 km scale.

The conversion of the Delaware Basin in west Texas to a seismogenic region was the result of widespread increases in pore fluid pressure by shallow SWD and HF operations. Both activities are part of the oil and gas extraction that intensified in the region during the period 2007–2009 and has accelerated since. However, earthquakes are the result of slip on existing faults and require more than just a pore pressure increase to occur. Consequently, the finer details of earthquake location and magnitude are determined more by the characteristics of the network of faults that has been reactivated since ~2007, than simply by the presence and distribution of pore pressure perturbations. For the two fault systems analyzed here, there is a correspondence between the timing, intensity, and character of HF and SWD proximal to the fault systems and the timing and character of induced seismicity. This implies a causal link between oil- and gas-related fluid injection and the seismic response of reactivated faults. Specifically, we conclude that in a regime of broadly increasing pore pressure acting on faults that are generally well oriented for reactivation in the ambient stress state:

• Earthquakes are triggered by pore pressure perturbations causing existing fault segments to reactivate.

• Over time, the patterns of both HF and SWD proximal to a fault system are reflected in the pattern of seismic energy release on the system.

• The spatial distribution of earthquakes along fault systems is determined by the preexisting faults’ segmented structure. More and larger earthquakes will occur on larger, hard-linked segments, which are capable of hosting larger slip events without the need to overcome local impediments to slip.

• The resulting pattern of cumulative slip magnitude (inferred from the number and magnitudes of earthquakes) along a fault system will faithfully represent the segmented nature of the fault system established during its initiation and growth prior to anthropogenic reactivation.

• Knowledge of a fault system’s architecture (displacement maxima/minima, degree of segmentation and linkage, orientation within ambient stress state) can provide a predictive insight into the distribution of seismic energy release if the fault system is reactivated.

• Both slip tendency and critical pore pressure analyses of previously interpreted faults in the two locations analyzed here are fully consistent with the finding that total seismic moment for fault system 2 is greater than that for fault system 1.

Files containing well completion data, seismic event details, and InSAR data used in this work are available from the Texas Data Repository.

The authors declare that they have no conflicts of interest regarding the publication of this article.

This work was conducted as part of the authors’ employment by the University of Texas Bureau of Economic Geology.

This work was supported by the sponsors of the Center for Injection and Seismicity Research (CISR) at The University of Texas Bureau of Economic Geology. The 3DStress^(R) [62] software was used for stress analyses. We thank J.-P. Nicot and J. Ge for access to the results of the Delaware Mountain Group pore pressure model. Robert Reedy and K. Yut compiled the well database. Comprehensive and constructive reviews that greatly improved this manuscript were provided by Tim Davis and Scott Wilkins. We also thank David Ferrill, guest editor of this Lithosphere Special Edition.

Appendices 1 and 2 contain graphs of each of the seventeen, 4-month long time steps for fault systems 1 (Appendix 1) and 2 (Appendix 2). Each graph is an along-strike (NW to SE) section along the fault system. The solid line on each graph is the aggregated M_L (as a proxy for displacement) for each along-strike 500 m interval. Circular dots on the graphs are the locations along-strike and by elevation of earthquakes occurring within the basin-fill and within the area designated as contributing to fault system seismicity for each four mo time step. Events are color coded by M_L. Dates posted on the graphs are the end of each 4-month time step in decimal years. Time runs from top to bottom and left to right.