Multivariate Statistical Appraisal of Regional Susceptibility to Induced Seismicity: Application to the Permian Basin, SW United States

Induced earthquake sequences are typically interpreted through causal triggering mechanisms. However, studies of causality rarely consider large regions and why some regions experiencing similar anthropogenic activities remain largely aseismic. Therefore, it can be difficult to forecast seismic hazard at a regional scale. In contrast, multivariate statistical methods allow us to find the combinations of factors that correlate best with seismicity, which can help form the basis of hypotheses that can be subsequently tested with physical models. While strong correlations do not necessarily equate to causality, such a statistical approach is particularly important for large regions with newly emergent seismicity comprising multiple distinct clusters and multifaceted industrial operations. Recent induced seismicity in the Permian Basin provides an excellent test‐bed for multivariate statistical analyses because the main causal industrial and geological factors driving earthquakes in the region remain highly debated. Here, we use logistic regression to retrospectively predict the spatial variation of seismicity across the western Permian Basin. We reproduce the broad distribution of seismicity using a combination of both industrial and geological factors. Our model shows that the proximity to neotectonic faults west of the Delaware Basin is the most important factor that contributes to induced seismicity. The second‐most important factor is salt‐water disposal at shallow depths, with hydraulic fracturing playing a less dominant role. The higher tectonic stressing, together with a poor correlation between seismicity and large‐volume deep salt‐water disposal wells, indicates a very different mechanism of induced seismicity compared to that in Oklahoma.

