Insights into Temporal Evolution of Induced Earthquakes in the Southern Delaware Basin Using Calibrated Relocations from the TXAR Catalog (2009–2016) Open Access

Over the past decade, the Permian basin, west Texas, and southeastern New Mexico, United States, has become the world’s leading petroleum producer using hydraulic fracturing. The basin is divided into the western Delaware Basin (DB), central Basin Platform, and eastern Midland Basin. The increasing seismicity rate correlates with the rise in cumulative injection volumes for the whole DB. Seismicity rates, the cumulative number of M 2+ earthquakes summed by month, first increased in the southern DB, herein defined as the study region south of the Grisham fault zone inclusive of Reeves and Pecos counties (Fig. 1a) (e.g., Savvaidis et al., 2020; Skoumal et al., 2020). Figure 1b shows that the seismicity rate steadily increased alongside rising injection volumes before 2017. After 2017, a fivefold surge in saltwater disposal coincided with at least a sixfold increase in seismicity rate, though it decreased to twice the pre‐2017 rate after 2020 (Fig. 1b). Over 95% of shallow injection in the entire DB occurs within the study area. All wells inject into shallow strata of the Delaware Mountain Group and Bone Spring Formation, typically at the mean depth of 1 km below mean sea level (MSL) with no injection deeper than 2 km (Fig. 1c,d). Hydraulic fracturing commenced in the study region after 2015 at depths >2 km below MSL (Fig. 1c).

Establishing causality between oil and gas operations and seismicity beyond temporal correlations is challenging due to inaccurate location data, incomplete catalogs obscuring patterns, and limited temporal data on operational parameters. Here, the basin’s intricate geology, combined with variations in seismic network coverage over time, introduces considerable uncertainties in the hypocentral depth. The spatial or temporal gaps in waveform and reflection data complicate the identification of seismogenic faults subject to anthropogenic stress changes. Analyzing the regional EarthScope Transportable Array, Frohlich et al. (2020) resolved 10,753 epicenters in the southern DB that show seismicity began in 2009. They used P waves traveling across the Lajitas TXAR array to determine the station‐to‐source direction and TX31 S−P differential times to estimate distance (hereafter TXAR catalog; 2009–2016). The majority of TXAR events have epicentral uncertainties of 10–15 km and depth estimates from stations located 75–450 km away, with uncertainties exceeding 5 km. Reducing spatial uncertainties in the TXAR catalog and tracing temporal changes in induced earthquakes back to 2009, we seek to improve our comprehension of the original causal factors behind seismicity in the southern DB.

We employ the hypocentroidal decomposition (HD) algorithm (Bergman et al., 2023) specifically designed for calibrated relocation of regional catalogs using local recordings. Our local dataset comprises a core cluster of 116 post‐2020 events sourced from the TexNet catalog (Savvaidis et al., 2019). Through multistep clustering and data calibration, we establish station–phase–event connections between pre‐ and post‐2017 seismicity using common regional stations. Like other multievent relocation procedures, the HD method reduces bias caused by unaccounted subsurface structures along lengthy event–station paths using differential times that represent shorter paths within the cluster (Fig. 2). This method also utilizes travel‐time residual statistics for outlier removal and data weighting. We relocate 73 events from the TXAR catalog using catalog phase times and adding cross‐correlation differential times of paired events at regional stations. The relocation results help examine spatial and temporal connections between pre‐ and post‐2017 seismicity, mapped faults, and oil and gas operations in the southern DB.

Travel‐time methods are particularly effective in pinpointing shallow and small earthquakes, especially in scenarios with sparse data, suboptimal station geometry, and uncertain subsurface structures. The widespread adoption of waveform cross‐correlation for calculating differential time has bolstered source characterization accuracy across diverse scales and tectonic settings (e.g., Aziz Zanjani et al., 2021; Aziz Zanjani and Lin, 2022; Trugman et al., 2022). However, akin to certain travel‐time methods, the cross‐correlation requires pairing earthquakes at seismic stations and lacks simultaneous enhancement through bridging local and regional catalogs across various time frames.

