On the Depth of Earthquakes in the Delaware Basin: A Case Study along the Reeves–Pecos County Line

The Delaware basin of west Texas and southwestern New Mexico comprises the western portion of the Permian basin, and contains some of the deepest sedimentary deposits in the world. It is also one of the largest oil and gas fields in the world and, as of 2020, was producing over 4 million barrels per day of oil (U.S. Energy Information Administration Drilling Productivity Report 2/2021). The transformation of the Delaware basin from a field in decline to a giant is due entirely to production from shale source rocks, developed using horizontal drilling and “multiphase” hydraulic fracturing technology. As development proceeded, the southern portion of the Delaware basin in Texas experienced increasing rates of seismicity beginning in about 2009 that continue to the present (Frohlich et al., 2020; Skoumal and Trugman, 2021). Seismicity is not homogeneously distributed within the basin, with most events in the southern Texas portion of the basin (Savvaidis et al., 2019; Fig. 1). Although it is commonly accepted that the earthquakes are induced by oilfield activity, the mechanisms responsible for Delaware basin seismicity are under debate. Both hydraulic fracturing during well stimulation and wastewater injection have been proposed (Lomax and Savvaidis, 2019; Deng et al., 2020; Skoumal et al., 2020; Dvory and Zoback, 2021; Zhai et al., 2021).

Seismic monitoring within the basin significantly improved starting in 2017 with the installation of the TexNet array (Savvaidis et al., 2019). However, due to unique station siting challenges, the large area to be monitored, and complex wave propagation within the basin, the reported hypocenter depths vary over a wide range, from very near the surface to deeper than 12 km, with uncertainties as large as 4–5 km (Lomax and Savvaidis, 2019). Such ambiguity poses a challenge for determining the mechanisms responsible for inducing the seismicity and consequently for identifying effective mitigation measures.

This study aims to resolve the large depth uncertainty, with a focus on a subarea in the Delaware basin where the U.S. Highway 285 S crosses the Reeves–Pecos County line (Fig. 1). This small area, approximately 15 km on a side, is among the most seismically active in Delaware basin. In the four years spanning 2017–2020, the Advanced National Seismic System Comprehensive Earthquake Catalog reports 21 $ML$ ≥3 earthquakes within 5 km of the intersection and 30 within 10 km. The study area was selected owing to the presence of two stations—PB04 and PB16—that sit on top of the earthquakes. They enable a detailed examination of hypocentral depths. In addition, sonic logs and formation depths from three deep wells in the study area provide detailed control on the velocity structure from the surface to the crystalline basement (Fig. 2). All the depths discussed in this manuscript are measured relative to the ground surface.

Continuous seismic waveforms for the ten closest publicly available TexNet stations (Fig. 1b, Table S1) in operation from 1 August 2019 to 23 July 2020 are analyzed. The study period begins on the start date for station PB16, located near the center of the study area. Installation of this station gives us two near‐source stations, the other being PB04. The waveforms were analyzed using EQTransformer (EQT, Mousavi et al., 2020), an attentive, deep‐learning model designed to both detect earthquakes and time P and S arrivals. EQT analyzes seismograms in a manner analogous to a human analyst, looking for amplitude and frequency characteristics of an earthquake in addition to picking P and S arrivals. We set relatively low probability thresholds of 0.3 for event detection and 0.1 for both P and S phase picking, because EQT is robust to false positives (Mousavi et al., 2020). EQT detected over 3000 earthquakes, with the requirement that each event was detected by at least three stations. The majority (∼85%) of these detections correspond to earthquakes outside of the study area and were not considered further.

Within the study area the Delaware basin contains over 6 km of sediments. Stratigraphic information was obtained from DrillingInfo (now Enverus), along with technical details on all wells in the area. The main stratigraphic units of interest are the Upper Permian Ochoan Series evaporites (salt), located between approximately 1.0 and 1.5 km depth below the local land surface; the Delaware Mountain Group, the primary wastewater disposal formation, between 1.5 and 2.5 km depth; the Bone Springs Formation between 2.5 and 3.25 km, which contains a few producing horizontal wells; and the Wolfcamp Formation between 3.25 and 4.6 km depth where over 90% of the production wells are located, all bottoming above 3.6 km below surface.