industrial activities, such as reservoir compaction (e.g., van Thienen-Visser & Breunese, 2015). Although these studies generally provide plausible mechanisms of physical triggering, due to the high level of model complexity required, they tend to focus on small regions featuring individual clusters of seismicity over subbasin scales. Such modeling-based studies also often fail to provide information on why some regions experiencing predicted stress increases remain aseismic (e.g., Keranen et al., 2014), possibly partly because the stress increment required to nucleate seismic slip in a given preexisting stress field is poorly understood.
Induced seismicity must also be modulated by natural geological controls (e.g., Keranen & Weingarten, 2018). Such geological factors may include the presence of faults that are sufficiently large (e.g., Rubinstein & Mahani, 2015) to produce earthquakes that can be felt or cause damage, and whether these faults are optimally oriented for failure (e.g., Sibson, 1990). The depth of injection/extraction may also contribute to seismicity due to variable pore-pressure diffusion and stress transfer in different lithologies, depth-dependent background stress frictional strength level (e.g., Kim, 2013), or formation overpressure effects (Ries et al., 2020). Thus, the development of effective methods to assess the susceptibility for induced seismicity, and hence seismic hazard, at a regional scale remains a key research challenge.
Statistical approaches provide a way to address this challenge. For HF, efforts have been made to generate statistical models that can highlight the factors which increase the probability of induced seismicity (Pawley et al., 2018;Ries et al., 2020;Wozniakowska & Eaton, 2020). In these cases, correlating well activity with induced seismicity is relatively straightforward because HF-induced seismicity is often tightly clustered around the injection well and period in space and time due to the short-lived nature of HF jobs (days to weeks), generating a highly localized pore-pressure diffusion field Schultz et al., 2020;Skoumal, Ries, et al., 2018). In contrast, other subsurface operations in the oil and gas industry can have a much wider footprint. For example, it has been shown that earthquakes induced by multiple high-volume and high-rate SWD over periods of months to years can occur over much larger space (e.g., >5 km) and time scales (e.g., earthquakes lagging operations by several months; Cochran et al., 2018;Goebel et al., 2017;Norbeck & Rubinstein, 2018;Rubinstein et al., 2018). Although efforts have been made to correlate seismicity with operations at individual SWD wells Teng & Baker, 2020), any association is less straightforward than for HF seismicity because of the SWD's larger footprint, combined effect of multiple wells, and larger time lag. Hincks et al. (2018) used a Bayesian Network to find that the proximity of SWD wells to basement was a strong predictor of high seismic moment release in Oklahoma. This study, however, only assigned any residual correlations to a single geospatial correction parameter, so was not able to fully explain the gelogical factors controlling the spatial distribution of induced seismicity in Oklahoma. Instead, full analysis of geological datasets could help to better constrain predictive models because induced earthquakes are more likely to occur along faults that are optimally oriented in the background regional stress field, and relative to any perturbed stress change due to anthropogenic operations.
Multivariate statistical studies, such as logistic regression (LR) models, have proven useful in predicting the occurrence of other geohazards, such as landslide susceptibility (e.g., Jessee et al., 2018;Nowicki et al., 2014). However, no previous studies have exploited this approach to understand induced seismicity and to assess the role played by the different triggering factors mentioned above, and how these vary across, and between large regions. One disadvantage of such statistical methods is that they cannot necessarily determine the physical mechanisms of induced and triggered earthquakes because they do not include physics, such as earthquake-to-earthquake triggering. Yet, statistical analyses help us to test competing hypotheses regarding the key controls on the location and magnitude of induced seismicity, and they may ultimately enable us to forecast the potential for induced seismicity to occur in given areas. To the west of the Mentone area, there is another intense area of generally low-magnitude seismicity (the "Culberson cluster") which seems to lie a long distance from the densest area of well activity (Figure 1a; Skoumal et al., 2020). In contrast, the northern portion of the Delaware Basin, which lies in south-eastern New Mexico, appears mostly aseismic, as does the Central Basin Platform. Further east, in the western part of the Midland Basin, there is a distinct cluster of earthquakes close to the cities of Odessa and Midland, which we call the "Odessa cluster." The Odessa cluster appeared to emerge in late 2018 to early 2019.
In Oklahoma, near-basement injection clearly induces most seismicity; however, it remains highly debated as to whether earthquakes in the western Permian Basin are primarily influenced by deep SWD into the Ellenburger Group (Lemons et al., 2019;Skoumal et al., 2020;Tung et al., 2021), HF in the Wolfcamp Group , shallow SWD into the Delaware Mountain Group (DMG; Deng et al., 2020;Zhai et al., 2021), or conventional production of hydrocarbons at shallow depth in the DMG (Deng et al., 2020;Doser et al., 1991Doser et al., , 1992. Moreover, in much of the region, hydrocarbons have been produced from the DMG since the 1970s, which may affect the preexisting state of stress (Dvory & Zoback, 2021). Hypocenter depths remain uncertain over the region due to a combined effect of uncertain velocity models station distribution, unclear direct arrivals due to complex subsurface structure, and high anthropogenic noise levels due to industrial activity Skoumal et al., 2020). Operational earthquake locations from the TexNet network suggest that most earthquakes occur at near-basement depths of 6-8 km , although there is a lot of scatter ( Figure S1 in Supporting Information S1). However, more recent analysis for earthquakes close to the boundary between Reeves and Pecos counties suggests shallower hypocenter depths of 1-3 km depth, consistent with the depths of the DMG (Dvory & Zoback, 2021;Sheng et al., 2020). When we consider multifaceted operations across a range of depths, together with hypocenter depth uncertainty, and the lack of detailed hydrogeological models for specific areas, we find that it is not easy to separate the different factors and their effect on seismicity. Another potential key difference between induced seismicity in Oklahoma and the western Permian Basin is the presence of numerous north-south striking mid-to-late Quaternary faults lying 30-50 km west of the edge of the Delaware Basin ( Figure 1a). These faults belong to the West Delaware Mountain Fault Zone, which strike NNW-SSE have a normal sense of offset, and slip at a rate of less than 0.2 mm/year (Collins et al., 1996). These recently active faults may indicate stronger neotectonic stresses compared with the Oklahoma region.
Although there have been several recent studies on seismicity in the Delaware Basin (e.g., Skoumal et al., 2020Tung et al., 2021;Zhai et al., 2021), they have tended to focus on the individual clusters described above rather than taking a holistic view of seismicity across the whole area. Our focus in this paper is to use a statistical approach which can help to identify broad correlations that help to hindcast the spatial evolution of seismicity over the large area. The spatial variability of seismicity across the western Permian Basin, and particularly the Culberson cluster which occurs far away from SWD wells, provides us with a unique opportunity to carry out a large-scale statistically driven analysis of what the combined geological and industrial controls are surrounding induced seismicity in this region.
In this study, we use LR analysis to model the spatial distribution of earthquake occurrence between 2017 and 2020 across the western Permian Basin. We uniquely combine geological factors with industrial variables based on detailed data sets of fluid injection and extraction from hundreds of thousands of wells in the region. While correlation does not necessarily mean causation, and because our approach takes a time-integrated view and cannot model the effects of all factors that trigger earthquakes (e.g., earthquake-to-earthquake triggering), our resulting LR model demonstrates that the main zones of seismicity in the Permian Basin can be robustly modeled using a small number of both industrial and geological features.