The HD relocation hinges on two separate stepwise inversions (Fig. 2) (Bergman et al., 2023). It determines calibrated hypocenters with respect to a virtual point in space and time called the hypocentroid (HYP), here the centroid of a core cluster (Fig. 2b). The relocation analysis utilizes a locally appropriate velocity model and core‐cluster data to independently estimate the HYP using true travel‐time residuals at near‐source distances and avoids observations near the P_g/P_n crossover to minimize absolute location errors. Subsequently, the relative hypocenter for each event in the cluster is calculated with respect to the HYP (cluster vectors) (Fig. 2d), representing differential times instead of absolute times to mitigate bias introduced by the assumed velocity model (Bergman et al., 2023). After completing the calibration and relocation of the core cluster with the inverse weighting of phase data using station–phase residual statistics, we integrate the TXAR events. Regional station statistics assist in identifying and removing outliers within the TXAR catalog. We then link the locations of post‐ and pre‐2017 seismic events using shared regional stations (Fig. 2).

The HD algorithm provides quantitative measures of absolute and relative hypocenter errors. Statistical analysis iteratively employs empirical reading errors (EREs) to estimate observation uncertainties, remove outliers from arrival times, refine data weighting, and improve both the velocity model and focal depths. EREs are a measure of residual spread (⁠$σ$⁠) for each station–phase pair in the cluster. Larger EREs indicate increased scatter due to arrival‐time picking discrepancies, station variations, and timing irregularities. Data calibration ensures a normal distribution for each station–phase pair, with phase arrivals within $2σ$ of the empirical path anomalies, representing the mean of residuals (Bergman et al., 2023).

The selection of the core cluster based on data coverage in time and space ensures minimal location errors for the HYP. Its spatial dimension allows for smoothing out variations in the local structure around the HYP along multidirectional station–event paths (Fig. 3a). We limited the number of events in the inversion to a maximum of 200 to streamline the calibration process. The core cluster includes 116 earthquakes from the post‐2020 TexNet catalog with 2289 P_g and 2006 S_g within 1° radius, augmented with 290 S−P differential times within 10 km radius to constrain depths. Figure 3a illustrates the distribution of station–event paths for the core cluster within the distance limit of 0.6° for calculating HYP (Fig. 3b). Event selection for the core cluster is guided by a comprehensive analysis of post‐2017 seismicity, employing the HD algorithm (Aziz Zanjani et al., 2023). We ensure spatial coverage around the HYP through local stations and uniform event distribution, prioritizing dense local paths to minimize uncertainties (Fig. 3a). All events in the core cluster are M 1.8+ with minimum azimuthal coverage of 200°. The core cluster includes 4295 phase arrivals within a 1° radius with 81% identified as P_g and S_g arrival times. During data calibration using EREs, 14% of readings within this distance were removed from the inversion. Figure 3b presents the travel time for P_g and S_g (as direct local seismic phases), as well as P_n and S_n (as regional seismic phases), and includes S−P differential times for the entire cluster after calibration. The theoretical travel times calculated from the velocity model at HYP depth (1.8 km below MSL) (Fig. 3c), accurately fit the observed moveout of P_g and S_g arrival times (Fig. 3b,d). At the 0.6° cutoff distance, travel‐time residuals are within ±1 s, with mean residuals of −0.09 and 0.12 s for the P‐wave and S‐wave, respectively, and have a smaller spread around mean for the P‐wave compared to S‐wave (Fig. 3e,f). At the cutoff distance, S−P differential times have the smallest spread and a zero mean (Fig. 3g). The core cluster provides minimally biased absolute hypocenters, meaningful errors, and curated arrival times as a reference for calibrating the HYP and the TXAR regional catalog.

We relocate M 1.5+ earthquakes from the TXAR catalog with sufficient azimuthal coverage for a high‐quality location to explore the major loci of energy released before 2017. Only 5.4% of the TXAR catalog within our study area are M 2+; 62 of the 73 relocated TXAR events are M 2+. The selected TXAR events have at least three stations within 2° and are located within the 3D spatial volume of the core cluster. The calibrated pre‐2017 events include 1081 P_g and S_g arrival times with 20% within a 1° radius and an additional 798 P_n and 393 S_n readings. During data calibration using EREs, 23% of readings within a 1° distance were removed from the inversion. The general workflow includes the stepwise addition of events to the core cluster, flagging outliers, and defining focal depth. Events with multiple outliers at regional distances lose azimuthal coverage after calibration, contributing to their large uncertainty.