Between 2014 and the end of our observation period in mid‐2020, 77 horizontal wells were drilled, fracked, and put into operation (Fig. 1c). These wells produced over 1.5 million barrels of oil and five times that amount of water, almost all of which was disposed of by injection into the Delaware Mountain Group. Of the 77 production wells, 8 were completed during the period for which we have seismic data. These wells locate in the northeastern quadrant of the study area and are identified in Figure 1c. Twelve disposal wells are in the study area. They all inject into the Delaware Mountain Group in the depth interval from 1585 to 2440 m. We also obtained sonic logs for three vertical wells in the study area (Fig. 1c). The logged intervals span the depth range from 10 to 6700 m below land surface.

A local 1D layered velocity model was developed using sonic logs from three deep wells in the study area (Fig. 1c). Median velocities computed every 500 ft (152.4 m) for each well show similar depth profiles (Fig. 2a). The general trend of the velocities features a rapid increase with depth in the upper 1000 m before jumping to approximately 6 km/s between 1000 and 1500 m. This high velocity layer is the Ochoan salt. Beneath the salt, velocities precipitously drop to around 4 km/s and continue to decline over the next 2 km before rising in the deep part of the basin. The lowest velocities correspond to the Wolfcamp Formation—the principal unconventional production formation in the southern portion of the Delaware basin. Basement is encountered at about 6 km depth in the study area. For purposes of event location, average compressional velocities were extracted in intervals corresponding to the formation depths, simplifying the model to 13 homogeneous layers overlying the basement (Fig. 2a and Table S2). We were unable to find dipole sonic logs near the study area. As a consequence, the shear‐wave velocity is unknown. We assumed a constant $VP/VS$ ratio, which was estimated as 1.89 using EQT picks. The detailed analysis is given in the supplemental material available to this article (Text S1).

The presence of the salt layer exerts a major impact on wave propagation from sources located below it to the surface. Ray paths for sources located below the salt would bend toward the surface if the velocities increased monotonically with depth. However, the shallow high‐velocity salt layer bends rays incident from below toward the horizontal upon entering the salt. Depending on the incidence angle, the ray refracts horizontally within the layer. Figure 2b,c compare record sections for the $Mw$ 2.8 earthquake at 17:40 on 5 August 2020, with synthetic seismograms computed using the 1D model in Table S2. The circle on each trace marks the theoretical arrival time of the P wave. In both the data and the synthetics, an impulsive P‐wave arrival is observed at the four closest stations to the epicenter, all within 10 km distance. At four slightly more distant stations, between 17 and 26 km, the first‐arriving P‐wave emerges so weakly from the noise to render its timing unreliable. Because we will be using first‐arriving P waves to locate the earthquakes, we cannot rely on automatic picks made at epicentral distance beyond about 10 km and consequently restrict our hypocenter analysis to the four closest stations—PB04, PB14, PB16, and PB18 (Fig. 1b). Stations at greater distance are used to model low‐frequency seismograms for moment tensor analysis

We approach the problem of absolute depth determination using three complementary methods: constraints from S–P times; iterative hypocenter solutions computed using P‐ and S‐arrival times; and moment tensor modeling of low‐frequency body and surface waves.

The observation of the S–P time interval at a single station constrains the maximum focal depth for the event in a vertically stratified earth model when the velocity structure is known. The earthquake can be no deeper than the depth corresponding to the S–P time from a hypothetical source located directly below the station. We use this fact together with the 300 S–P times measured by EQT at stations PB04 and PB16 and the well‐controlled velocity model to explore the range of maximum focal depths. In Figure 3a, the cumulative distribution of S–P times, sorted from the shortest to the longest, are compared with the corresponding depth for a source vertically below the station. We find that 72% of the events can be no deeper than the base of the Delaware Mountain Group at 2.5 km below ground level. In addition, 98% of the events can be no deeper than the top of the Wolfcamp at 3.25 km. If the shear‐wave velocity in any of the layers were higher (lower $VP/VS$ ratio) than assumed, the maximum depths would deepen slightly. The results for a constant $VP/VS$ ratio of 1.73 are also given in Figure 3a for comparison.

It should be noted that only 1% of the measured S–P times require a hypocenter in the Ochoan salt, directly above the Delaware Mountain group, although more of such shallow events are possible, because the S–P time only restricts the maximum depth. To better constrain the depths, it is necessary to determine hypocentral locations for individual events.

The S–P time analysis shows that many of the earthquakes locate at shallow depth, within about 1 km of the base of the Ochoan salt at 1.5 km. As discussed earlier, horizontal refraction in the salt complicates accurate timing of the initial P wave for such events, and we use only the four closest stations to locate the hypocenters. We further restrict the events to those with at least five phases. This guarantees that at least one S phase is included in the solution.