Target Variable: Earthquake Occurrence
Our earthquake database comes from the TexNet catalog . The TexNet catalog starts in January 2017. We do not include pre-2017 seismicity from the USGS ComCat catalog due to inherent epicentral uncertainties prior to installation of regional seismic monitoring networks. Since our area of interest (hereafter, AoI) also covers south-eastern New Mexico, we supplement the TexNet data set with additional events reported by the New Mexico Tech Seismological Observatory (Pankow et al., 2019). Although the magnitude of completeness (M c ) in certain areas across the AoI reaches as low as M L 1.2 (Savvaidis et al., 2019), given that our AoI covers areas away from the main seismicity clusters and that we merge two catalogs, we opt for a more conservative M c of 2.2 based on these considerations, together with formal assessment of departures from a regional Gutenberg-Richter distribution. To start with, we generate a target variable based on seismicity, which represents the spatial distribution of earthquake occurrence. We divide up the AoI, use the NAD83/Texas State Mapping System (EPSG:3081) projected coordinate system, with a uniform spacing in the x and y directions of 5 km (Figure 1b). We choose this grid spacing because it remains larger than the typical one-sigma epicentral errors and it allows us to define a minimum search radius for industrial activity of 5 km (see Section 2.2).
We assign each grid node a value of one if at least one M L ≥ 2.2 earthquake has occurred within 5 km of the central point of each grid square since 2017. Otherwise, if no earthquakes locate within that grid node, we assign a value of zero. Using this distance kernel ensures a smoother distribution of seismicity that does not suffer from artifacts caused by epicenters lying close to the edge of a grid square. The resulting binary 2-D grid of the earthquake occurrence target feature is shown in Figure 1b and is similar to that used by Ries et al. (2020) used to study HF-induced seismicity in Oklahoma. One advantage of using this binary target variable is that it avoids the need to decluster the seismicity catalog. In addition, our approach negates earthquake depths, which remain highly uncertain in the region (e.g., Savvaidis et al., 2019), and do not correlate directly with depths of fluid injection and extraction ( Figure S1 in Supporting Information S1). In our target grid, 16% of the 2,655 grid points are assigned the value "1."

Candidate Predictive Features
Our model features include both industrial data provided by the IHS Markit database (see Data Availability Statement) and a range of geological data, as described below.

