Role of Deep Fluid Injection in Induced Seismicity in the Delaware Basin, West Texas and Southeast New Mexico

Rates of seismicity in the Delaware Basin of Texas and New Mexico increased from 10 earthquakes per year of local magnitude (ML) 3.0 and above in 2017 to more than 185 in 2022, coincident with increasing oil and gas production and wastewater re‐injection into strata shallow or deeper than producing intervals. Events of large magnitude—up to ML 5.4 to‐date—occur on faults extending into formations above the basement that have received more than four billion barrels of injection. Here, we demonstrate the link between injection geology, pore pressure evolution, fault stability, and induced seismicity in this region. We find that the injection targets are largely dolomitized platform carbonates with low (<5 vol.%) matrix porosity and fracture‐enhanced permeability with inherent heterogeneity in flow properties. A comprehensive, three‐dimensional geological model populated with reservoir properties is used for fluid flow modeling, with global calibration supplemented by dynamic injectivity data. Pore pressure changes with deep injection are up to 5 MPa from 1983 to 2023, increasing the native pore pressure state by 10% locally. Modeling results show that earthquakes occurring at distances of up to 30 km from deep injection have experienced small (<0.1 MPa) pore pressure increases, indicating that the faults hosting these earthquakes are highly sensitive to changes in effective stress and have lower frictional stability than the 0.6 generally assumed. These results serve as a critical step in understanding the stress changes that induce earthquakes in one of the most seismically active and geologically complex basins in the US.


Introduction
The Delaware Basin (DB) of West Texas (TX) and southeast New Mexico (NM) (Figure 1) is one of the most seismically active, geologically complex, and hydrocarbon-productive basins in the U.S.Although the number of conventional vertical wells drilled per year declined ∼75% from 2012 to 2022, horizontal drilling and hydraulic fracturing (HF) to recover hydrocarbons from low-permeability formations increased by over 300% in that same time period.More than 20,000 horizontal wells have been drilled in the basin, with drilling rates peaking in 2018-2019, and again in 2022, at ∼3,000 wells per year, targeting formations of Wolfcampian and Leonardian age (Figure 2).More than 5 billion barrels (Bbbl) of oil and 17 trillion cubic feet of gas have been produced from horizontal wells in the basin (S & P Global, 2023).
Unconventional hydrocarbon production in the DB results in the co-production of volumes of HF flowback and produced formation water that are more than double the volume of oil produced per well (Scanlon et al., 2020).This oilfield wastewater has high salinity (hence the term "saltwater disposal," or SWD) and other impurities including organic compounds, and it must either be recycled for reuse in HF or injected into formations with porosity and permeability sufficient to accommodate fluid storage.
Approximately 20 Bbbl of wastewater have been injected into the subsurface of the DB since 2010 for permanent disposal; 75% of that has gone into Guadalupian-age Delaware Mountain Group sandstone reservoirs (Ge et al., 2022;Lemons et al., 2019;Smye, Banerji, et al., 2021) that are shallower than producing shale intervals (Figures 2 and 3).This is termed "shallow injection" with a depth range of 450-2,400 m (1,500-8,000 ft), averaging 1,370 m (4,500 ft) in the DB.More than four billion barrels have gone into deeper carbonate-dominated formations of late Ordovician, Silurian, and Devonian age (Figures 2 and 3) (Gao et al., 2020), with injection depths ranging from 1,370 to 5,800 m (4,500-20,000 ft), averaging 3,570 m (11,700 ft).Although secondary in volume to shallow injection, deep injection is of primary importance in induced seismicity in the northern DB given that deep disposal is the most vertically proximal operational activity to seismicity (Figure 4), which occurs mainly in the crystalline basement.Vertical proximity of injection zones to the basement, or more generally, greater depth of injection, has been shown to be a risk factor for induced seismicity in many oilfield wastewater injection scenarios (e.g., Hincks et al., 2018;Skoumal & Trugman, 2021;Skoumal et al., 2018;Skoumal et al., 2020Skoumal et al., , 2021)).
Of the regions in the midcontinent of the U.S. that have little natural seismicity, the DB has seen the greatest increase in rates of seismicity in recent years, beginning in 2009 and further accelerating from 2016 (Ellsworth, 2013;Frohlich et al., 2020;Lomax & Savvaidis, 2019;Savvaidis et al., 2019).Since robust monitoring began in 2017, there have been more than 40 earthquakes of local magnitude (M L ) ≥ 4.0 in the DB recorded by the Texas Seismological Network (TexNet, 2023) and the New Mexico Tech Seismological Observatory (NMT).Earthquake rates increased from 10 per year of M L ≥ 3.0 in 2017 to more than 185 in 2022 (TexNet, NMT).In the southern DB earthquake zone (SDBEZ) (Figure 1), earthquakes frequently occur in the sedimentary section rather than in the crystalline basement (Figure 4), and shallow SWD and HF have been shown to be causal agents (Grigoratos et al., 2022;Savvaidis et al., 2020;Zhai et al., 2021).These earthquakes appear to be limited to M L ≤ 4 and the earthquake rate in this sub-region has declined in recent years (Hennings, Dvory, et al., 2021).However, in the northern DB earthquake zone (NDBEZ) of NM and TX addressed in this study and delimited on its southern end by the Grisham Fault Zone (GFZ) (Figure 1), earthquake rates have remained elevated.This region includes a dense area of seismicity spanning Culberson and Reeves Co., TX, called the Culberson-Mentone earthquake zone (CMEZ) (Figure 1), including the M L 4.9 "Mentone" earthquake that occurred in 2020 and the M L 5.4 "Coalson Draw" event that occurred in late 2022.Other notable events in the NDBEZ include the M L 4.0 "County Line" earthquake spanning Eddy and Lea Co., NM, which occurred in 2021 (NMT).
Hypocentral depths for earthquakes reported in the NDBEZ are consistent with events occurring in the crystalline basement, with average event depths of 7.7 km (25,000 ft) below the ground surface and 6.6 km (22,000 ft) below mean sea level (MSL) (Figure 4), and average depth uncertainty of 1.5 km (4,900 ft) for TexNet events and depths of 5 km reported for NMT events.Probabilistic associations between operational activities and earthquakes show a likely causal influence of deep injection on earthquakes in NM and in Culberson Co., TX (Grigoratos et al., 2022).However, many of these earthquakes occur at distances of 30 km (20 mi) from injection wells, and probabilistic associations with deep SWD are weak, leading previous authors to interpret cascading earthquake interactions initially induced by deep SWD.While deep injection has been implicated in induced seismicity in this system (e.g., Skoumal & Trugman, 2021), analysis of pore pressure changes associated with deep injection in these units using comprehensive models has not been carried out, and is a key step in causal assessment of induced seismicity regionally, and for analogous application elsewhere.
Here we seek to understand the role of deep injection in induced seismicity in the northern DB through modeling of pore pressure changes with injection using a detailed geologic model of the injection stratigraphy.The link between induced earthquakes and anthropogenic activities is complex in this region given the varied subsurface operational activities, involvement of multiple stratigraphic intervals, and the complex geologic character, necessitating mechanistic analyses to understand the role of specific operational activities in induced seismicity and to enable mitigation.Mechanistic analysis of induced seismicity begins with detailed geological characterization and assessment of stress change from pore pressure evolution modeling-we provide these contributions here and relate them to earthquake occurrence.

Geologic Background
The DB of West TX and southeast NM is part of the greater Permian Basin (PB) region, which encompasses all or part of 37 counties in Texas and NM and covers ∼90,000 km 2 (34,000 mi 2 ).Formations targeted for deep injection and studied here are largely pre-and syntectonic with respect to structural differentiation of the Delaware and Midland sub-basins and the intervening basement-rooted uplifted Central Basin Platform (CBP) (Figure 1).
The PB region extends across several basement terranes which assembled during the Proterozoic period from 1.30 to 1.07 Ga.These terranes include the southern crystalline area, the northern crystalline area, and the southern granite rhyolite province with an associated continental margin arc batholith extending from the PB to the eastern Fort Worth Basin (Adams & Keller, 1996;Denison et al., 1971;Ewing, 2019;Flawn, 1956;Lund et al., 2015;Muehlberger, 1965;Muehlberger et al., 1967).This region is also impacted by the NW-verging collisional event known as the Grenville Orogeny (∼1.2-1Ga) (Keller et al., 2005;Mosher, 1998;Mosher et al., 2008).The Grenville front (Figure 1) bisects the region and likely influences later patterns of deformation through tectonic inheritance, including basement-rooted fault systems (Adams & Keller, 1996;Davis & Mosher, 2015;Horak, 1985;Horne et al., 2021Horne et al., , 2022Horne et al., , 2024)).(Horne et al., 2021(Horne et al., , 2024) ) and regional faults (Ewing, 1991), study area (gray dashed polygon), maximum horizontal stress orientations (Dvory & Zoback, 2021), and Seismic Response Areas as determined by regulators in New Mexico and Texas.CMEZ = Culberson-Mentone Earthquake Zone; GFZ = Grisham Fault Zone; AOI = modeled area of interest, NDBEZ = Northern Delaware Basin Earthquake Zone (deep earthquakes in basement), SDBEZ = Southern Delaware Basin Earthquake Zone (some shallow earthquakes in the sedimentary section).Large-magnitude events discussed in text are noted.Above the crystalline basement is a thick (up to 2,400 m or 8,000 ft) section of pre-Permian Paleozoic rocks, with deposition beginning on a flooded continental margin with an extensive shallow-water carbonate platform that developed during Eocambrian or Middle-Late Proterozoic rifting, forming the Tobosa Basin (Derby et al., 2012;Hills, 1984;Ruppel, 2019a).The Ordovician-age Ellenburger Group (Figure 2) is likely the oldest of these rocks in the PB, apart from siliciclastics of unknown age that sit above the basement locally, filling Precambrian basement paleotopographic lows (Ruppel, 2019a).Along with its age-equivalents including the Arbuckle Group, the Ellenburger Group has been targeted for deep injection in other basins in the midcontinent of the U.S. including the Fort Worth Basin (Gao et al., 2021;Hennings et al., 2019;Hornbach et al., 2016;Smye et al., 2019), Midland Basin (Woo & Ellsworth, 2023), and in basins of Oklahoma (Walsh & Zoback, 2015;Weingarten et al., 2015).Carbonates of the Ellenburger Group are dominated by karsted limestone and dolomite of peritidal and shallow subtidal origin (Loucks & Kerans, 2019).Karst features are widespread and likely regionally extensive (e.g., Sanchez et al., 2019;Smye et al., 2019), but correlation of these exposure surfaces in the subsurface has proved challenging.