We select program VELEST (Kissling et al., 1994) to determine the hypocenters. VELEST traces rays in 1D models with embedded low‐ and high-velocity zones. Table S3 lists the 101 best-determined earthquakes, with the selection criteria including short epicentral distance (minimum station distance ≤ 2 $×$ depth), sufficient observations (≥5 phases), and small root mean square error (≤ 0.05). Standard errors for the east, north, and depth components of the solutions averaged 160, 110, and 180 m, respectively (Table S3, the fifth, fourth, and third columns from the right). The depth distribution of the events was estimated from the locations using a nonparametric kernel density estimator. Estimates made with the bandwidth set to the mean depth error or using the “rule‐of‐thumb” bandwidth estimator (Silverman, 1986) produce comparable results, the former being smoother. Approximately 75% of the density lies within the Delaware Mountain Group, between 1.5 and 2.5 km depth, with the remainder approximately equally distributed between the shallower Ochoan salt and the deeper Bone Springs formation (Fig. 3b).

Independent estimates of focal depth are obtained by modeling the waveforms of both long‐period body and surface waves. We select nine earthquakes larger than $ML$ 2.5 occurring between August 2019 and August 2020, and determined their focal mechanisms and hypocentroids using the generalized cut and paste (gCAP) method (Zhu and Ben‐Zion, 2013). The gCAP method cuts the seismograms into $Pnl$ and surface‐wave segments for inversion, and allows flexible time‐shifts to account for imperfect Green’s functions and event locations. Because we are modeling small events, we filter the waveforms into relatively high‐frequency bands to reduce the low frequency noise. We choose a 0.3–0.6 Hz frequency band for $Pnl$ wave and 0.2–0.4 Hz for surface wave.

We solve for the earthquake hypocentroid location, in addition to the moment tensor. We set the initial epicenter independent of the location determined from arrival times by first estimating the earthquake back azimuth at nearby stations (PB04, PB14, PB16, and PB18) through analysis of the P‐wave particle motion. The back azimuths are used to define a search region divided into regular grid with a spacing of 0.2 km. The grid spacing is smaller than the shortest wavelength considered. We take each grid center as a potential epicenter to perform moment tensor inversion. At each grid point, we vary the focal depth from 1 to 3.2 km at a 0.2 km interval. The preferred hypocenter has the smallest misfit between the observed seismograms and the synthetics. We limit stations considered in the inversion to those within 30 km from the epicenter for surface‐wave modeling and those within 12 km for the $Pnl$ phase. The selected $Pnl$ and surface‐wave phases are weighted equally.

Table 1 summarizes the results for the nine earthquakes. Details on the solutions appear in Figures S2–S10, including the comparison between observed and synthetic waveforms for the preferred solution and moment tensor, and misfit as a function of depth at the preferred location appears in Figure S11. Figure S12 illustrates that the centroid depths are stable with respect to varying epicenters. All of the events have normal‐faulting mechanisms striking northwest–southeast, with one high‐angle dip between 63° and 82° (Table 1 and Fig. 4). Centroid depths range between 1.4 and 2.6 km, agreeing well with the hypocenter depths discussed previously. All of the centroids are consistent with faulting in the Delaware Mountain Group within their uncertainties, although it should be noted that the two shallowest events (events 4 and 5 in Table 1) formally locate in the Ochoan salt (Fig. 3b).

On average, the centroid depths are very similar to the hypocentral depths. The 5 August 2020 $Mw$ 2.8 earthquake shown in the record section of Figure 2b has a hypocentral depth of 1.85 km, agreeing well with the moment tensor centroid depth of 2 km. The S–P time at the nearest station PB16 is 0.54 s, which rules out a hypocenter as deep as the Wolfcamp formation. If the event were vertically below station PB16, it would locate near the base of the Bone Springs at 3.2 km depth. A depth of 3.2 km at that epicenter predicts an S–P time 0.27 s longer than what is observed at station PB04.

Three independent analyses of earthquake depths in our study arrive at the same conclusion: the earthquakes occur at shallow depth, principally within the Delaware Mountain Group, where disposal of co‐produced formation water occurs by injection. We found no events in the basement or deeper than the level of production wells, and none of the best‐determined hypocenters locates in the Wolfcamp Formation. Epicenters of the earthquakes cluster along the prominent northwest‐trending subsidence lineation determined from Interferometric Synthetic Aperture Radar (Fig. 4; Pepin et al., 2021), as do the centroids determined in the moment tensor analysis. Focal mechanism solutions show normal‐faulting events on northwest–southeast‐striking fault planes (Fig. 4), consistent with the newly mapped shallow faults (Hennings et al., 2021). Their orientation and sense of slip agree well with the maximum horizontal stress orientation in this area (⁠$SHmax$⁠; Lund Snee and Zoback, 2018), which aligns with the fault‐plane strike directions.