Geological Features
We use a set of geological features that includes depth-to-basement (Figure 2a), because basement depth has been found to be an important factor seismicity induced by SWD and HF (Hincks et al., 2018;. We also use information about preexisting faults and their stress state. In the western Permian Basin, there is a normal faulting stress regime, as indicated by an Anderson fault shape parameter (Simpson, 1997), A Φ = 0.5-0.9 (Lund Snee & Zoback, 2016, 2020; Figure 2b). Further to the east, the stress regime becomes a mixture of strike-slip and normal faulting. Overall, the predominance of normal faulting means we can use orientations of the maximum horizontal compressive stress (S Hmax ) from the same studies as above (Figure 2c), to derive a simple proxy of how optimally oriented faults are to failure. Such a simplified approach is needed because we do not have adequate constraints on fault dip angles to derive a more quantitative slip tendency value. We use a variety of fault data sets to derive smoothly varying fields of fault strike that allow us to compute the cosine of the angular difference with S Hmax . S Hmax orientations align quite strongly with lineations of seismicity in the Delaware Basin (Figure 2c), supporting the use of this simple angular difference metric.
For subsurface faults within the Permian Basin, we use the basement-rooted fault database of Horne et al. (2021). For areas outside the study area of Horne et al. (2021), we include faults in the Woodford formation from Comer (1991) (see Figure 2d). We also include a database of Precambrian basement structures from Ewing et al. (1983) and Comer (1991) (see Figure 2e). These maps of subsurface faults are largely derived by well correlations and stratigraphic thickness changes rather than direct imaging. For each grid point, we compute both a minimum distance to nearby faults and a smoothed mean fault strike and S Hmax azimuth. We then compute the cosine of the resulting angular difference between S Hmax and fault strike with lower values indicating a more optimal fault orientation in the regional normal faulting stress regime. While such fault maps are likely to be spatially incomplete, our smoothing approach captures the broad variation in fault orientations across the AoI governed by control points where fault data are available. Finally, given the proximity of the active Rio Grande rift system, ∼150 km to the west of the western boundary of the Delaware Basin, we use the USGS Quaternary Faults database ( Figure 2f; Collins et al., 1996) to compute the distance between each grid point and the closest recently active fault, as a simplified proxy for preexisting tectonic stress.

Industrial Features
We use well metadata and corresponding injection and extraction time series data for wells in the study area of the western Permian Basin from the commercial Production-Allocated and Well database provided by IHS Markit (see Data Availability Statement). The well metadata used comprise well locations, total vertical depths, and formation top at total depth information. The time series data comprise fluid (oil, gas, water) extraction and SWD injection volumes given to a month-level resolution. These data are originally based on operational data reported by the regulator, the Railroad Commission of Texas (RRC). The locations of injection and extraction wells are shown in Figure 3. We supplement this data set with HF data reported by the FracFocus database (Dundon et al., 2015;Figure 3b). We include data for all wells for the geographic AoI covering ∼70,000 km 2 (Figure 1), with monthly injection and extraction data covering the period January 2017 to July 2020 to ensure sufficient regional geographic data completeness across all wells, and so that this timeframe temporally coincides with the seismicity data set. We tested whether the inclusion of operational data from prior to 2017 affected results; however, these data are highly correlated with the post-2017 data, so were not included as a separate feature. This spatial correlation over different time periods means that recent industrial activity is representative of operations over a longer period, even if contemporary seismicity may be driven by a build-up of industrial operations over time. The approach of using post-2017 data broadly represents recent fluid injection and extraction rates, the former of which has been shown to correlate with seismicity in Oklahoma (e.g., Keranen et al., 2014).
For fluid injection and extraction volumes, we compute total cumulative monthly volumes within a given circular radius from the center point of each grid point. We consider three possible search radii of 5 km (i.e., the grid spacing), 10, and 25 km. These length scales were chosen to be representative of distances relevant for nearfield direct pore-pressure effects and longer-distance poroelastic triggering due to different industrial operations. For HF, we use a single distance of 5 km, consistent with pore-pressure effects from short-lived injection (Ries et al., 2020;Schultz et al., 2020). For SWD, we also use the larger distances of 10 and 25 km to encompass the possibility of coupled pore pressure and longer-range poroelastic effects (e.g., Goebel et al., 2017). Based on the interpreted well formation tops at total depth from IHS's PRODFIT database, we divide up the injection and extraction volumes into three different groups: "Shallow," "Wolfcamp," and "Deep." We do this because of the variable injection and extraction volumes at different depths in the Permian Basin ( Figure S1 in Supporting Information S1). The Shallow group includes the upper Permian and Triassic formations, such as the DMG and Bone Spring Group in the Delaware Basin. Wolfcamp refers to the Lower Permian Wolfcamp Shales, which are hydraulically fractured for hydrocarbons. Deep corresponds to deeper, sub-Permian formations, such as Carboniferous, Devonian, and Silurian formations, such as the Ellenburger Group. For SWD, we include a further feature of total injection occurring at greater than 2,000 m depth and where the injection depth is within 1,000 m of the top-basement, since that is a factor believed to strongly drive seismicity in Oklahoma (Hincks et al., 2018). The locations of wells which inject into deep and shallow formations are shown in Figures 3c and 3d, respectively. For HF, we compute both the total number of frac jobs within 5 km each grid point as well as the total volume injected.

LR Workflow
Binary variables, like our target feature of earthquake occurrence, can be fit using models of LR (e.g., Cox, 1958;James et al., 2013). LR is a simple yet powerful statistical function using a linear model to predict the log-odds ratio of a binary outcome ("target variable"). The linear model is expressed as linear combination of multivariate input data ("features"). A LR approach is adopted because we need to consider a large number of factors that have already been linked to induced seismicity in the Permian Basin (e.g., Dvory & Zoback, 2021;Savvaidis et al., 2020;Skoumal et al., 2020;Tung et al., 2021;Zhai et al., 2021). Also, it can be seen from Figures 3 and 4 that there is no strong relationship between individual industrial operations or geological features and seismicity, showing that a multivariate approach is required.
We standardize and normalize all input features using the Yeo-Johnson power transform (Yeo & Johnson, 2000;hereafter, YJ) to reduce data skewness and to make it more Gaussian-like, which is important for variables such as injection volumes which cover many orders of magnitude, as well as zero values. This feature normalization thus allows us to interpret relative differences between feature coefficients. Examples of raw and transformed features are shown in Figure 4.
As some industrial activities strongly overlap in space, they cannot be effectively distinguished due to high collinearity. We assessed feature similarity using a clustered matrix based on Pearson correlation coefficients computed on the transformed features ( Figure 5). Where features are strongly colinear, we removed, grouped, and renamed certain features to account for this. One key example is SWD and production from shallow layers, which share a correlation coefficient of 0.9 ( Figure 5). We therefore group these features, use a single feature, and rename them to "shallow injection/extraction." We generate a LR model with a forward stepwise approach, in which we iteratively add in feature variables while ensuring that the p-value of each new feature remains below 0.05 (Ries et al., 2020;Teng & Baker, 2020) to ensure that each feature is statistically significant at the 95% confidence level (e.g., Nowicki et al., 2014). This process eliminates insignificant features but still does not account for any remaining high collinearity between variables. For example, in the case where the spatial distribution of fluid injected into the Wolfcamp Shales for HF strongly correlates with the location of oil produced from the same formation, so we compute the Variance Inflation Factor (VIF; Jessee et al., 2018;Midi et al., 2010) for each feature and iteratively remove the most colinear variable over a given number of steps. The number of steps we choose to remove features for is based on manual  Figure S2 in Supporting Information S1. assessment of several fit-quality parameters. These are described as follows: (a) the F1 score which accounts for model precision and recall, especially due to the large class imbalance (more true negatives than true positives) for our target feature (e.g., Chicco & Jurman, 2020). (b) The AUROC (area under the receiver operating characteristic) score which is a metric for the model's ability to discriminate between classes (e.g., Fawcett, 2006). (c) The model's pseudo-r 2 score ( 2 ; McFadden, 1973), which measures the amount of variance explained by the LR model. (d) Akaike's Information Criterion (AIC; Akaike, 1974). (e) The Bayesian Information Criterion (BIC; Schwarz, 1978). AIC and BIC assess model fit by penalizing the inclusion of additional variables. Smaller values of AIC and BIC indicate a better-performing regression model based on log likelihood and complexity. (f) Moran's I-number for the model residuals (MI r ; Moran, 1950), which measures the extent to which the model residuals are spatially correlated. We also ensure that the maximum feature VIF does not exceed 5, a commonly used cutoff in regression studies (e.g., Stine, 1995). Visual inspection of the feature correlation matrix ( Figure 5) verifies that this approach successfully eliminates highly colinear variables with a correlation coefficient of >0.8.

Results
In this section, we consider the performance of different model classes in terms of the values of various fit-quality parameters described in the preceding section (F1 score, AUROC score, 2 , AIC, BIC, and MI r ). These parameter results from the different models considered in this section are shown in Table 1.

Null Hypothesis: Univariate Regression Using SWD
As a first step to explore which operational parameters best fit the earthquake occurrence target variable, we first test the null hypothesis that deep SWD is a primary influence on seismicity, as has been inferred for Oklahoma (e.g., Goebel et al., 2017;Hincks et al., 2018), and has been proposed for the Delaware Basin Skoumal et al., 2020Tung et al., 2021). We therefore perform univariate regression using deep SWD feature only to test this hypothesis.
Deep SWD is a very poor predictor of the spatial distribution of earthquake occurrence, as it is unable to model any of the observed locations of earthquake occurrence (F1 score = 0.0), and the 2 is very small (0.02; Table 1). This result arises because in the study area the highest SWD volumes mostly lie to the north of the main seismicity clusters, as seen through both the distribution of individual high-volume disposal wells (Figure 3c) and the summed contribution from multiple wells (Figure 4e). This finding suggests that we can confidently reject the null hypothesis that, on its own, deep SWD volume predicts the spatial distribution of seismicity throughout the western Permian Basin. Even the possibility of a near-field aseismic zone and a longer-wavelength triggering hypothesis, associated with deep SWD wells (Goebel et al., 2017;Guglielmi et al., 2015), does not seem to be supported by the location of high-volume disposal wells relative to the main clusters of seismicity.
Other univariate regressions were also performed (Figure 4), and no individual feature demonstrated strong predictive performance on its own.