Figure 4a illustrates the relocation results with relative errors for the pre‐and post‐2017 events. Our calibrated dataset includes 3039 P_g, 2954 S_g, 1974 P_n, 529 S_n arrival times, 384 single‐station S−P differential times, and 484 differential times from waveform cross‐correlation of paired events at TX31 (co‐located with TXAR). The calibrated HYP is located at coordinates 31.155° N, −103.349° W, and 1.8 km below MSL, with errors of 0.04 s in origin time and 110 m in latitude and longitude, providing a measure of absolute errors for the entire cluster (Fig. 4a). The majority of post‐2017 seismicity is shallower than 2 km below MSL—the bottom of the Bone Spring Formation (Fig. 4b–d). The depth distribution of pre‐2017 events is similar with larger depth errors and 88% shallower than 2 km below MSL (Fig. 4b–d). The deeper pre‐2017 events systematically show larger depth errors (Fig. 4b). The mean value of depth error for TexNet and TXAR events are 300 m and 1.8 km, respectively (Fig. 4d). The epicentral uncertainties for relocated post‐2017 events are negligible (<300 m in latitude and longitude), and error ellipses are less than $2km2$ with a 95% confidence interval (Fig. 4a,e). The shifts in epicentral locations relative to initial locations for the TXAR events underscore the significance of biases introduced by 1D velocity model using a single‐event algorithm (gray lines in Fig. 4a). Each event in the core‐cluster features depth constraint from seismic stations within 10 km and S−P differential times with errors within ±500 m (Fig. 4d). Pre‐2017 error areas are 3–6 times larger due to the lack of local phase data and poor station coverage (Fig. 4e). The larger east‐northeast–west‐southwest errors for relocated TXAR events reflect the direction of poor azimuthal coverage from regional stations (Fig. 4a,e).

Seismicity in the southern DB has been linked to shallow injection (e.g., Pepin et al., 2022; Sheng et al., 2022) or a combination of saltwater disposal and hydraulic fracturing (e.g., Savvaidis et al., 2020; Hennings et al., 2021; Grigoratos et al., 2022). Shallow normal faults in this region are primarily located within the Delaware Mountain Group, extending into the overlying Ochoan salts or underlying Bone Spring Formation, and they are optimally oriented in the direction of maximum horizontal stress (Hennings et al., 2023). Unlike deeper basement faults requiring larger pressure changes for reactivation, these shallow faults are in a near‐critical stress state before anthropogenic triggering, with uniformly low deterministic fault‐slip potential (Hennings et al., 2021). Therefore, they are susceptible to seismic slip under slight increases in pore pressure and the commensurate decrease in normal stress. A pore‐pressure increase of 2.5 MPa destabilizes ∼71% of the fault length (Dvory and Zoback, 2021). Using injection data, Ge et al. (2022) suggest pore‐pressure changes of up to 2.8 MPa from saltwater disposal into the Delaware Mountain Group, reaching 3.5 MPa near high‐volume saltwater disposal, beyond the slip‐triggering limit for these faults.

Our relocation does not include the TXAR catalog in its entirety, which restricts a comprehensive analysis of pre‐2017 seismicity, but focusing on the largest events unveils noticeable spatiotemporal patterns. Figure 5a–d compares the spatial correlation between the 73 relocated pre‐2017 events and active seismic lineaments (hereafter ASLs) from relocated post‐2017 seismicity using the HD method (Aziz Zanjani et al., 2023). Except for seismic cluster 1 in Figure 5a, all relocated pre‐2017 events spatially associated with post‐2017 ASLs, interpreted as northwest‐trending shallow normal faults in former studies. Figure 5 also shows the correlation between seismicity and spatially smoothed cumulative injection volumes using the inverse distance weighting deterministic method adopted by Hornbach et al. (2016). The earliest 2009–2012 seismicity, clusters 1–3 and event 4 in Figure 5a,b, occur 5–20 km from high‐volume injectors at the time (locations A and B in Fig. 5a). However, event 5 is located at least 60 km away from injections. As time progresses, the loci of seismicity shift from north to the southeast of the study area, with clusters 8–14 occurring less than 20 km from injector locations C, D, and E (Fig. 5d). We suggest that this seismicity is induced by near‐field effects of pore‐pressure changes due to increases in cumulative injection volume or new saltwater disposal at locations C, D, and E (Fig. 5d).