In a companion study, Pepin et al. (2021) attributed surface deformation in the study area to slip on normal faults in the Delaware Mountain Group. Their final model of vertical and east–west deformation included three parallel, high‐angle normal faults. Two form a graben structure aligned with the prominent depression and seismicity in Figure 4, and the third is associated with the centroid locations of events numbers 4, 5, and 9 (Table 1) to the southwest of the proposed graben. Slip in their model is confined to the Delaware Mountain Group, between 1.5 and 2.5 km below ground level. Dip‐slip magnitudes reach up to 27.5 cm, greatly exceeding the slip we can attribute to the earthquakes. Indeed, only about 3% of the moment of the deformation model can be accounted for by the earthquakes. The implication is clear: most of the fault slip occurs aseismically. Pepin et al. (2021) also noted the coincidence of saltwater disposal wells actively injecting wastewater into the Delaware Mountain Group near the activated faults.

Dvory and Zoback (2021) investigated the pore pressure and geomechanical constraints on induced seismicity in the Delaware Basin. They found the Delaware Mountain Group is in a state of normal‐faulting equilibrium for high‐angle planes aligned with the $SHmax$ direction. Pressure changes resulting from disposal of wastewater by injection into the Delaware Mountain Group are capable of inducing seismic activity, and are likely responsible for the observed deformation and seismicity in the study area. Charzynski et al. (2019) analyzed seismic reflection data in a nearby location in the Delaware basin. They identified high‐angle, rootless faults as shallow graben features in the Delaware Mountain Group and hypothesized that these features could be responsible for both seismic and aseismic slip.

The coincidence of the seismicity, aseismic deformation (Dvory et al., 2021; Pepin et al., 2021), and wastewater injection to a common, restricted depth interval, where geomechanical analysis predicts northwest‐striking, high‐angle normal faults are critically stressed, strongly suggests that preexisting faults are being activated by reducing the effective normal stress via pore pressure elevation. Although the study area contains multiple horizontal production wells (∼3.6 km below the ground surface), all located below the depth of seismicity, we found no evidence of changes in seismicity coincident with the hydraulic fracturing of the eight wells (shown in Fig. 1c) completed during the study period. As noted earlier, they are all located in the northeastern corner of the study area, where we found no earthquakes. Consequently, we feel confident that hydraulic fracturing played no role in inducing these earthquakes.

We utilize three independent analyses and confirm that the earthquakes along the Reeves–Pecos county line are shallow. The shallow seismicity is triggered by wastewater injection into the Delaware Mountain Group. Other studies of seismicity in these counties, including Savvaidis et al. (2020), Skoumal et al. (2020) and Skoumal and Trugman (2021), also concluded from a variety of lines of evidence that most of the seismicity is driven by injection, without localizing it to the shallow injection zone. Seismicity temporally associated with hydraulic fracturing operations also occurs (Lomax and Savvaidis, 2019; Skoumal et al., 2020; Trugman and Savvaidis, 2021). The question remains open, however, if seismicity elsewhere in Reeves and Pecos Counties is as limited in depth extent as in the study area or extends to greater depths and, particularly, into the basement. The challenge for ongoing research will be to determine accurate focal depths. As we learned in this study, this is possible despite challenges presented by wave propagation conditions, when stations are located at close epicentral distances and when the crustal model includes local well control.

The seismic data are from the Texas Seismological Network (DOI: https://doi.org/10.7914/SN/TX) and is publicly available on Incorporated Research Institutions for Seismology (IRIS). The well data are obtained from DrillingInfo (now Enverus, https://www.enverus.com, last accessed July 2020). The supplemental material for this article includes descriptions on estimating the $VP/VS$ ratio, tables for the seismic stations, the 1D velocity profile, the hypocentral solutions, and TexNet earthquake catalog, and figures for the $VP/VS$ ratio analysis and the moment tensor solutions.

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

The authors thank Ruijia Wang, Rob Porritt, and two anonymous reviewers for their constructive comments. This work was supported by the Stanford Center for Induced and Triggered Seismicity (SCITS) and Yixiao Sheng was also supported by the Department of Energy (Basic Energy Sciences; Award DE‐SC0020445).