Optimum Multivariate Regression Model
Given the clear rejection of the null hypothesis that deep SWD correlates with seismicity, we compute the best fitting LR model using a combination of geological and industrial features for a more thorough exploration of the solution space. The results of this model are shown in Figure 6.
We test the data-dependent robustness of this model through a bootstrap approach in which we randomly remove 20% of the input data grid points and repeat this process 2,000 times to generate uncertainties in the resulting model coefficients. This test allows us to assess the sensitivity of the model to small variations in the distribution of earthquake occurrence. These uncertainties are shown as histograms in Figure 6b and show that all significant features in the optimum model are robust, with very little contribution from the features that do not feature in the final, optimum model. Although we do not have a sufficiently large seismicity data set to generate a training, testing, and validation data subsets, this Bootstrap resampling shows that we are not overfitting the data. We also compute the standard errors in the output model using the delta method, which uses a Taylor series expansion of the inverse logit function of the regression. The upper and lower confidence interval in the binary model space is shown in Figure 6c.
Our model reproduces the main cluster of seismicity in the Delaware Basin, and within upper and lower confidence intervals (Figures 6a and 6c) Note. The preferred and best-performing model is italicized. AUROC, area under the receiver operating characteristic.

Table 1 Model Performance of the Four Model End-Members Tested
the main clusters of seismicity. The normalized coefficients of the model features are shown in Figure 6b. The model's of 0.37 indicates that the model has fair predictive power (approximately equivalent to a traditional r 2 value of ∼0.65 [Domencich & McFadden, 1975]). We discuss the relative importance of each of these features in Section 4. By testing a smaller AoI covering the Pecos and Reeves seismicity clusters in the Delaware Basin, we find that these relative feature coefficients are stable and insensitive to the exact definition of the AoI (Figure S3 in Supporting Information S1). The distance to quaternary faults is the most important feature of the LR model. With Quaternary Faults lying 30-40 km from the western edge of the Delaware Basin (Figures 1a and 2f), indicating a higher neotectonic stressing rate than in Oklahoma, we thought it plausible that the distance to the closest quaternary fault feature would help to reproduce the seismicity at the western edge of the Culberson cluster. However, our LR model shows that we are still unable to reproduce the westernmost edge of the Culberson cluster, highlighting that we are perhaps missing features, such as inter-earthquake triggering, which may also be specific to certain regions. Alternatively, a fault zone with very high diffusivity could allow long-range earthquake triggering in such a way that is not captured by our approach of assuming a spherically constant search radius. The second-most important industrial feature is SWD into shallow formations. We are unable to recover the Odessa cluster of seismicity which indicates that the causative factors of seismicity across the Permian Basin vary spatially. In other words, the factors causing seismicity in the Midland Basin are different to those in the Delaware Basin. We are also unable to recover small clusters of seismicity in the northern parts of the Delaware Basin and Central Basin Platform, especially in the New Mexico-Texas border region.