Other processes can change the relationship between shear and normal stress to trigger events on critically‐stressed faults. Zhai et al. (2018) demonstrate that although poroelastic stress changes have a limited effect on seismicity in Oklahoma, their combination with pore‐pressure changes can substantially increase the seismicity rate by 2–6 times. Goebel et al. (2017) underscore the critical role of poroelastic stress change in triggering earthquakes 10–40 km away from saltwater disposal wells in Oklahoma. The rate of fluid pressurization impacts effective normal‐stress unloading, leading to frictional sliding destabilization with abrupt and a large increase in injection rates when stabilizing under gradual or small injection rates (Alghannam and Juanes, 2020). High‐volume injection to the northern part of the southern DB, combined with changing rates, could theoretically increase the stress perturbation in geologic units below the shallow injectors and/or at distances of 10s of kilometers, leading to the increased number of earthquakes in the southeast section of our study area over time. We favor pore‐pressure diffusion as the primary triggering mechanism for near‐field earthquakes, with any other effects considered secondary.

Far‐field pressure diffusion is documented in the DB and other basins (e.g., Hornbach et al., 2016; Peterie et al., 2018; Savvaidis et al., 2020; Jin et al., 2023). In the Dallas–Fort Worth area of north‐central Texas, it is shown that the faults near high‐volume saltwater disposal destabilize first, followed by delayed triggering of distant faults up to 40 km away at lower pressure‐change levels (Hornbach et al., 2016). Peterie et al. (2018) associated increased seismic activity in southern Kansas in 2013 with heightened pressure diffusion from distant high‐volume saltwater disposal. Ogwari et al. (2018) documented similar effects for the Dallas–Fort Worth Airport earthquakes, where swarm‐like seismicity near a well within days of initial injection migrates along the fault and away from the well years after injection ceases. Skoumal et al. (2020) infer the far‐field effects of saltwater disposal inducing earthquakes within a distance greater than 25 km in the southern DB, south of the Grisham fault zone. The southeast shift we observe in the loci of seismicity spatially correlates with porosity and permeability changes in the subsurface strata, as documented by Smye et al. (2021), indicating heightened fluid capacity and effective pressure diffusion along highly porous and permeable sandstone within the Delaware Mountain Group. The clockwise northwest–southeast migration of the sandstone depocenter in the study area correlates with strata‐bounded fault trends and shallow seismicity, where faults serve as permeable and anisotropic conduits for pressure migration (Hennings et al., 2023). Here, we assume the Grisham fault zone effectively compartmentalizes the hydrogeologic system. Sensitivity analyses by Ge et al. (2022) show that pressure build‐ups in these basins are sensitive to horizontal permeability and compressibility, supporting directional pressure diffusion in our study region.

Based on frictional properties of cores sampled from the Delaware Mountain Group, Bolton et al. (2023) favor aseismic creep for shallow faults in the DB, but other studies have associated both seismic slip and aseismic creep to faults in the southern DB (e.g., Pepin et al., 2022; Sheng et al., 2022; Hennings et al., 2023). Modeling of Jeong et al. (2023) shows that the level of pore‐pressure perturbation positively correlates with aseismic slip, characterized by low‐stress drops, which can either advance or delay seismic slip on faults. Eyre et al. (2019) propose a seismic‐slip initiation model where an aseismic creep front interacts with fault regions exhibiting dynamic and slip‐rate‐weakening behavior. Although intricate spatiotemporal seismicity patterns in the southern DB provide valuable insights into triggering mechanisms of induced earthquakes over the initial eight years of seismic activity, the earthquakes cannot speak to the balance between aseismic and seismic processes. Analyzing larger magnitude events in the TXAR catalog, we suggest that pre‐2017 seismicity was primarily triggered by the interplay between near‐field and far‐field effects from saltwater disposal on optimally oriented faults. Although the involvement of hydraulic fracturing cannot be dismissed, it appears relevant to only a small percentage of events occurring at depths exceeding 2 km below MSL in our dataset.