Geochemistry, Geophysics, Geosystems
Global sea level fall in the late Early Ordovician to early Middle Ordovician resulted in the Ellenburger Group being unconformably overlain by the Simpson Group (Figure 2), which contains increased siliciclastic content related to reflooding of the shelf.Source areas for Simpson Group sandstones were likely to the north and west of the PB region (Harrington, 2019).Late Ordovician and Early Silurian deposition on the shallow-water platform, including that of the Montoya Group and Fusselman Formation, was dominated by carbonates with local terrigenous clastics sourced from highlands in present-day northern NM (Harrington & Ruppel, 2019).Based on cores and cuttings (Harrington & Ruppel, 2019;Thomas & Liu, 2003), outcrops (Howe, 1959;Pope, 2004Pope, , 2014)), and subsurface studies (Ball, 2003;Behnken, 2003), Montoya Group lithology is interpreted as dolomitic in the northern part of the basin and cherty in the more distal parts of the basin to the south, with prevalent limestone elsewhere (Harrington & Ruppel, 2019).The Fusselman Formation represents the continued development of this  (Ewing, 1991), basement-rooted faults (Horne et al., 2021(Horne et al., , 2024) ) and shallow faults (Horne et al., 2022)  shallow-water platform and consists of carbonate facies that reflect the generally east-west trending platform margin (Canter et al., 1992) with shallower shelf conditions to the north.
Continued deepening and northward shift in depositional environments occurred throughout the Middle and Late Silurian and Early Devonian, including deposition of the Wristen Group and Thirtyone Formation.Distribution of lithofacies reflects not only deepening and proximal-distal positioning but also differential subsidence and the formation of a Silurian carbonate platform margin (Figure S1 in Supporting Information S1) that bisects the footprint of the DB and appears to control sediment distribution patterns through the Early Devonian (Ruppel, 2019c;Ruppel, Hovorka, & Barnaby, 2019).Perhaps importantly, this platform margin turns south in the western DB, such that Wristen Group rocks to the north and west are largely shallow water platform carbonates, and to the south and east are dominantly mudrocks.Facies distributions of the Early Devonian Thirtyone Formation are controlled by the location of the underlying platform margin, and by the transition from dolomitized platform carbonate deposits in the north and west (proximal) to deepwater chert and carbonate deposits to the south (distal).Differentiation of Silurian and Devonian rocks in the subsurface is challenging due to their lithologic similarity, and to truncation and removal by Middle Devonian uplift and erosion (Ruppel, Hovorka, & Barnaby, 2019) that is evidenced by karsting of Silurian and Early Devonian rocks.
Following marine flooding, the Late Devonian Woodford Formation was deposited on this karsted Silurian and Early Devonian surface, locally infilling topographic lows.The lithology of Woodford Formation is organic-rich (up to 25%) mudrock (Ruppel, Rowe, et al., 2019).A sea-level highstand in the Mississippian resulted in deposition of the Mississippian Limestone associated with an extensive carbonate platform, followed by deepwater mudrock deposition of the Mississippian Barnett Formation.
The Pennsylvanian was marked by major tectonic activity associated with the final collision of the Laurentian and Gondwanan plates, resulting in the formation of the structural features that define the geometry of the PB, including the CBP separating the Midland and Delaware Basins (Ewing, 1990(Ewing, , 1991;;Horne et al., 2021;Yang & Dorobek, 1995).Sedimentation patterns in the Pennsylvanian reflect this structural differentiation, with early carbonate and mixed carbonate-siliciclastic shelf deposits followed by Late Pennsylvanian basinal facies in the DB, surrounded by carbonate platforms to the west, north (Northwest Shelf) and east (CBP) (Wright, 2011).Early (Wolfcampian) and Middle (Leonardian) Permian shale development intervals are dominantly mixed carbonatesiliciclastic systems.Importantly, these units exhibit a degree of overpressure that has been maintained for millions of years (Sinclair, 2007;Luo et al., 1994), indicating that they have very low permeability and serve as an effective no-flow boundary at the top of the modeled domain.
Two distinct and seemingly kinematically disconnected sets of faults occur in the DB (Figure 3).The first set are characterized as contractional basement-rooted faults, which form the general structure of the basin and separate it from the N-S-trending uplifted CBP (Horne et al., 2021).These faults are rooted in the Proterozoic basement and tip out vertically in lower portions of the Wolfcamp Formation, therefore impacting the deep injection intervals of Silurian and Devonian age, and appear to host most of the deep earthquakes in the DB (Hennings, Dvory, et al., 2021;Horne et al., 2021).Notably, the E-W-trending GFZ (Figure 1) consists of an important set of basement-rooted faults that bisect the DB and fully cut the deep injection units.The second set is characterized as a shallow, extensional fault system that is mostly strata-bound to the Delaware Mountain Group (Horne et al., 2022).These faults trend parallel to the azimuth of maximum horizontal stress (S Hmax ), dip steeply, host linear arrays of shallow earthquakes, and frequently displace the ground surface, as observed in surface deformation maps from the Interferometric Synthetic Aperture Radar (Hennings, Dvory, et al., 2021;Horne et al., 2022;Pepin et al., 2022;Staniewicz et al., 2020).
The state of neotectonic stress in the PB region transitions from a normal faulting regime in the west (DB) to a normal/strike-slip regime in the east (CBP and Midland Basin) based on in situ measurements, geologic and seismogenic data (Lund Snee & Zoback, 2018).In the normal faulting stress field of the DB, vertical stress (S v ) is the greatest principal stress (S 1 ), and its magnitude is 0.0239-0.0240MPa/m (1.06-1.1 psi/ft) (Smye, Hennings, & Horne, 2021).Large variations in the azimuth of the maximum horizontal compressive stress (S Hmax ) have been documented, with rotations from roughly north-south in the northern DB (NM) to east-southeast-west-southwest in the south (Lund Snee & Zoback, 2018), likely due to elevated pore pressures and small differences between the horizontal stresses, as has been observed in other regions (e.g., Moos & Zoback, 1993), as well as structural complexity related to late Paleozoic tectonism from the combined Ancestral Rocky Mountain and Ouachita-Marathon orogenic events (Horne et al., 2021).Morris et al. (2021) used well density data (S v ), borehole breakouts (S Hmax azimuth), qualitative estimates of faulting stress regime from available fault slip and focal mechanism data (A φ ) to provide representative stress tensors for 20 sub-regions of the DB and provided an uncertainty analysis as well as a most-likely stress tensor for each stress domain to determine stress criticality of basement-rooted faults.

Hydrogeology
The general features of groundwater flow in the DB under natural conditions are well-known, and can be described by five flow units, the deepest of which is the topic of this study (fifth unit) (Figure 2).The first unit is composed of fresh or brackish water aquifers hosted by a veneer of mostly post-Permian sediments overlying an evaporite succession of Late Permian age (second unit) that forms a thick confining layer.The third flow unit consists of the siliciclastic Delaware Mountain Group that has been evaluated previously as the shallow injection interval in the DB (Ge et al., 2022).It overlies the thick producing units of the Bone Spring and Wolfcamp Formation (fourth unit), which forms another confining system containing low-permeability shales.All five units have received secular recharge on their western and/or northern boundaries and the general flow of the entire system is to the east/southeast, following the general topographic slope (Engle et al., 2016).However, if the first flow unit is dynamic at the human scale, the aqueous phase of the other units can be considered stagnant, especially when compared to and dominated by injection-driven flow.

Geologic Characterization
Characterization of the stratigraphic architecture and reservoir properties of deep injection intervals provides the basis for hydrogeologic modeling.While the model domain is restricted to north of the GFZ, geologic characterization extends beyond the model boundaries to provide control on the stratigraphic and lithologic character of injection formations.

Stratigraphic Interpretation
The stratigraphic architecture of deep injection intervals was interpreted through correlation of formation tops in 3,151 wells that penetrate sub-Wolfcamp strata and have digital geophysical logs for further analysis (Figure S2 in Supporting Information S1).Formation tops from Geologic Data Services (IHS Energy, 2009) and those identified by operators and provided to commercial databases (S & P Global, 2023) were quality checked and included in the data set.Key units were correlated using a lithostratigraphic approach, meaning that in the absence of major lithologic changes in the sedimentary section, formations were grouped into "disposal units."The top of the model boundary is the base of the Pennsylvanian-age interval of the Wolfcamp Formation (base of Wolfcamp D) (Figure 2).The units underlying the Wolfcamp D are Pennsylvanian-age carbonates of the Strawn Group and Atoka Formation in the northern DB study area, which grade into mudrocks in the central DB.The Pennsylvanian-age Strawn Group, Atoka Formation and Morrow Formation were not differentiated, and are treated as one Pennsylvanian-age disposal unit.
Below the Pennsylvanian-age formations, the Barnett Formation was correlated with 1,636 wells.The Barnett Formation is interpreted by an increased gamma ray (GR) log signature (>100 API) indicative of the presence of clay in the dominantly mudrock lithology.The contact between the Barnett Formation and the overlying Morrow Formation is gradational, so correlation also utilizes the density and neutron porosity logs, with decreased density indicative of increased total organic carbon and porosity, both of which decrease the bulk density log response.
Underlying the Barnett Formation is a thin limestone of lower Mississippian age, the Mississippian Limestone, which is identified by its low GR and high resistivity log signature.Underlying the Mississippian Limestone, the Woodford Formation is interpreted by its very high GR log response, frequently exceeding the recorded values in well logs, associated with high clay and total organic carbon volumes.While the Woodford Formation was identified in the stratigraphic analysis, reservoir characterization was limited to fewer wells due to poor log quality in this interval, with lithologic interpretation requiring repair of the GR log response using the density log.
Main injection units of Silurian and Devonian age underlie the Woodford Formation and are correlated by the sharp decrease in GR below the base of the high-gamma ray Woodford in 1786 wells.Because these formations are challenging to differentiate with petrophysical logs alone, and because of the absence of clear lithologic variation and flow barriers, the Montoya Formation, Devonian Thirtyone Formation, and Fusselman Formation were interpreted as one Silurian-Devonian disposal unit (Figure 2).This approach is consistent with previous work suggesting that this dominantly carbonate section was deposited consistently without interruption apart from minor changes in sea level (Frenzel et al., 1998).
The base of the Silurian-Devonian disposal unit is the Simpson Group, which is identified in 952 wells by its increased GR and porosity log responses compared to over-and underlying formations.The contact between the Simpson and Ellenburger Groups was identified by a decrease in GR and increase in resistivity log responses in 1,105 wells that penetrate the Ellenburger Group.However, most of these wells are located where the formation is shallowest along the CBP margin.Although the lithology of the crystalline basement and the nature of the basement-sediment interface are likely important for understanding the connection between injection and basement-hosted earthquakes, characterization of the basement was not possible given few (∼35) basement well penetrations in the study area, most of which are lacking robust log suites for petrophysical characterization.

Reservoir Property Interpretation
Initial flow properties of the hydrogeologic model include porosity and permeability, which were interpreted by analysis of subsurface petrophysical log data for 46 wells that record the deep injection intervals of interest (Figure S2 in Supporting Information S1).Initial processing included conversion of logs reflecting porosity, that is, those measuring bulk density, neutron porosity, and compressional sonic transit time, to limestone porosity units using the Basic Log Function module of Interactive Petrophysics ® (IP) (2021) software.Photoelectric curves were combined with bulk density curves to derive the cross-section/cc.Initial estimates of clay volume from differences in neutron and density porosity logs, and from rescaled GR logs were computed and compared.
The Mineral Solver tool in IP® was used to compute mineral volumes and porosity.The optimized solver approach allows for solving a linear system of equations while minimizing the error between input log variables and the forward-modeled version of these variables computed from the solver results.Input logs included density and neutron porosity, GR, curves derived from photoelectric effect logs, and clay volume derived from neutron porosity logs.In the absence of a photoelectric log, which is key for differentiating limestone from dolomite, a log proxy N (Burke et al., 1969) was computed from porosity-independent crossplots of density and neutron porosity, where N is the slope of a line connecting 0% and 100% porosity end points.This approach has been useful in other carbonate-dominated systems for identification of lithology (e.g., Smye et al., 2019 for the Ellenburger Group of the Fort Worth Basin), and was calibrated to mineral solver model results with N ≤ 0.55 for dolomite and N > 0.55 for limestone, respectively, to reflect dolomitic and calcitic mineralogy.Volumes of quartz, calcite, clay, dolomite and porosity were computed in the optimized log solver.Where sonic logs were available, the sonic log response in the Ellenberger Group carbonates was used to aid in validating relative amounts of calcite and dolomite computed by the solver.
Multi-formulae scripts were used to estimate the mechanical properties of lithologies encountered by the well while accounting for the impact of overpressure (Eastwood & Smye, 2021), which may exist in clay-rich Pennsylvanian and Mississippian strata.Estimation of mechanical properties begins with the M elasticity modulus, which requires bulk density and compressional transit time logs, from which Young's Modulus (E) and bulk modulus (K) are computed by use of Poisson's ratio (ν) in well-known relationships between elastic moduli (Zoback, 2007).The algorithm for computation of ν may underestimate ν in lithologies that exhibit total porosity greater than 0.08 to 0.10 v/v, such as sandstones of the Simpson Group.
Mineral volumes interpreted by petrophysical analysis were grouped into rock types or lithologies based on dominant mineralogy.Four generalized rock types were interpreted: limestone (calcite-dominated), dolostone (dolomite-dominated), siliciclastic (quartz-dominated) and mudrock (clay-dominated) for a total of 120 wells (Figure S2 in Supporting Information S1).For wells in which mineral volumes could not be computed using an optimized mineral solver model due to an under constrained system, lithology was assigned using simplified log models (Table S1 in Supporting Information S1) calibrated to results for wells with mineral solver model results.

3D Geologic Model
We use the correlation of lithologic and stratigraphic layers, general structural control, and petrophysical analysis to develop a 3D geological model of the northern DB, within which reservoir properties are distributed.The geomodel spans basement through Pennsylvanian strata and encompasses all or parts of four counties in TX and two (Eddy and Lea) in southeast NM, covering an areal extent of roughly 28,000 km 2 (10,800 mi 2 ).The model is bounded to the east by the western structural boundary of the CBP and to the south by the GFZ, both of which nearly or completely offset the intervals of interest.These boundaries are assumed to be no-flow.The Northwestern Shelf and Diablo Plateau (Figure 1) define the northern and western boundaries of the geomodel, respectively; while no-flow conditions may not be strictly applicable to the northern and western boundaries where there is no evidence of faults offsetting the intervals of interest, these areas are distal from operational activity, specifically deep injection.
The geomodel was constructed using PETREL (v.2019.4) with a horizontal cell size of 500 m × 500 m (1,640 × 1,640 ft), and an average vertical cell size of 5 m (10 ft), totaling 137 × 10 6 cells.Of the wells utilized in the regional stratigraphic interpretation, 1,513 wells are contained in the study area and used in the geomodel, 35 of which penetrate the Precambrian basement.The model is compartmentalized into eight lithostratigraphic zones, with a lower model boundary defined at 600 m below the basement-sediment interface, assuming a permeability decrease associated with crystalline and metamorphic basement examples observed elsewhere (Gao et al., 2021).To reproduce the stratal layering geometry and maintain documented lateral fluid continuity among zones, we used an onlap grid for the Ellenburger, Simpson, and Woodford zones, a truncation grid for the Silurian-Devonian zone, and proportional zones for the Base Wolfcamp, Barnett, Mississippian Limestone, and Precambrian base zones.
Although basement-rooted reverse and shallow extensional faults are present in the study area (Horne et al., 2021(Horne et al., , 2022)), they were not modeled explicitly.Apart from those faults bounding the model, that is, the GFZ on the southern boundary, their throws are minor compared to the thickness of the modeled units, and their properties are uncertain.For the purpose of this model, the top (base Wolfcamp) boundary is assumed to be no flow because of the significant contrast in matrix permeability, as evidenced by generation and millions-of-years retention of overpressure in the overlying shale intervals (e.g., Luo et al., 1994;Sinclair, 2007).
Lithofacies, porosity and permeability were upscaled and distributed in the geomodel.The vertical resolution of these properties is 15 cm in each well.The lithofacies model was built on a by-facies basis one zone at a time, allowing separate retention of data for each facies in each zone.The first step in the reservoir modeling process was to create a generic subdivision of the reservoir based on five generalized rock types (shale, siliciclastic, limestone, dolomite, and granite for Precambrian basement).This lithology classification allows representation of the reservoir heterogeneity and calculation of statistics for reservoir modeling using the stochastic sequential indicator simulation, which allows for reservoir properties to be distributed within the constructed 3D structural grid.Lithofacies were key in conditioning the matrix porosity model, which was simulated using the Gaussian random function.Porosity-permeability transforms were used as a first approach to model the matrix permeability field by assuming similar rock fabrics and utilizing core petrophysical values in Ruppel (2019b) and Ruppel, Hovorka, and Barnaby (2019) in porosity-permeability transforms from Lucia (1993).

Model Setup
The modeling package CMG-STARS v.2019.10(Computer Modelling Group CMG, 2019) was used to evaluate the pore pressure history.The flow model has 20 layers, from Pennsylvanian to Precambrian basement, constrained by the stratigraphic interpretation.The model corner-point grid is 270 × 220 km 2 (168 × 137 mi 2 ) with 123 × 116 × 20 cells of approximately 1.6 × 1.6 km 2 (1 × 1 mi 2 ) and includes 285,360 active cells.The average cell thickness is 103 m (337 ft) and the thickness of 85% of the cells around the median ranges from 24 to 183 m (80 to 600 ft).Grid cells are aligned in the N-S direction.Fluid properties such as viscosity and density were computed as a function of temperature, pressure, and salinity.
Data including monthly injection volumes, average wellhead pressures, bottomhole pressure (BHP) estimated from surface pressure, and temperature data were obtained from the Railroad Commission of Texas (RRC) H-10 forms, and from commercial databases including S&P Global IHS Markit and B3 Insight.The data were qualitycontrolled based on previous experience (Gao et al., 2021;Ge et al., 2022;Lemons et al., 2019).Physical characteristics of the SWD wells, such as perforated intervals, open hole or cased status, and tubing dimensions were sourced from S&P Global IHS Markit.
Porosity and permeability derived from petrophysical analyses and distributed in the geomodel are best understood as a representation of the rock matrix properties, while dynamic fluid flow in karsted carbonate systems is mainly a function of fracture and fault systems.Populated and distributed permeability data from log and core data were used as a "permeability index" and supplemented with larger-scale permeability values estimated from injection rate and SWD well head pressures (Figures S3 and S4 in Supporting Information S1) and step-rate tests obtained from the RRC.
These field permeability values best represent the regional permeability that accommodates fluid injection (Figure S5 in Supporting Information S1).Injectivity-derived permeability was estimated from the quality controlled history injection data from SWD wells by assuming steady-state flow conditions, constant water viscosity of 1 cP (1 kPa • s), constant wellbore radius of 0.076 m (0.25 ft) and constant effective radius of 304.8 m (1,000 ft).These calculated permeability values were added to the flow model at their original SWD well locations, with the assumption that SWD wells are usually located in higher permeability zones.Injectivity-derived permeability locations are more numerous than log-derived permeability locations and were critical in areas with a lack of petrophysical data (Ge et al., 2022).
Injectivity-derived permeability from 175 data points shows a median of 11.2 md and mean of 21.8 md (Figure S3a in Supporting Information S1).These data were complemented by permeability values from 35 field data points taken from historical reservoir records that show a median of 10.0 md and mean of 40.6 md (Figure S3b in Supporting Information S1), resulting in combined permeability of 10.6 and 24.0 md for median and mean, respectively (Figure S3c in Supporting Information S1).Field and injectivity-derived permeability data vary by formation but are most abundant in the Silurian-Devonian intervals where many of the SWD wells are located (Figure S4 in Supporting Information S1).
The field and injectivity-derived permeability, together with the log-derived permeability data points, were combined as control points in the geomodel to create upscaled well logs and generate variograms.The variograms were used in the petrophysical modeling process in the geomodel to generate realizations of the 3D permeability field using a sequential Gaussian simulation method, as is common in reservoir characterization studies, with a lithofacies-specific permeability realization.This uncalibrated horizontal permeability (Kh) field was then exported to the flow model to be optimized during the calibration process (Section 3.3.3).Vertical permeability values were assumed to have a constant Kv/Kh ratio per layer, with Kv/Kh set to 0.1 (base case).However, in clay-rich layers including the Woodford and Barnett, a Kv/Kh of 0.01 was used to account for the low vertical permeability conditions.
Compressibility is the other parameter controlling diffusivity in addition to permeability.Compressibility values were computed using sonic and density logs.Initially, the M elasticity modulus was computed as where V P is the compressional wave velocity and ρ is density, or ρ/ΔT C 2 , where ΔT C is the compressional wave transit time.The shear modulus, G, is then computed by G = M (1 2ν)/(2 2ν), where ν is Poisson's ratio.The bulk modulus K can be calculated by K = M 4 G/3, and is the inverse of compressibility.Log-derived compressibility values range from 4 to 27.6 × 10 6 (1/psi).Model parameters used here are comparable to those used for modeling of other deep carbonate injection systems in other basins, for example, Gao et al. (2021).

Initial Conditions
Initial pressure and temperature for the flow model were estimated using initial shut-in pressure measurements obtained from drill-stem tests for sub-Wolfcamp formations.The pressure data points (n = 1,034) are consistent with an initial pressure gradient of 0.46 psi/ft (Figure S7a in Supporting Information S1) and temperature data points (n = 207) and are consistent with a temperature gradient of 1.05 F/100 ft (Figure S7b in Supporting Information S1) (Table S2 in Supporting Information S1), in line with values noted from well headers during petrophysical analysis.When assigned to injection stratigraphy, initial pressure gradients are 0.479 psi/ft (n = 472) for Pennsylvanian and Mississippian formations, 0.457 psi/ft (n = 227) for Silurian and Devonian, and 0.463 psi/ft (n = 124) for Ordovician.However, prior to injection, the flow model is at equilibrium (no flow, density-stratified), which is reached by running a long "multi-century" transient simulation before injection begins (Ge et al., 2022).
Formation water salinity ranges from 20,000 to 200,000 ppm; pore water is relatively fresh toward the west of the DB, where the sedimentary cover is thinner and secular recharge occurs, but it increases to high values against the CBP and with depth (Blondes et al., 2018).Injection water salinity was estimated to have a constant value of 120,000 ppm as measured from produced water total dissolved solids (TDS) measurements; however, the TDS is known to vary across the Basin (Nicot et al., 2020).Water is modeled as a mixture of two geochemical components: saline water of 20,000 ppm at the shallowest depth and saline water of 200,000 ppm at the deepest part of the model.The initial salinity field is modeled as containing variable fractions of these two components.

Model Calibration
Calibration of the flow model is carried out globally (across the entire model) and relies on analyzing wellhead injection pressures and converting them to equivalent BHPs.The wellhead flowing pressures range from ∼1 psi to >3,500 psi (0-24.1 MPa), with an average of ∼1,100 psi (7.6 MPa).Surface pressures were converted to equivalent BHPs using an algorithm by considering fraction losses (0-1,543 psi), that is, the energy lost as the fluid travels through the wellbore due to friction between the fluid and pipe, which are estimated by assuming a relatively low absolute roughness coefficient of 1.3 × 10 5 in for the steel casing (Bourgoyne et al., 1986).Values Geochemistry, Geophysics, Geosystems 10.1029/2023GC011260 of all or large sections of the permeability field were globally altered by iterating with the geomodel until an acceptable match with the 2,873 BHP data points was found (Figure S8 in Supporting Information S1).

Sensitivity Studies and Scenario Testing
Results for selected sensitivity cases are reported as the change in differential pressure from model base case (ΔP p ) results and Root-Mean Square error (RMSE, Table S3 in Supporting Information S1).The sensitivity of pore pressure changes to model parameters was investigated by scaling porosity, horizontal permeability, vertical permeability, and compressibility up and down by a factor of two, and by changing the relative anisotropy in horizontal permeability.It is well established, especially for stiff and fractured reservoir rock, that permeability can be enhanced in the direction parallel to the azimuth of maximum horizontal stress (S Hmax ) in a normal faulting regime such as the DB (Barton et al., 1995;Hennings et al., 2012).We therefore include consideration of horizontal permeability anisotropy in our uncertainty analysis by keeping the overall permeability constant but varying its directional proportions.The model was split into two areas: one in which S Hmax orientations are largely in the N-S direction, and one in which S Hmax is closer to E-W (Figure 1) (Dvory & Zoback, 2021;Lund Snee & Zoback, 2018;Lundstern & Zoback, 2020).Two scenarios were tested: in the first scenario, KI = 2Ko, KJ = 0.5Ko where S Hmax azimuth is N-S, and KI = 0.5Ko, KJ = 2Ko where S Hmax is E-W, with KI = permeability in the x direction, KJ = permeability in the y direction, and Ko = original permeability field.In the second scenario, KI = 4/3Ko, KJ = 3/4Ko where S Hmax azimuth is N-S, and KI = 3/4Ko and KJ = 4/3Ko where S Hmax is E-W.
Comparisons of flow model results and the ΔP p required to reach slip criticality were performed for a subset of the basement-rooted faults in the CMEZ region (Figure 1) that have hosted recent seismicity.Three parameters were modified to achieve a match between modeled ΔP p and slip criticality: local increase in permeability, increase in overall injection rate and volumes through addition of injection wells, and adjustments in rock compressibility (Table S2 in Supporting Information S1).

Analysis of Seismicity and ΔP p
To enable quantitative comparison of ΔP p and earthquakes, an analysis of trends in seismicity in the northern DB region was carried out.Earthquakes occurring in the study area consist of 1,559 from NMT Seismological Observatory with local magnitudes (M L ) ranging from 0.3 to 4.03, and 9,098 from TexNet of M L of 0.1-5.4 from January 2017 to December 2022, with a magnitude of completeness for the TexNet catalog of M L 1.55 (Huang et al., 2022).Due to high uncertainty in earthquake hypocentral depths (averaging 1.5 km or 4,900 ft for TexNet earthquakes analyzed), the vertical dimension was not considered for sampling of ΔP p .However, earthquake relative depths within clusters were compared because the events are spatiotemporally limited.
Pore pressure models were sampled at earthquake epicenters using Python 3 and Jupyter Notebook libraries to extract and refine high-resolution pressure data from the reservoir simulation output generated in CMG.The middle layer of the Silurian-Devonian injection interval was selected as representative of "deep injection strata" because it is the zone most frequently injected into and represents the approximate maximum ΔP p on basementrooted faults extending into the deep injection interval.The extraction algorithm interpolates between time steps prior to and after earthquake events to deliver a pressure matrix corresponding to each earthquake.Earthquake location-specific ΔP p is achieved using a K-Nearest Neighbor approach, with pressure data interpolated using an inverse distance-weighted algorithm.
Clustering analysis was performed to identify groups of earthquakes that may represent individual faults and to enable characterization of the spatiotemporal evolution of earthquake clusters in the presence of increased ΔP p .Earthquakes were clustered using the density-based spatial clustering of applications with noise (DBSCAN) algorithm (Ester et al., 1996) in MATLAB.This algorithm partitions observations into clusters based on a threshold for a neighborhood search radius and the minimum number of neighbors required to identify a core point and is frequently used in the analysis of earthquakes to assess spatiotemporal evolution (e.g., Schoenball & Ellsworth, 2017;Skoumal et al., 2019).Time and depth were not used as input in the clustering analysis to enable assessment of behavior of earthquake clusters over time with pore pressure evolution, and because the uncertainty in depth is the greatest among spatial locations.Because two earthquake catalogs were used with varying event densities, the algorithm was applied at two levels of spatial resolution using parameters based on estimated epsilon nearest neighbor analysis for each earthquake catalog.Regional clusters were identified using a neighborhood search radius of 2,000 m and at least 10 earthquakes required to identify a cluster.The search radius and minimum number of earthquakes required to identify a cluster were modified to 500 m and 20, respectively, for the CMEZ due to higher event density.
Regional clustering of earthquakes in the study resulted in the identification of 17 clusters, with an average of 85 earthquakes per cluster and 988 unclustered events (Figure S9 in Supporting Information S1).Most events in the NMT catalog are unclustered due to the distance between earthquakes as currently located in horizontal space.In the regional clustering approach, earthquakes in the CMEZ form one large cluster due to the close spacing of events.Decreasing the neighborhood search radius resulted in the identification of 30 clusters in this sub-region, with an average of 237 earthquakes per cluster and 2,243 unclustered events.

Assessing Stability of Recently Seismogenic Basement-Rooted Faults
We quantified the likely ΔP p required to push faults to criticality and characterized their native stability in a static stress analysis using the 3DStress® software tool (https://www.swri.org/3dstress-geology-software).Previously, the stability of mapped faults in the DB was characterized using the most likely stress tensor in the relevant stress domain (Morris et al., 2021).Here, we map and study five representative basement-rooted fault segment traces that are known to host large-magnitude events and have a range of fault orientations, and we characterize the range of stability given all possible stress states in each region.Utilizing recent earthquake hypocentral locations and fault plane solutions along with data described in Horne et al. (2021), basement-rooted fault segments were mapped as traces representing the approximate hanging wall intersection between the basement-rooted fault and the top of the Ellenburger Group.These map-view fault traces were extruded to fault planes by applying an average dip angle normal to the trend of each fault line following the methodology of Horne et al., 2022.Geologically native stability of these four faults can be described by slip tendency (T s , or the ratio of resolved shear to normal stress on any plane surface as defined by Morris et al., 1996), which represents the likelihood that a fault will slip under a given set of conditions.Each fault is treated as a set of triangular facets to capture variability in slip tendency over its area, as faults are not generally planar.We use triangulated surfaces of the faults because, although not unique 3D surface representations, they are an efficient and widely used method for such representations.We opted to generate non-planar surfaces to capture whatever non-planar fault geometry exists, which is an important component in determining how, where, and over what rupture areas seismic events may occur.Areas of the fault with high T s will be more likely to slip, as will faults with a higher average T s .Fault stability is also characterized by the pore pressure from the ambient state that will cause a fault surface to reach criticality (ΔP fc ) (Morris et al., 2016(Morris et al., , 2021;;Zoback, 2007), which is equal to the critical value of friction.Here, we assume a critical friction of 0.6, providing a consistent representation of relative fault stability and allowing comparison to modeled pore pressure changes with deep injection.

Geology of Deep Injection Strata
More than 80% of deep injection in the DB is into Late Ordovician, Silurian, and Devonian age strata including the Montoya, Fusselman, Wristen, and Thirtyone Formations; this injection interval is bounded by the lowpermeability Woodford Formation at its top and the Simpson Group at its base.Elevation at the top of the deep injection interval ranges from 0.5 km above MSL in the western DB and on the CBP to more than 5 km below MSL in the basin center, deepest proximal to the CBP (Figure 5a).The present-day structure of these units reflects basement-rooted reverse faults bounding the DB to the east and offsetting the injection units, and the eastwest trending GFZ at the southern boundary of the model domain.
The thickness of this deep (Siluro-Devonian) injection unit ranges from 75 m (250 ft) on the western and southern boundaries of the study area to ∼700 m (2,250 ft) toward the northeast in Lea Co., NM (Figure 5b).The average thickness in the study area is 380 m (1,250 ft) and is a product of localized erosion in the northern part of the CBP.Apart from this localized thinning, the greatest depositional thickness of the injection units is north of the Silurianage platform margin, as inferred from thickness trends and facies distributions.Northwest of this inferred margin and covering the modeled area, dolomite is the dominant mineral, averaging from 50 to 72 vol.% in the logged deep (Silurian-Devonian) section (Figures 6 and 7a).Dolomite transitions to limestone toward the southeast, where the dolomite fraction is less than ∼20 vol.% on average.Locally, quartz volume is up to 50 vol.%,particularly in the southern DB, and is related to chert deposition in deep water off the Silurian platform margin.
Of secondary importance, the Ordovician-age injection interval consisting of the Ellenburger and Simpson Groups (Table S4 in Supporting Information S1) is defined at its base by the top of the Precambrian basement, which ranges from 1.5 km below MSL to the north and west to more than 6 km below MSL to the east, with similar structural trends to those observed for the Siluro-Devonian section.Thickness of this injection interval ranges from 150 m (500 ft) in the northwest to more than 1,000 m (3,000 ft) in the southeast (Figure S10a in Supporting Information S1), with a regional average thickness of 560 m (1,860 ft), and an average thickness of 290 m (960 ft) in the study area.Ellenburger Group mineralogy is dominantly dolomite and limestone, with the dolomite fraction being the most abundant (Figures 6 and 7a).The overlying Simpson Group contains increased clastic material, frequently having >50 vol.% quartz, with subordinate clay and limestone (Figure 7a).Because very few wells log the complete stratigraphic section from Precambrian basement to Simpson Group, characterizing the spatial trends in lithology of the Ordovician disposal unit is challenging and mineral proportions averaged by well are not predictable across the basin.
Injection into Pennsylvanian-age units accounts for 16% of deep injection in the DB.The depth of Pennsylvanianage units ranges from 600 m (1,900 ft) below MSL to the west in Culberson Co., TX and Eddy Co., NM to deeper than 3,000 m (9,800 ft) below MSL in Loving Co., TX (Figure S11 in Supporting Information S1).The Pennsylvanian disposal is 150 m (500 ft) thick toward the south and north to more than 1,000 m (3,500 ft) to the east (Figure S10b in Supporting Information S1), and averages 520 m (1,700 ft) in the study area.The Pennsylvanian disposal unit is a mixed carbonate-siliciclastic system and is calcite-dominated throughout most of the study area, with less abundant quartz and clay transitioning to quartz-and clay-dominated in the southwest.
Porosity and permeability values computed here are low-representing the rock matrix-compared to fractureenhanced permeability values in measured data.Matrix porosity values as interpreted by petrophysical modeling are uniformly less than 5 vol% in carbonate-dominated units including those of Ordovician, Silurian and Devonian age (Figure 7b, Figure S12a in Supporting Information S1).Porosity is slightly elevated where dolomite is the dominant mineral.Fractures are evidenced by infrequent spikes in petrophysical logs, including in the caliper-measuring wellbore diameter-and density logs.Fracture-enhanced permeability is observed at low porosity values based on an analysis of published core data combined from Ruppel (2019b), and Ruppel, Hovorka, and Barnaby (2019) for Silurian and Devonian rocks in the PB, with porosity of <5% corresponding to permeability values of 10 md or more (Figure S12b in Supporting Information S1).These permeability values are greater than can be predicted by relationships based on rock fabric classes (Lucia, 1993), even for grainstone and large-crystal dolostone (Class 1).Limestone and dolomite lithologies both exhibit increased permeability to differing degrees, highlighting the importance of lithology-controlled distribution of flow properties in the geomodel.Because log-based methods for determining porosity-and therefore permeability-do not reliably capture fracture porosity and permeability, flow model calibration using injectivity and field permeability is a key step in determining the overall permeability field.

Flow Modeling of Deep Injection
Results of flow modeling are shown as ΔP p from an initial state in January 1983.Although pressure can be extracted from any of the 20 layers in the model, we take the middle of the Silurian-Devonian injection unit as representative of the deep injection interval because in most cases it corresponds to the maximum ΔP p throughout the sedimentary section where fluid injection volumes are greatest, and has received more than 80% of the fluid injected into deep strata in the basin (Table S4 in Supporting Information S1).Pore pressure associated with injection into deep strata increased by 0-6.5 MPa (0-950 psi) across the model domain from January 1983 to January 2023 (Figure 8).A maximum annual increase in the rate of change in pore pressure was reached in 2021-2022 when ΔP p had increased by more than 3.4 MPa (500 psi) (Figures 9a and 9b).In the main deep injection interval, pore pressure increased by an average of 0.6 MPa (87 psi) as of the beginning of 2023, with average annual increases of 0.07 MPa (∼10 psi) per year in 2019-2021, and an average pressure increase of 0.2 MPa (30 psi) in 2022 (Figure 9a).Pressure increase is greatest proximal to injection wells, with smaller perturbations away from injection centers due to pressure diffusion.
Pressure buildup also shows high vertical variance (Figures S13a-S13c in Supporting Information S1) due to preferential injection into formations with greater injection capacity and better reservoir quality, and due to pressure diffusion being greatest in the downward direction.For deeper strata, for example, the Ordovician-age Ellenburger Group directly above basement, ΔP p ranges from 0 to 5.5 MPa (0-805 psi), averaging 0.55 MPa (80 psi) (Figure S13a in Supporting Information S1), which is on the same order as that of the main Silurian-Devonian injection interval (Figure S13b in Supporting Information S1).Up-section, ΔP p in the Pennsylvanian-age formations ranges from 0 to 1.85 MPa (0-268 psi), averaging 0.34 MPa (50 psi) (Figure S13c in Supporting Information S1), which is 40% lower than in deeper strata.This is likely related to low Kv/Kh ratios modeled for Devonian and Mississippian shale strata, which separate injection intervals from overlying Pennsylvania strata, and the degree of lithologic similarity and hydrologic connectivity between carbonates of Ordovician, Silurian and Devonian age.Results of the sensitivity analysis (described in Text S1, Figures S14a-S14h, Table S3 in Supporting Information S1) show that horizontal permeability has the largest impact on pressure buildup, with a larger RMSE than other parameters.
History matching of modeled BHP and BHP estimated from field data for a few wells (Figure S15a-S15d in Supporting Information S1) shows a good match (calculated RMSE ranging from 2.8 to 4.1 MPa, 400-600 psi).Cases run to assess the impact of anisotropy in the permeability field aligned with S Hmax show only a modest difference in ΔP p trends (Figure S16 in Supporting Information S1).Cases run to attempt to push faults to criticality based on fault stability analysis included doubling all current injection volumes and decreasing compressibility by an order of magnitude (Table S3 in Supporting Information S1); these resulted in increased pore pressure along faults of interest but not of the magnitude to push them to slip based on current characterization of faults and stress.
While little downhole pressure data is available to assess the accuracy of the model, BHP measurements in the Eddy-Lea County line area of NM track an increase in pore pressure gradient of 0.68 kPa/m (0.03 psi/ft) from 2018 to 2022, which is equivalent to 3.45 MPa (500 psi) at average Silurian-Devonian injection depths of ∼5,300 m (17,500 ft) in this region (Figure S17a in Supporting Information S1).Although comparison of modeled ΔP p to single BHP measurements is challenging because BHP measurements are taken for one location at one point in time, the trend of increasing regional modeled ΔP p of several hundred psi (Figure S17b in Supporting Information S1) is consistent with an increase in BHP measurements in the northern DB.

Pore Pressure and Seismicity
Sampling of ΔP p at each earthquake epicenter and time in the deep injection interval (Figures 9d and 9e) shows that earthquakes proximal to deep injection, that is, in NM, occur at an average modeled ΔP p of 0.80 MPa (116 psi), and range from 0 to more than 4 MPa (0-640 psi) (Figure 9e).In contrast, earthquakes in the CMEZ that occur at greater distances from deep injection, including the largest magnitude events, were induced at modest pore pressure increases averaging 0.06 MPa (8.5 psi), with range from 0 to more than 3 MPa (0-462 psi) as currently modeled (Figure 9d).The maximum earthquake magnitude in this region increased over the study period, with the largest event to-date occurring in late 2022 (Figure 9c).
Modeled ΔP p sampled on five recently seismogenic faults (Figure 10a) shows a wide range of ΔP p associated with earthquake occurrence, from as little as 7 kPa (1 psi) to more than 1 MPa, indicating variability in stability among the faults and portions of faults hosting earthquakes in this region.Earthquakes occurring on two faults proximal to deep injection (NM1 and NM2) were initiated at modeled ΔP p of 0.2-1.7 MPa (30-230 psi), averaging 0.65 MPa (95 psi) for events cataloged by TexNet (Figures 10b and 10c).Fault NM1 experienced the greatest ΔP p , and at earthquake sequence onset it had experienced ΔP p of 0.78 MPa (113 psi) for events cataloged by NMT, and 0.87 MPa (127 psi) for events cataloged by TexNet (Figure 10b).Most seismicity on this fault, including the largest magnitude events, occurred several months after local maxima in ΔP p and quieted throughout 2022 as pressure along the fault subsided.This fault has experienced the greatest temporal variability in stressing rate of the five faults studied, likely due to its proximity to deep injection wells.Seismicity hosted on NM2 occurred coeval with that on NM1 at lower ΔP p following a slow and steady increase in pressure over the 3year period from 2017 to 2020, but occurred as the rate of change in ΔP p began to increase from 2020 to mid-2021 (Figure 10c).The onset of this sequence was initiated at ΔP p 0f 0.21 MPa (30 psi).
Faults CMEZ1, CMEZ2, and CMEZ3 are greater distances (11, 6, and 21 km, respectively) from deep injection and provide a comparison to earthquakes hosted on faults nearer to deep injection discussed above.Two early (2018) isolated events of M L < 1.5 occurred on CMEZ1 (Figure 10d), but the sequence began in earnest in 2020  after steadily increasing ΔP p , with little temporal variability in stressing rate, until it reached 0.05 MPa (7 psi).This fault hosted several episodes of seismicity associated with larger magnitude events, one in early 2020 related to the M L 4.9 "Mentone" event, and another in late 2021 associated with an M L 4.0 event.These temporally delineated episodes of seismicity consist of aftershocks associated with the larger magnitude events, and one possible explanation for their isolated occurrence is local variations in stressing rate associated with variability in injection rate from the nearest deep injection wells to the northwest (Figure 10a).Earthquakes on CMEZ2 occurred after a steady increase in ΔP p reached ∼0.3 MPa (40 psi) as sampled from deep injection strata at earthquake epicentral locations (0.37 MPa at mapped fault center as shown in Figure 10e) and continued at a low rate over several years after an initial several-month period of seismicity (Figure 10e).Seismicity on CMEZ3 (Figure 10f) began with an isolated pulse of earthquakes in mid-2018, increasing in frequency from 2020, and further increasing in 2022, resulting in the M L 5.4 "Coalson Draw" event, which is one of the largest earthquakes recorded in this region.Nearly 1,000 earthquakes occurred on faults in this local CMEZ3 area through the end of 2022, with magnitudes ranging from 0.9 to 5.4.Clustering analysis results in the identification of two clusters; however, improvements in earthquake catalog resolution and in the precision of locations would likely increase the number of distinct clusters.
Analysis of the spatiotemporal evolution of earthquake clusters shows that in some cases, migration of seismicity along strike (Figures 11a and 11b) and to greater depths (Figures 11c and 11d) occurs.In these cases, the earliest earthquakes occurred along portions of the fault that experienced the greatest ΔP p (Figure 10a).Notably, seismicity experienced on NM2 (Figures 11a and 11c) as cataloged by TexNet occurs in two spatio-temporally isolated phases.The earlier documented seismicity occurred on the NE portion of the fault, and later seismicity occurred on the southwestern portion of the fault.As presently mapped, this fault consists of one continuous segment with a slight jog in trend.This jog in the fault may represent a breached relay structure between two formerly isolated faults, which merged as deformation progressed during fault formation.In this case, the northeasterly portion of NM2 is activated first, followed by the southwesterly segment, and seismicity does not extend across the relay in either instance.In the case of CMEZ1, migration of seismicity to greater depths is consistent with the greatest ΔP p occuring where faults intersect the deep injection intervals and with pressure diffusion along the basement-rooted fault.

Fault Stability
Stability of the five representative faults as 3D surfaces calculated using all possible stress tensors applicable in their locations can be visualized as the cumulative fault area that has a given pore pressure required to push it to fault criticality (ΔP fc ) (Figure 12), with faults plotting toward the left having a greater area that may slip at lower pore pressure changes.The faults most proximal to deep injection (NM1 and NM2) on the Eddy and Lea County line, trend NE-SW and dip toward the SE (fault characteristics described in Table 1).The stress state in this area and its uncertainty are represented by nine possible stress tensors, which are points on a continuum of possible stress states (Figure S18 in Supporting Information S1, Table 1), with S Hmax azimuth ranging from 045°to 080°a nd A φ of 0.325-0.5 (Morris et al., 2021).The stability for both faults NM1 (Figure 12a) and NM2 (Figure 12b) has a narrow range even with uncertainty in the stress field A φ and the azimuth of S Hmax .As expected from earthquake occurrence and its experienced ΔP p , NM2 has a greater fault area that is unstable at lower ΔP p than NM1.Faults CMEZ1, CMEZ2, and CMEZ3 are located within the CMEZ.These faults have hosted earthquakes with recorded magnitudes ranging from 0.1 to 5.4, including the M L 4.9 "Mentone" event in March 2020 and the M L 5.4 "Coalson" event in November 2022.The stress state in this area has possible S Hmax azimuth ranging from 063°t o 103°and A φ of 0.55-0.6(Figure S19 in Supporting Information S1, Table 1) for CMEZ1 and CMEZ2, and S Hmax azimuth ranging from 085°to 130°and A φ of 0.55-0.65 for CMEZ3 (Figure S20 in Supporting Information S1, Table 1).Comparison of the calculated stability of CMEZ1 and CMEZ 2 (Figures 12c and 12d) with NM1 and NM2 shows a wider stability range with a strong dependence on uncertainty in S Hmax azimuth.Because the strikes of CMEZ1 and CMEZ2 differ at 064°and 101°, respectively, they are least stable under differing S Hmax azimuth conditions.Modeled ΔP p along both CMEZ1 and CMEZ2 compared to ΔP fc (Figures 12c and 12d) shows that either these faults are much less stable than currently understood or that the actual ΔP p is higher than  (Hennings, Nicot, et al., 2021).Sequences NM1 and NM2 are shown with earthquakes from both TexNet and New Mexico Tech earthquake catalogs without duplicate events removed.currently modeled.As currently interpreted, CMEZ3, the fault hosting the M L 5.4 "Coalson Draw" event appears to be stable in all possible stress tensors (Figure 12e) but nevertheless has slipped at modest modeled ΔP p .

Injection Geology
The pore pressure response of the strata we model is governed conceptually by three overlapping flow systems: (a) original matrix and matrix altered by diagenetic processes, (b) structural elements including fractures and faults, and (c) dynamic modification of the first two systems due to processes coupled to the pore pressure change.This third flow system includes poroelastic modification of the matrix and enhancement of permeability along fractures and faults, governed by the slip criticality of each, which may change dynamically with ΔP p .The pore pressure response of the matrix (flow system one) is controlled by distributions of thickness, lithology, and flow properties, which are a function of depositional environment and post-depositional diagenetic processes including exposure during relative sea-level lowstands, dolomitization, and karstification.Reservoir quality-namely fluid storage (porosity) and flow (permeability) properties-is directly related to the abundance of dolomite, with dolomite having increased porosity and permeability relative to limestone, as observed both in this study and others (Harrington & Ruppel, 2019;Ruppel, 2019b).This finding also parallels what Smye et al. (2019) found in their geological characterization of lower Paleozoic carbonate strata used for deep injection in the Fort Worth Basin.The dolomite fraction in Silurian and Devonian injection strata is greatest in the northern and northeastern DB; given its prevalence in shallower water platform areas, this is likely a function of more frequent exposure of the platform, and potentially prolonged (∼20 m.y.) exposure during the Middle Devonian (Entzminger & Loucks, 1992).
Silurian platform margin paleogeography influenced thickness and facies distributions from the Late Silurian to the Devonian.The position of this platform margin appears to be controlled by the Grenville orogenic front, where a basement domain boundary likely acted as a zone of crustal weakness controlling the subsidence rate.   1, and the stability of NM1 and NM2 was calculated using the stress tensors of Figure S18 in Supporting Information S1, CMEZ1 and CMEZ2 calculated using Figure S19 in Supporting Information S1, and CMEZ3 calculated using Figure S20 in Supporting Information S1.Shaded lines represent the fault stability for each stress tensor, with tensor 5 being the most likely stress state to be extant in the relevant domain.The range of modeled ΔP p experienced by each fault from 2017 to 2023 is shaded in gray for comparison to the ΔP fc .(Anderson, 1951;Simpson, 1997) (McKerrow, 1979;Ruppel, 2019a).
Although injection strategies are driven by operational considerations including proximity to produced water sources and acreage positions, deep injection in the DB occurs in areas where reservoir quality-namely, porosity and permeability-of these intervals is likely highest.However, characterization of matrix porosity and permeability is insufficient to fully understand the flow of injected fluid and associated perturbation of pore pressure because it does not account for the role of larger-scale fracture networks and transmission along faults, necessitating careful model calibration using measured pressure data from injection wells.Furthermore, the nature of the basement and basement-sediment interface is unknown in this region due to the lack of deep well penetrations; this adds uncertainty to attempts to characterize pressure diffusion in basal sedimentary and basement layers in the model.If, for example, a porous and permeable siliciclastic layer is present at the basement-sediment interface, it could provide a pathway for migration of fluids that is not represented in the geologic model presented here.
While conceptual flow system two (faults and fractures) is not directly addressable here because we lack data allowing for robust interpretation of fault systems, this enhanced permeability system is approximated in the flow modeling process, as discussed in Sections 3.3 and 4.2.We do not implement conceptual flow system three (dynamic modification of matrix and structural permeability due to pore pressure changes, including poroelastic stressing) here due to what we believe is a prevailing inability to parameterize the governing poroelastic parameters for the rock volume we study, which is characterized by heterogeneous and ubiquitous diagenetic modification, fractures and faults.However, our flow modeling process attempts to replicate effective permeability by calibration to field data.
Basin-scale discrete fracture and fault permeability enhancement to deep carbonates targeted for injection has been attempted in other studies (e.g., Gao et al., 2021 for the Ellenburger Group of the Fort Worth Basin).In that work, the authors found that varying the permeability tensor represented by the faults of semi-regional extent, even by a factor of 1,000, has a limited impact on the final differential pressures.This is because the permeability of large faults, estimated from local measurements and from a general understanding of karstic aquifers, is already very high (tens of darcys).Faults are an important factor in the flow system, but variations around the high fault permeability values do not control the dynamics of the regional flow system as long as the fault permeability stays high.However, the matrix permeability and the stochastic fault permeability, which includes the effect of fractures on the matrix through enhanced permeability, have a clear impact on the resulting differential pressures.The matrix and mesoscale fault horizontal permeability controls differential pressure and its highest values because injection wells are typically located away from major faults, as shown in this study.Based on this prior work, we conclude that it would be of value to implement stochastic discrete fracture and fault permeability systems in the flow model presented here when sufficient information is available to populate and calibrate the model.Then, the implementation of the dynamic coupling of evolving fault stress criticality and fault permeability could be considered.

Pore Pressure Evolution, Fault Stability, and Seismicity
Deep injection in the northern DB has increased the native pore pressure state by more than 10% locally (more than 5 MPa ΔP p in places), and has resulted in seismicity in two distinct geologic habitats, that is, in proximity to deep injection in NM and distal from (10s of km away from) injection in CMEZ.Rates of earthquakes of both >M L 2.0 and M L 3.0 increased steadily from 2017 through 2021 along with increases in pore pressure associated with deep injection.In 2022, while both injection and earthquake rates increased, the rate of increase slowed.Within the CMEZ region, the monthly rate of earthquakes is not directly related to the monthly rate of injection (Figure 13a), with monthly injection rates of 8 or more million barrels corresponding to the full range of monthly earthquake rates, from <50/month of M L 2.0+ to nearly 300.In this system, once pore pressure increases from regional deep injection were sufficient to cause fault rupture, earthquake rates remained elevated even when injection rates for the closest wells were decreased.In contrast, when considering all deep injection in the northern DB region, the monthly rate of injection and the monthly rate of earthquakes are exponentially related, with injection rates of 50 MMbbl/month and above corresponding to strong increases in earthquake rate.
Earthquakes proximal to deep injection, that is, in NM, have occurred at pore pressure increases that average 0.56 MPa (82 psi), and range from 0 to 4.1 MPa (0-600 psi) (Figure 9e).Most of these events occur in areas without mapped basement-rooted faults, but even events occurring early in the earthquake monitoring history (e.g., in 2017) occurred at ΔP p of up to 1.4 MPa (200 psi).The frequency of events in NM has increased with increasing regional pore pressure associated with deep injection and a greater geographic region impacted by ΔP p ; however, most earthquakes are still occurring at modeled ΔP p of 0-2.0 MPa (0-300 psi), which is evidence of the presence of basement-rooted faults in the impacted area that all likely have similar slip tendency and ΔP fc .These findings parallel those of Hennings, Nicot, et al. (2021), who found that seismogenic rupture of welloriented faults in the Fort Worth Basin in areas near concentrated injection occurred at ΔP p of 0.3 MPa.
In contrast, seismicity in the CMEZ region has occurred at modest pore pressure increases of <10 psi on average (Figure 9d).Most events occur in crystalline basement, with the average earthquake hypocentral depth at 3,350 m (11,000 ft) below the basement-sediment interface.As in NM, event frequency has increased in the CMEZ; however, more earthquakes are occurring at greater pore pressure increases compared to in early years (Figure 9d), while the overall geographic footprint of seismicity has not increased.This indicates that highly sensitive faults or portions of faults were pushed to criticality at modest (<0.1 MPa) pore pressure increases, and with continued pore pressure increase in recent years, some faults or portions of faults that were stable at ΔP p of In CMEZ, a decoupling of the rate of injection and the rate of earthquakes has occurred in recent months, as injection rates have decreased dramatically but rates of seismicity have been slower to fall.In the northern DB region, injection rates above 50 MMbbl/month are associated with sharply increasing earthquake rates.
<0.1 MPa have been made unstable at ΔP p of 1.4 MPa (200 psi).This can be observed in the evolution of earthquake sequences (Figure 11), which in some cases persist over the course of years with seismicity migrating along strike and deepening, assuming that earthquake relative depths are valid for sequences that are relatively spatiotemporally isolated.This finding also parallels Hennings, Nicot, et al. (2021), who found that seismogenic rupture of well-oriented faults in the Fort Worth Basin in areas distal to concentrated injection occurred at ΔP p of 0.05 MPa, and is consistent with the findings of Tung et al. (2021), who model pore pressure change of 154 kPa (22 psi) for the M5.0 "Mentone" event in 2020.
Many of the earthquakes in the CMEZ occur on mapped basement-rooted faults, for which native stability can be assessed (Morris et al., 2021) and compared to ΔP p .At ambient pore pressure, the least stable fault orientations are those with dips greater than 50°that are oriented optimally with respect to the most likely stress tensor field.
As currently characterized in terms of their geometry and coefficient of friction, faults in the CMEZ region require pore pressure increases to reach criticality that are on the order of at least 2× those currently modeled for deep injection (Figures 12c-12e).Either this fault system is much more sensitive than previously understood; ΔP p is underpredicted, as would be the case with permeability fields that are too low or compressibility fields that are too high; or stressing from poroelastic coupling serves to drive faults toward failure.
While poroelastic stressing, rather than direct pressure diffusion, is found to be a dominant mechanism for induced seismicity in some systems (Goebel & Brodsky, 2018), we find it unlikely that poroelastic coupling could play a dominant role in promoting fault rupture in CMEZ due to consideration of the location and distribution of ΔP p , the orientation of seismogenic faults within the in situ stress field, and the observed rupture style.As per Huang et al. (2022), the earthquake style in CMEZ is strongly dominated by normal-slip ruptures on steeply dipping planes with WNW strikes.These rupture planes parallel the interpreted faults in the CMEZ area and are subparallel to the local azimuth of S Hmax (Dvory & Zoback, 2021) (Figure 8).The dominant strike of seismogenic faults in the CMEZ is tangential to areas of significant ΔP p to the north northeast, and the azimuth of S hmin trends directly toward the areas of greatest nearby ΔP p .Therefore, poroelastic stressing from ΔP p is likely to cause an increase in the magnitude of S hmin , which would shift the faults away from Coulomb failure.This is consistent with the modeling work of Tung et al. (2021), who show that for the M5.0 "Mentone" event, deep injection serves to destabilize faults, and that the pore pressure component of the change in Coulomb failure stress is an order of magnitude larger than and opposite in sign to the poroelastic stress component.This proposition could be further tested numerically using the input data and results of our geomodel and flow model, as provided in the Supplementary Data.
Pore pressure model scenarios were run to assess whether ΔP p values could be made to approach ΔP fc for specific faults using any reasonable input parameters.These scenarios included local increases in permeability (×100, up to >40 darcys) trending toward CMEZ, increases in overall injection rate and volumes through addition of injection wells (doubling present injection), and decreasing rock compressibility (/10).While some areas of the CMEZ had a pressure buildup of ∼0.34 MPa (50 psi) as a result, none of these scenarios resulted in the ΔP p required to push faults to criticality given assumptions about fault geometry and stress state.
High permeability along critically stressed faults or fracture networks allowing for transmission of fluids or pressure (Townend & Zoback, 2000) to greater distance could also play a role in abundant seismicity in the CMEZ area.This phenomenon was proposed by Gao et al. (2021) and Hennings, Nicot, et al. (2021) as operative in the Fort Worth Basin allowing systems of critically stressed normal faults to transmit ΔP p up to 40 km from areas of significant injection to faults that became seismogenic.Given that the geomodel reflects the rock matrix, even when the permeability field in the flow model is adjusted globally using injectivity data, fault and fracture permeability conduits could locally play a role in far-field fluid and pressure migration.
Another possibility is that the faults in CMEZ are much more sensitive than previously understood, with lower frictional stability, and slip at ΔP p that is on the order of modeled results.Multiple lines of evidence indicate that observed earthquake sequences are frequently triggered by small perturbations to the Coulomb stress, on the order of tens of kilopascals (Keranen et al., 2014;Peterie et al., 2018;van der Elst et al., 2013;Yeck et al., 2016); these pore pressure increases are also on the order of that observed for triggering of more distal earthquake sequences in the Ellenburger of the Fort Worth Basin (Hennings, Nicot, et al., 2021) and within the range observed here.To assess the potential impact of frictional stability deviating from 0.6, we tested two cases: (a) ambient critical friction (μ c ), in which the principal faults and the host rock are considered to have the same value of friction (our default case is μ c = 0.6), and (b) fault friction (μ f ), in which we assume that the host rock has a friction value of 0.6, but the principal faults may experience a lower effective friction value (Figure S21 in Supporting Information S1).In case (a) the estimate of ambient friction affects the magnitude of the differential stress extant in the rock mass, decreasing the ambient friction can both increase (e.g., faults CMEZ 1 and 2) and decrease (e.g., faults NM1 and 2, and CMEZ3) the stability of faults, depending on the orientation of the fault within the stress field.However, if the principal faults experience effective friction that is lower than the surrounding host rock (nature of fault gouge, lack of asperities; Bedford et al., 2022), reductions in fault friction decrease ΔP fc values for all five faults considered here, rendering them less stable than predicted by applying the ambient friction to the whole system.Low-friction, principal faults may not reduce the rock mass strength by incremental slipping, as predicted by critical friction theory (Townend & Zoback, 2000), because of the presence of fault cohesion that has developed over the millions of years since the faults were last active.In this scenario, small events within the fault zone are likely to have been generated by early anthropogenic pore pressure increases localized on the fault, triggering small cohesion-reducing events.This would have the cumulative effect of destroying cohesion and linking slippable surfaces, culminating in larger events.
Notably, there are several regions (e.g., northern Culberson Co. TX, Eddy Co., NM) with mapped basementrooted faults that are not slipping seismically at elevated pore pressures of several These faults appear to either be not optimally oriented with respect to S Hmax for fault reactivation and/or are not as natively sensitive to rupture as compared to CMEZ and are therefore more stable with pore pressure increases associated with deep injection.

Summary
Deep injection in the DB primarily targets platform carbonates of Silurian and Devonian age that are up to 700 m (2,250 ft) thick with fluid storage and flow properties that are a function of lithology-specifically, proportions of limestone and dolomite.Lithology and facies distributions are controlled by the development of a Silurian platform margin, the location of which is consistent with a zone of crustal weakening along the Grenville front bisecting the DB.The strata subject to injection are highly complex-heterolithic, diagenetically modified, fractured and faulted.
Pore pressure increases with deep injection modeled using detailed injection geology ranging from 0 to 6.5 MPa (0 to 950 psi), with greatest increases proximal to injection centers.Across the model domain, the average pressure increase from January 1983 to January 2023 is 0.6 MPa (87 psi), with the average annual increase greatest (∼0.7 MPa or 10 psi per year) in 2019-2021.Sensitivity analyses show that the model is most sensitive to horizontal permeability, which is expected as this parameter controls pressure diffusion proximal to injection centers.
Earthquakes in NM have experienced pore pressure increases averaging ∼0.55 MPa (80 psi), as currently modeled.Earthquakes in CMEZ are initiated at much lower ΔP p (<0.1 MPa), indicating that either local anisotropy in permeability is present or that the seismogenic basement-rooted faults are much more sensitive than currently understood.The calculated stability of faults shows a strong dependence on their orientation relative to the azimuth of S Hmax , and faults hosting notable seismicity appear to be near-optimally oriented for slip.
The comprehensive models and data provided here can be used for further mechanistic studies of fault rupture in the northern DB and to understand other systems of seismicity induced by injection in the strata overlying seismogenic basement, such as those currently experienced in the Midland sub-basin of the greater PB.
Subsurface and operational data were obtained from S&P Global IHS Enerdeq and B3

Figure 1 .
Figure 1.Permian Basin regional map showing earthquakes cataloged by TexNet and New Mexico Tech (NMT) and sized and shaded by magnitude, basement-rooted faults in the basins(Horne et al., 2021(Horne et al., , 2024) ) and regional faults(Ewing, 1991), study area (gray dashed polygon), maximum horizontal stress orientations(Dvory &  Zoback, 2021), and Seismic Response Areas as determined by regulators in New Mexico and Texas.CMEZ = Culberson-Mentone Earthquake Zone; GFZ = Grisham Fault Zone; AOI = modeled area of interest, NDBEZ = Northern Delaware Basin Earthquake Zone (deep earthquakes in basement), SDBEZ = Southern Delaware Basin Earthquake Zone (some shallow earthquakes in the sedimentary section).Large-magnitude events discussed in text are noted.

Figure 3 .
Figure 3. Delaware Basin injection (SWD) wells for shallow (yellow) and deep (purple) sized by injected volume from 2015 to 2022.Area of most dense seismicity indicated by the North Culberson-Reeves Seismic Response Area (NCR SRA) 4.5 km radius; regional faults(Ewing, 1991), basement-rooted faults(Horne et al., 2021(Horne et al., , 2024) ) and shallow faults(Horne et al., 2022) are also shown.Inset shows monthly injection volumes into shallow and deep strata (B3 Insight) and monthly rate of earthquakes M L 2.0+ for the DB using TexNet earthquakes for Texas and New Mexico (NM) Tech earthquakes for NM.NDBEZ (AOI) = Northern DB Earthquake Zone and Area of Interest (this study); SCDBEZ = Southern DB Earthquake Zone.
Figure 3. Delaware Basin injection (SWD) wells for shallow (yellow) and deep (purple) sized by injected volume from 2015 to 2022.Area of most dense seismicity indicated by the North Culberson-Reeves Seismic Response Area (NCR SRA) 4.5 km radius; regional faults(Ewing, 1991), basement-rooted faults(Horne et al., 2021(Horne et al., , 2024) ) and shallow faults(Horne et al., 2022) are also shown.Inset shows monthly injection volumes into shallow and deep strata (B3 Insight) and monthly rate of earthquakes M L 2.0+ for the DB using TexNet earthquakes for Texas and New Mexico (NM) Tech earthquakes for NM.NDBEZ (AOI) = Northern DB Earthquake Zone and Area of Interest (this study); SCDBEZ = Southern DB Earthquake Zone.

Figure 4 .
Figure 4. Relative depths of injection and production wells and earthquakes from (a) north to south showing shallowing of seismicity in the southern Delaware Basin (DB) and (b) west to east showing extensive deep seismicity in the Culberson-Mentone earthquake zone region.Shallow injection (amber) and deep injection (purple) wells are sized by injection volume (B3 Insight), and hydraulically fractured wells (blue squares) (S&P Global) are represented by the horizontal midpoint.Earthquake hypocentral depths are from TexNet (solid points) and New Mexico Tech (open circles) earthquake catalogs.The average uncertainty in TexNet hypocentral depths is 1.5 km (4,900 ft).Geographic features shown in the inset map include NM = New Mexico, TX = Texas, CMEZ = Culberson-Mentone Earthquake Zone, GFZ = Grisham Fault Zone, SDBEZ = Southern DB Earthquake Zone, and CBP = Central Basin Platform.

Figure 5 .
Figure 5. (a) Elevation (m below mean sea level) of the deep injection interval (representing base of Woodford Formation and top of Silurian-Devonian); (b) Thickness (m) of deep (Siluro-Devonian) injection interval with a dashed line denoting location of interpreted platform margin, north of which deep injection strata are comprised mainly of dolomitized platform carbonates, and south of which coeval formations are mostly basinal limestone and chert.The inferred location of the Grenville Front (Mosher, 1998) is denoted as a solid black line.

Figure 6 .
Figure 6.Interpreted reservoir properties of deep strata for a well in Culberson Co., TX, showing raw input logs including gamma ray, deep and shallow laterologs (LLD and LLS, respectively), porosity and mineral volumes computed using a multi-mineral optimized solver approach, and dominant lithology.Deep injection mainly targets dolomites below the Woodford Formation.

Figure 8 .
Figure 8. Modeled change in pore pressure (ΔP p ) from January 1983 to January 2023 in the middle of the deep (Silurian-Devonian) injection interval shown with TexNet and New Mexico Tech earthquakes sized and shaded by magnitude, basement-rooted faults (modified from Horne et al., 2021), location of deep injection wells active anytime in the period 2010-2023, and S Hmax azimuth (Dvory & Zoback, 2021).The region shown in detail in Figure 10 inset is delineated by the dashed black line.

Figure 9 .
Figure 9. (a) Pore pressure increases in the model domain from 2010 to 2023 showing median ΔP p and 90th percentiles for all cells in the sampled deep injection interval (middle of the Silurian-Devonian), (b) Pore pressure increases in the model domain from 2017 to 2023 plotted with the rate of earthquakes of M L 3.0+ per year.(c) Earthquake occurrence and magnitude for events in the northern Delaware Basin from TexNet and New Mexico (NM) Tech earthquake catalogs, with notable events marked.(d) Modeled deep injection ΔP p sampled at earthquake epicenters for earthquakes in the Culberson-Mentone earthquake zone area (TexNet catalog).(e) Modeled deep injection ΔP p sampled at epicenters for earthquakes in the New Mexico (NMT).

Figure 10 .
Figure 10.(a) Five representative seismogenic faults (NM1, NM2, CMEZ1, CMEZ2, CMEZ3) shown with ΔP p in deep (Silurian-Devonian) injection interval as of January 2023 at the fault center.Temporal evolution of ΔP p experienced by each fault (left axis) shown with associated EQ sequences (right axis) and the modeled ΔP p at sequence onset denoted by dashed line for (b) NM1, (c) NM2, (d) CMEZ1, (e) CMEZ2, and (f) CMEZ3.Gray shaded area is ΔP p experienced by seismogenic basement-rooted faults in the Fort Worth Basin(Hennings, Nicot, et al., 2021).Sequences NM1 and NM2 are shown with earthquakes from both TexNet and New Mexico Tech earthquake catalogs without duplicate events removed.

Figure 11 .
Figure 11.Temporal evolution of seismicity in horizontal space (along strike) for (a) NM2 and (b) CMEZ1, and evolution of seismicity with depth for (c) NM2 and (d) CMEZ1.Earthquakes are shaded by date of occurrence and sized by magnitude and use the TexNet earthquake catalog, which reports depths for each event, enabling analysis of the temporal evolution of seismicity.

Figure 12 .
Figure 12.Fault stability comparison using nine possible stress tensors for the five representative faults (a) NM1, (b) NM2, (c) CMEZ1, (d) CMEZ2, and (e) CMEZ3.Plots show the cumulative fault area as a function of the change in pore pressure required to push it to criticality (ΔP fc ) with lower values of ΔP fc corresponding to less stable fault areas.The range of stress state inputs are noted in Table1, and the stability of NM1 and NM2 was calculated using the stress tensors of FigureS18in Supporting Information S1, CMEZ1 and CMEZ2 calculated using FigureS19in Supporting Information S1, and CMEZ3 calculated using FigureS20in Supporting Information S1.Shaded lines represent the fault stability for each stress tensor, with tensor 5 being the most likely stress state to be extant in the relevant domain.The range of modeled ΔP p experienced by each fault from 2017 to 2023 is shaded in gray for comparison to the ΔP fc .

Figure 13 .
Figure 13.Comparison of the monthly rate of deep injection and the monthly rate of earthquakes of M L 2.0+ for (a) Culberson-Mentone earthquake zone (CMEZ), and (b) the northern Delaware Basin (DB) deep injection region of New Mexico and TX.In CMEZ, a decoupling of the rate of injection and the rate of earthquakes has occurred in recent months, as injection rates have decreased dramatically but rates of seismicity have been slower to fall.In the northern DB region, injection rates above 50 MMbbl/month are associated with sharply increasing earthquake rates.

Table 1
Characteristics of Mapped Seismogenic Faults Selected for Stability Analyses, the Upper and Lower Bounds and Most Likely Stress State To Be Extant in the Vicinity of Each Fault as Defined by S , and Δσ vertical , the Native Pore Pressure, and the Modeled Pore Pressure as Sampled From Deep Injection Strata at Locations of Earliest Earthquakes of deepwater facies (chert and mudrock) to the southeast of this boundary during Silurian and Devonian times is evidence of drowning of the platform, likely with a tectonic control given the lack of evidence of global sea level rise Deposition Insight and can be accessed with a subscription; however, well data are also publicly available by individual well search through the Railroad Commission of Texas Oil and Gas viewer (https://www.rrc.texas.gov/resource-center/research/research-queries/)and New Mexico Energy, Minerals and Natural Resources Department Oil Conservation Division Online (https:// ocdimage.emnrd.nm.gov/imaging/LogFileCriteria.aspx).Software utilized for stratigraphic interpretation (IHS Petra), petrophysical modeling (Interactive Petrophysics), geomodeling (Schlumberger Petrel), and fluid flow modeling (CMG STARS) are available with license agreement.