Need for a Hybrid Industrial and Geological Model
To demonstrate that a mixture of industrial and geological factors is needed to explain the spatial distribution of seismicity in the western Permian Basin, we consider two end-member models in which we consider industrial features only and geological features only. The results from these two models are shown in Figure 7. The prevalence of industrial features such as shallow injection/extraction and Wolfcamp HF remains as per our optimal hybrid model. These industrial features on their own manage to describe the spatial distribution of the main Pecos cluster of seismicity. Industrial features are a stronger overall predictor of seismicity, which is in line with the general hypothesis that earthquake occurrence in the Permian Basin has a significant induced component. However, the proximity to Quaternary faults feature is not prevalent in the geological-only model, implying that this feature works in tandem with the industrial features to reproduce the spatial pattern of seismicity. Although the industrial-only model performs better than the geological-only model, both perform significantly less well than the hybrid model, as shown by the corresponding 2 , F1 score, AUROC, AIC, BIC, and MI r values (Table 1). Nevertheless, these different end-member models help us to understand how different combinations of features help to explain predicted seismicity distributions.
The Pecos and Odessa clusters correlate quite strongly with HF activity (Figure 4a) and shallow injection/extraction (Figure 4d). The lack of HF on the Central Basin Platform likely drives the lack of predicted seismicity here. Deep SWD appears to be quite prevalent in the Odessa region, but this is not a large enough part of our model space to become a significant feature. The basement-rooted faults also appear less well optimally oriented in the Central Basin Platform (Figure 4c). Depth-to-basement is prevalent in the geological-only model; however, this overpredicts seismicity in the south-eastern part of the Delaware Basin where the basin remains deep (Figure 2a), which is why this feature is not significant in our hybrid model.

Discussion of the Hybrid LR Model
We find that a combination of industrial and geological features is needed to accurately predict the distribution of seismicity in the western Permian Basin. Although our approach does not allow us to derive definite physical mechanisms of induced seismicity, it permits us to perform hypothesis testing with each of the statistically significant features in our preferred model, considering known triggering mechanisms both in the Permian Basin specifically, and more generally, from the published literature. We discuss each of the features below in order of their significance in the LR model.