Our new catalog offers improved depth resolution for pre‐2017 events compared to previous regional studies (e.g., Frohlich et al., 2020; Skoumal et al., 2020). Skoumal et al. (2020) associates the majority of pre‐2017 seismicity with saltwater disposal and link only 5% to hydraulic fracturing. Our approach is based on very distinct volumes and depth intervals of hydraulic fracturing versus shallow saltwater disposal in our study region (Fig. 1c). The initial depth distribution for both TexNet and TXAR events is presented in Figure 5e. Note the default depth of zero reported for 33 events in the TXAR catalog has shifted the mean value (red curve in Fig. 5e). The mean depth value for the TXAR events after relocation is 1 km below MSL (Fig. 5f). A comparison of the depth distribution of 62 common relocated events between this study and Savvaidis et al. (2022) also reveals mean values at 2 km below MSL (Fig. 5g). Shifts in mean depth values from initial catalogs to shallower depths after relocation (Fig. 5e,f) leads to improved correlation with the average injection depth, typically 1 km below MSL in our study area (Fig. 5h). The shallow depths (≤2 km) for most earthquakes (85%) are consistent with studies analyzing post‐2017 seismicity (Savvaidis et al., 2022; Trugman et al., 2022; Aziz Zanjani et al., 2023) (Fig. 5g). Moment tensor solutions confirm shallow depths for post‐2017 seismicity in this region and indicate a direct link with shallow saltwater disposal (Sheng et al., 2022; Aziz Zanjani, et al., 2024). The pre‐2017 events remain below moderate magnitudes (M < 3.5) like post‐2017 shallow seismicity, and their depth distribution suggests the reactivation of high‐angle normal faults within the Delaware Mountain Group and the Bone Spring Formation. Because the relocated catalog still contains a small number of events below the Bone Spring Formation with larger errors, we cannot rule out the role of hydraulic fracturing in triggering some pre‐2017 events in 2015 and 2016. The debate for deeper seismic activity within the crystalline basement also remains inconclusive because our core cluster did not include many deeper events. Further study of the TexNet catalog would best explore the potential for reactivation of deep basement‐rooted faults. Smye et al. (2021) posit direct triggering is unlikely due to the limited vertical extent of saltwater disposal into the Delaware Mountain Group and intervening layers, but that reactivation through poroelastic effects exists.

We compare 73 calibrated hypocenters from the TXAR catalog to the active faults identified from post‐2017 seismicity, shallow saltwater disposal, and hydraulic fracturing in the DB, south of the Grisham fault zone. Analyzing local–regional seismic phase data of 116 post‐2020 earthquakes in the TexNet catalog with the HD method, we minimize epicentral errors of 73 pre‐2017 M 1.5+ earthquakes in the regional TXAR catalog to less than 5 km and report new depths. After relocation, pre-2017 events spatially correlated with post-2017 seismicity reveal the reactivation of critically stressed faults within the Delaware Mountain Group and Bone Spring Formation, typically at a mean depth of 1 km below MSL. The depth distribution of relocated events strongly correlates with evolving shallow saltwater disposal to the south of the Grisham fault zone since 2009. Pre-2017 seismicity patterns suggest the main loci of seismicity migrated from the north to the southeast of the study region. Because pore-fluid pressure is diffused from the wells, operational changes became evident, with new injectors emerging to the southeast and an increased volume of shallow injectors in the north. The southeastward migration of seismicity from 2009 to 2016 supports the previous research suggesting an efficient pressure diffusion along the permeable faults and lineaments and thickening of highly porous sandstone in the Delaware Mountain Group. Local pore-pressure changes and the poroelastic stress changes from distant shallow saltwater disposal contribute to increased shallow seismicity rates to the southeast from 2009 to 2016. In summary, we show that the initial seismicity in the southern DB occurred due to near‐ and far-field stress changes associated with shallow saltwater disposal south of the Grisham fault zone on shallow normal faults in the Delaware Mountain Group and Bone Spring Formation, many of which remain active over a decade later.

The TexNet catalog is from the Texas Seismological Network (doi: 10.7914/SN/TX). Injection data were obtained from the Railroad Commission of Texas (RCC; https://www.rrc.texas.gov, last accessed March 2024). Hydraulic fracturing data were downloaded from FracFocus (https://fracfocus.org/, last accessed March 2024). The MLOC (v10.6) is available at https://seismo.com/mloc/ (last accessed March 2022).

The authors acknowledge that there are no conflicts of interest recorded.

The authors thank Brian Stump, Keith D. Koper, Tom Green, and two anonymous reviewers for their constructive comments and feedback. The authors gratefully acknowledge funding from the State of Texas through TexNet and thank the Railroad Commission of Texas (RCC) and FracFocus for providing operational data to the public. In addition, the authors thank Katie Smye for her assistance with injection data and Eric Bergman for the Multiple Relocation (MLOC) open‐source code. The authors thank Lili Binetti for providing regional TXAR cross‐correlation data.