Proximity to Shallow and Recently Active Faults
The feature with the largest coefficient in our LR model is a geological parameter: the proximity to shallow and recently active (i.e., Quaternary) faults. This feature particularly helps to reproduce the seismicity at the southwestern edge of the Delaware Basin. This main difference can be seen clearly by comparing the industrial-only and geological-only LR models (Figure 7).
The Delaware Basin is bounded to the west by a series of NNW trending normal faults belonging to the West Delaware Mountain Fault Zone, trending subparallel to the Rio Grande rift zone (Collins et al., 1996;Muehlberger et al., 1978). This area has hosted moderate-to-large earthquakes before, such as the M w 6.3 Valentine earthquake in 1931 (Doser, 1987;Dumas et al., 1980;Storchak et al., 2013), and a M w 5.7 earthquake in 1995.
Given that active deformation occurs along the Rio Grande rift region (Murray et al., 2019), we speculate that the stress associated with this deformation may extend to the Quaternary faults in West Texas. Overall, our results suggest that one key difference between induced seismicity in Oklahoma and West Texas is the latter's proximity to recently active faults, and hence might have stronger preexisting tectonic stress that may modulate induced earthquakes. The role of these recently active faults has not yet been considered in studies of induced seismicity for the Delaware Basin.

Shallow Injection/Extraction
Injection/extraction (search radius = 10-25 km) in the shallow DMG and Bone Spring Group is the second-most important feature of our LR model. Yet, there is currently little evidence that seismicity is occurring in such shallow layers, although the hypocentral depth distribution of seismicity in the western Permian Basin remains debated. It has been recently suggested by Zhai et al. (2021) that shallow wastewater injection may increase stress on near-basement faults through poroelastic stress transfer. Our LR result provides some additional evidence that such shallow industrial activities may cause seismicity in the Delaware Basin. Moreover, Deng et al. (2020) and Staniewicz et al. (2020) proposed that by modeling InSAR observations, shallow extraction produces long-term surface subsidence, which aligns with seismicity in the Pecos area, and shallow injection produces signals of uplift.
High colinearity between pre-2013 shallow oil production and post-2017 shallow water injection (r = 0.88) means that we cannot confirm the hypothesis of Dvory and Zoback (2021) that shallow injection into the same formations that were produced from earlier production inhibits seismicity.

HF in the Wolfcamp Shale
The third most dominant feature in our LR model of earthquake occurrence is the number of HF jobs carried out in the Wolfcamp Shales. Even in the model derived from industrial operational features only (Figure 7), and in a univariate sense (Figure 4), this feature consistently has a good correlation with earthquake occurrence. Since this feature is highly correlated (r > 0.85) with both the oil and water production from the Wolfcamp Shale ( Figure 5), we are unable to fully distinguish between numbers of HF jobs, HF volumes, or the resulting production amounts, although in terms of statistical significance, the number of HF jobs is marginally favored. HF helps to fit clusters of seismicity in the Delaware Basin and in the Odessa region (Figure 4a).
HF can trigger seismicity due to localized pore-pressure increases on nearby faults (Schultz et al., 2020). In many documented cases, HF-induced seismicity occurs close (<5 km) to the injection well and within the same depth range as the shale target formations. However, HF-induced seismicity has also been observed in the basement lying several kilometers beneath the shale target formations (Lei et al., 2019). Using a probabilistic distance-time likelihood association of seismicity with well activity, HF was also shown by Savvaidis et al. (2020) to be one of the main factors causing seismicity in the Delaware Basin. Moreover, Dvory and Zoback (2021) showed that pore-pressure perturbations from HF could feasibly trigger shallower seismicity in the DMG and Bone Spring Groups. Earthquake depths remain highly uncertain in the Permian Basin, so it is difficult to ascertain whether some earthquakes occur within, above, or below, the Wolfcamp Shales. Even if the seismicity might not be directly related to HF, it remains an important question as to whether production from the Wolfcamp Shale may affect the state of the stress in the deep formations beneath or the shallow formations above. Nevertheless, given the significance of HF in our LR model, we speculate that HF activity affects the stress state in either the Wolfcamp Shale or the surrounding formations that encourages seismicity.

Deep SWD
Another feature of the model is the statistical significance of SWD volumes into deep formations, at a radius of 25 km. As opposed to the other features discussed above, this one has a negative coefficient, implying that SWD might impede seismicity within 25 km, although this distance has a high uncertainty since similar results can be obtained using the 10 km radius. As can be seen in Figures 3 and 4, the largest volume SWD wells are located far to the north of the main Pecos cluster of seismicity, in southern New Mexico. Therefore, there is no straightforward spatial correlation between SWD and seismicity. The negative contribution of SWD most likely helps to reproduce the absence of seismicity on the Central Basin Platform. Since HF activity extends further west than the regions of high-volume SWD, the combined factors help to predict a concentration of seismicity in the Pecos and Mentone regions. Therefore, deep SWD is a spatially clustered feature that helps to replicate the observed pattern of seismicity when combined with other features. Consequently, we cannot say that SWD physically inhibits seismicity across the western Permian Basin. Instead, if anything, it might relate to a unique state of stress in the Central Basin Platform area. Alternatively, this negative feature may indicate that our model is lacking another feature to explain this seismicity.
Recent studies have made the link between the 2020 M w 4.8 Mentone earthquake and deep SWD Tung et al., 2021). However, there are few studies (Skoumal & Trugman, 2021) that link deep SWD to widespread seismicity in the western Permian Basin. Our result of a negative, albeit relatively small, coefficient for deep SWD might imply an aseismic region surrounding high-volume wells, which might be caused by aseismic slip along faults (Guglielmi et al., 2015), a dominance of long-distance poroelastic triggering over nearfield pore-pressure effects (Goebel et al., 2017), or the distance between large SWD injectors and seismogenic faults.

Optimal Orientation of Basement-Rooted Faults
It is reasonably clear that many of the earthquakes that occur in the Delaware Basin do not occur on faulting structures mapped within the basin (Figure 2). For example, in the Pecos region, most mapped basement-rooted faults in the basin strike approximately WNW-ESE, but the seismicity occurs along NW-SE lineations. However, given that the theoretical optimal orientation of these faults (assuming an extensional stress regime, and based on estimates of S Hmax ) appears as a statistically significant feature, our results show that the orientations of these faults have mild predictive power and they increase the propensity of induced seismicity. Although the mapped faults do not exactly align with the lineations of seismicity (Figure 2), the alignment is close enough to drive the observed correlation in our LR model. It is likely that the faults that host many of the earthquakes in the Permian Basin, which are typically magnitude 4 and less, might be subseismic in scale, and hence do not cause the large changes in stratigraphic thickness typically required to be recognized in subsurface data sets.

Summary and Limitations of Our LR Model
Overall, our LR results show that when combined, geological and industrial factors produce a robust correlation with the spatial distribution of seismicity. Although correlation does not equal causation, this method has allowed us to test some hypotheses for causal factors of induced seismicity which we discussed and can be tested in the future with physical models that factor in these different proposed mechanisms. Our approach is highly adaptable to different regions and to different data sets. For example, our method could be expanded to include geodetic maps of deformation. Based on our trained model, we could then expand that AoI to test the stability of model features for a wider region of seismicity in the Permian Basin. Moreover, our method is straightforward to update in the event of a newly emergent cluster of seismicity, updated industrial operational data, or new geological data. Eventually, our approach may help guide regional probabilistic maps of seismogenic potential.
Our approach currently considers the spatial, time-integrated distribution of seismicity and industrial data. Including time-varying features and targets would vastly increase the complexity of such a regression model but could improve the model fits to certain clusters of seismicity. For example, interearthquake triggering (e.g., Coulomb stress transfer) may account for seismicity in the Culberson area, which is not predicted by our model.

Conclusions
Using a spatial LR method, we have provided new insights into multivariate statistical correlations with induced seismicity over a large portion of the Permian Basin, Texas, which includes multiple subbasins. Our work is the first study that we are aware of that considers multiple seismicity clusters across the Permian Basin in a single model. Based on thorough analysis of industrial and geological data, our results demonstrate that a combination of industrial and geological factors is required to hindcast seismicity in the Permian Basin.
We find that distance to the closest Quaternary fault is a key predictor of the spatial distribution of seismicity. This result indicates that, in contrast to induced seismicity in Oklahoma, a higher rate of background tectonic stress is important in determining the spatial pattern of triggered seismicity in the Permian Basin. This higher rate of background tectonic stress should therefore be accounted for when assessing seismicity rates and hazard due to anthropogenic activities in West Texas. Because earthquake depths in the region remain uncertain, we cannot determine whether direct stress changes from shallow SWD or the indirect changes due to extraction, which could extend to the deep basin and basement, are the main cause of seismicity. Also, in contrast to Oklahoma, it appears that deep SWD is not a dominant factor in determining where induced seismicity occurs. Deep SWD plays a weakly negative role in our model and it likely relates to the lower rate of seismicity on the Central Basin Platform, rather than taking a primary role in reducing seismicity overall.
Our modeling approach could be applied to other regions or our model adapted to broader regions, for example, regions where high-volume SWDs exist, but induced seismicity does not occur (Rubinstein & Mahani, 2015), to identify the reasons for the apparent aseismicity. Such regions include the Williston Basin in North Dakota (Frohlich et al., 2015) and along the Gulf Coast (Weingarten et al., 2015).