Small Magnitude Events Highlight the Correlation Between Hydraulic Fracturing Injection Parameters, Geological Factors, and Earthquake Occurrence

Hydraulic fracturing (HF) operations are widely associated with induced seismicity in the Western Canadian Sedimentary Basin. This study correlates injection parameters of 12,903 HF stages in the Kiskatinaw area in northeast British Columbia with an enhanced catalog containing 40,046 earthquakes using a supervised machine learning approach. It identifies relevant combinations of geological and operational parameters related to individual HF stages in efforts to decipher fault activation mechanisms. Our results suggest that stages targeting specific geological units (here, the Lower Montney formation) are more likely to induce an earthquake. Additional parameters positively correlated with earthquake likelihood include target formation thickness, injection volume, and completion date. Furthermore, the COVID‐19 lockdown may have reduced the potential cumulative effect of HF operations. Our results demonstrate the value of machine learning approaches for implementation as guidance tools that help facilitate safe development of unconventional energy technologies.

Mitigation Area (referred to as Kiskatinaw, BCOGC, 2018), which requires mandatory publication of detailed injection parameters for individual HF stages. A comparison of earthquake density to active well locations using the dense and nearly full-azimuthal seismic station coverage in this area ( Figure 1a) shows that only a small fraction of wells or stages (∼15%) induce earthquakes above the detection threshold. Therefore, identifying the combination of relevant geological and HF parameters and quantifying their correlation with earthquake occurrence is key to mitigating hazard associated with induced earthquakes (Ghofrani & Atkinson, 2020).
Here, we use a supervised machine learning algorithm to correlate injection parameters from 12,903 HF stages with detectable seismicity in the Kiskatinaw area. We first enhance the existing earthquake catalog using an AI-based phase-arrival picker combined with template matching to lower the magnitude of completeness (M C ) and increase the number of detected earthquakes. We then add a binary target variable to the HF-stage data set based on whether a particular stage may be spatio-temporally associated with an earthquake from the enhanced catalog. Next, we use tree-based models (random forests and gradient boosted trees) combined with a tree-explainer to interpret the models. Our interpretation focuses on the importance of single injection parameters, parameter interdependence, and their distribution in different geological subformations. We show that the target formation thickness and distance to the Precambrian basement exert primary influence on seismicity occurrence. The majority of individual stage injection parameters exert a comparatively weaker influence, with the exception of fluid volume and treatment date. Our results highlight the benefits of machine learning approaches in determining the relative importance of the available set of HF-stage parameters for generating induced earthquakes. They also highlight the potential to decipher mechanisms driving induced fault activation where more detailed HF-stage parameters are available. This study also underlines the significance of open-data access to enable induced seismicity research aimed at mitigating seismic hazard and associated risks.
In the following, we first describe the catalog enhancement, followed by the machine-learning workflow. We then discuss the results for each stage parameter on a collective basis and the results for individual stages, while comparing our results to findings from previous studies.

Method
The machine learning approach determines the relative importance of both geological and operational parameters on a per-stage basis and their potential association with earthquakes in an enhanced catalog. Catalog enhancement uses a multi-station matched filter (MMF, Figure S1 in Supporting Information S1) on an initial earthquake catalog  combined with the AI-based phase arrival picker EQTransformer (EQT, Figure S2 in Supporting Information S1, Mousavi et al., 2020) from 12 July 2017 to 31 December 2020. We locate the detections and calculate magnitudes as done for the initial catalog and relocate earthquakes using hypoDD (Waldhauser, 2000). We then obtain detailed HF-stage injection information from the BCOGC database (BCOGC, 2022) for 865 HF wells within the spatial area shown in Figure 1 that were actively stimulated during the catalog period. We extract the start-/end-times and injected fluid volumes for individual stages, and associate geological parameters to each stage, such as the target formation and thickness. We follow BCOGC (2012) by separating the Montney Formation into Lower and Upper units. Davies et al. (2018) (2012). We then spatio-temporally associate earthquakes ≥M C in the enhanced catalog to individual stages (horizontal event-stage distance ≤5 km, 5-day window from the HF treatment, Figure 1b, Figure S4 in Supporting Information S1). The 5-day time-window limit is intended to account for cases of delayed triggering and could potentially lead to misassociation. However, as later discussed in the Results, 82% of events occur within 3 hr of a specific stage. The Supporting Information also provides further details of catalog enhancement, stage information retrieval, and earthquake-stage association.
Following earthquake-stage association, we apply the machine learning approach with workflow shown in Figure 1c. The approach sets the 13-parameter stage parameter data set as independent variables and the condition of whether or not the stage has an associated earthquake as the dependent variable. We use a stratified split to generate a training (80%) and test (20%) set, where the stratified split ensures training and test sets contain equal numbers of stages with and without earthquakes. We leverage a repeated k-Fold cross validator with five splits and 10 repetitions on the training set to obtain a validation set and estimate confidence intervals. We then follow Howard and Gugger (2020) using a random forest classifier to remove insignificant columns in the independent variable to reduce the likelihood of overfitting the data with excessive numbers of independent variables (Breiman, 2001, Figure 1c, Figure S5 in Supporting Information S1). We use the gradient boosted tree classifier XGBoost (Chen & Guestrin, 2016) to increase accuracy and support missing values present in the  (a1) shows the study area within the Montney Formation (BC and Alberta) in blue, as reported by BC Oil and Gas Commission. Black points denote earthquakes in the enhanced catalog, with colorbar indicating the number of earthquakes in each 2 × 2 km grid. Hydraulic fracturing (HF) wells with/without associated earthquakes (hexagons/squares) with respective well trajectories (gray lines) demonstrate a variable seismic response that is resolved by the dense seismic station coverage (blue triangles). Mapped faults (barbed lines) follow Norgard (1997); Berger et al. (2009);Cui et al. (2017) and Davies et al. (2018). (b) Schematic HF well (gray square -well head, gray lines -well trajectories, diamonds -HF stages) with associated seismicity as points color-coded according to their respective associated stages. (c) Supervised machine learning approach. Dark rectangles at the top denote the model's independent (input) and the dependent (output) variables. Left column shows stage parameters used as independent variables, with parameters of negligible importance shown in lighter gray (see main text and Supporting Information for details). Underlined parameters denote operational (in contrast to geological) parameters. The XGBoost (Chen & Guestrin, 2016) model is trained as described in the main text and the tree-explainer  calculates the SHapley Additive exPlanations values. well reports ( Figure S6 in Supporting Information S1). A test of other gradient boosted tree classifiers or neural networks demonstrated no significant increase in accuracy. We use the optimization framework Optuna (Akiba et al., 2019) to optimize the hyperparameters of XGBoost. We measure the performance of the model using the Area under the Receiver Operating Characteristic Curve (ROC-AUC) score, which is a good measure of separability (Fawcett, 2006). Text S4 in Supporting Information S1 provides extensive details regarding the machine learning workflow and ROC-AUC score.
We use SHapley Additive exPlanations (SHAP), Lundberg et al., 2020) to explain the model and its predictions. For each stage, we calculate a SHAP value as a measure of the contribution of each single parameter (input) to the output of the model, where the sum of the total resulting SHAP values from all parameters indicates the model prediction. A positive or negative SHAP sum corresponds to a stage being likely or unlikely to induce an earthquake, respectively. We also compare SHAP values of various HF parameters over all stages to interpret their relative importance as well as whether they are positive or negative (i.e., increase or decrease the likelihood of inducing an earthquake).

Results
Earthquake catalog enhancement adds an additional 21,404 events using MMF and 10,260 events using EQTransformer to the initial catalog of 8,382 events . The following discussion considers additional detected and relocated events with a maximum horizontal error of 10 km (median horizontal error: 0.4 km). Catalog enhancement reduces the magnitude completeness M C from ∼1.6 to ∼0.2 ( Figure S3 in Supporting Information S1). We infer the observed magnitude-frequency distribution likely results from a summation of multiple earthquake clusters related to individual HF wells with variable M C and b-values. However, a more detailed investigation of clusters would be required to draw any conclusions. We make no further interpretations of M C and the b-value here, as it is not the focus of this work.
A total of 39,156/40,046 (∼98%) enhanced catalog events associate with an individual stage, and 4,502/12,903 stages (∼35%) have at least one associated earthquake. Seismicity occurs predominantly in swarm-like clusters, as noted in Section 1 ( Figure S3b in Supporting Information S1), where 82% of associated events occur within 3 hr of the start a specific HF stage. Delayed earthquakes are comparatively few (16%), and 83% of all delayed events occur within 3 hr of the end of a stage ( Figure S4 in Supporting Information S1). Considering events exclusively within a specific treatment period reduces the number of stages with associated seismicity by ∼8% to 4,163 but does not change the results of the machine learning models. The reduction of the initial 13 parameters to five in the independent variable ( Figure 1c) using a random forest classifier (Section 2, Figure S5, Text S4 in Supporting Information S1) does not decrease the accuracy of the model. The fitted XGBoost model ROC-AUC score of 0.89 for both the validation and test sets suggests the model robustly generalizes the relation between stage parameters and earthquake occurrence with negligible overfitting ( Figure S6 in Supporting Information S1).
The SHAP-value distribution in Figure 2a for individual stage parameters shows all predictions calculated using all stages and the trained XGBoost model. Figure 2a1 shows the mean absolute SHAP values for all stages separated by their target formation. Both the individual and mean values demonstrate the importance of formation thickness, hydraulic-fracturing date, fluid volume, distance to the basement, and target formation. The following discussion describes the dependence of SHAP values with respect to the four most important stage parameters (formation thickness, date, fluid volume, and distance to basement). The HF-treatment date also correlates positively with the SHAP values for all formations up to the COVID-19 lockdown ( Figure 2c). The last completed stage in the database before the lockdown was on 2020-03-29 and the first stage completed after the lockdown occurred on 2020-08-06. All formations correlated negatively with SHAP values following the suspended operational period, particularly the Lower Montney. The SHAP values exhibit a more variable range following the lockdown, where fewer HF stages were completed comparatively to the period before. We further discuss the influence of the injection in Section 4.   Table S1 in Supporting Information S1 provides all SHAP values. (b) Formation thickness dependence with points colored by the same formation colors as in (a). Note: stages in the undifferentiated sediments have no calculated thickness. Colored bars above the points from a locally weighted scatterplot smoothing (LOWESS) for stages in the same formation (color as in (a)). LOWESS lines correspond to the y-axis on the right. The 2σ confidence intervals are calculated from 100 random samples. (c) Dependence plot for completion date with colors and symbols as in (b). The gray rectangle outlines the COVID-19 lockdown period. (d) Dependence plot for the total injected fluid volume with colors and symbols as in (b), but with distance to the Precambrian basement shown by the respective colorbar. (e) Dependence plot for stage distance to the Precambrian basement with colors and symbols as in (b). sediments show stronger dependence on the injected fluid volume compared to the Lower Montney. They are also negative for volumes ∼≤500 m 3 and positive for larger volumes, showing a slightly positive correlation with fluid volume. SHAP values for the Lower Montney are roughly zero with a slight positive correlation at higher volumes, and are positive for volumes ∼≥1,000 m 3 . Considering both the distance to the basement together with injection volume also highlights variable behavior (colorbar, Figure 2d). For example, the increase of SHAP values at 500 m 3 is more pronounced for stages located ≥2,000 m from the basement.
The stage distance to the basement negatively correlates with the SHAP values for distances ≥1,700 m, and SHAP values are primarily negative for distances ∼≥2,000 m ( Figure 2e). The SHAP values are highest at distances of ∼1,700 m in the Lower Montney, and decrease sharply with decreasing distance. Figure 3a summarizes the correlation of the SHAP-value sums with the maximum observed magnitude (MOM) for all stages, as well as the mean formation SHAP value as a function of magnitude. Panel (b) shows the composition of the SHAP sum for seven representative stages in (a), including four stages (1-4) with an associated earthquake, and three (5-7) without. Two (1-2) of the four stages with an associated earthquake (also the largest magnitude earthquakes in our study area) are deemed likely to occur according to the model ("true positives"), while two (3-4) are deemed unlikely to occur ("false negatives"). The three stages without an associated earthquake include one stage (5) where the model predicts no earthquake ("true negative") and two stages (6,7) where the model predicts an earthquake ("false positives"). SHAP sums show no clear correlation with magnitude, but the sum is positive for stages with a magnitude ≥0.5 for all formations. In addition, the locally-weighted scatterplot smoothing curves show the highest predictability for the Lower Montney (Figure 3a). In contrast, the number of false positives is also higher in the Lower Montney (see boxplots in Figure 3a).

Discussion
This study employs a machine-learning model to quantify the dependence of seismic propensity on detailed HF-stage parameters in the Kiskatinaw area using a significantly enhanced earthquake catalog from multi-year continuous monitoring. We find that geological parameters (formation thickness and vertical distance to basement) and operational parameters (HF treatment date and total injection volume) are key parameters for generating seismogenic HF stages. The remaining operational parameters tested appear to be of negligible importance. Furthermore, we find that stages targeting the Lower Montney formation are more likely to induce earthquakes than stages targeting shallower formations (Upper Montney, undifferentiated sediments). We emphasize that feature importance does not entail a causal relationship with predictive value. Instead, the model and subsequent interpretations are intended to provide a starting point for further analyses that put operational parameters in geological context. Our findings largely agree with similar studies focused on HF wells targeting the Montney Play. Wozniakowska and Eaton (2020) estimated the Seismogenic Activation Potential based on 6,466 HF wells (rather than specific stages) using logistical regression on published earthquake catalogs between 2006 and 2020 in Kiskatinaw and the North Peace Ground Motion Monitoring Area (NPGMMA), north-west from Kiskatinaw. They established injection depth and distance to the Cordilleran thrust belt as key parameters in their model. The importance of injection depth is also observed in Oklahoma (e.g., Hincks et al., 2018). The local focus on Kiskatinaw in this study implies that distances to the thrust belt are all similar for all HF stages considered, and are therefore not evaluated here. While injection depth is also essential here, we mapped the injection depth to specific digital geological models, as the geological units exhibit a significant NE-SW oriented dipping trend across Kiskatinaw ( Figure S7 in Supporting Information S1). Amini and Eberhardt (2021) employed a similar approach using tree-based machine learning models (as XGBoost, among others) with detailed information for 95,807 HF stages in the NPGMMA and Kiskatinaw with a combination of published catalogs containing 9,949 events, similar to Wozniakowska and Eaton (2020). They found that geological parameters are more important than operational parameters, similar to our results. In particular, both studies found a negative correlation with the distance to the basement, as demonstrated in Figure 2e, indicating that injection depths closer to the crystalline basement are more likely to induce earthquakes. Wang et al. (2022) follow a similar approach to this one (XGBoost model with SHAP values) while focusing on the Northern Montney Play (NMP), which includes the NPGMMA. They found that the distance to the Cordilleran thrust belt (similar to Wozniakowska and Eaton (2020)), and the total injected fluid volume (similar to this study and to ) are most important. They also found that treating pressures and the vertical distance to the Debolt formation are key parameters. However, our results suggest that The colorbar shows the point density for each subset of points. The lines colored by their respective formation result from a locally weighted scatterplot smoothing (LOWESS) over the MOM based on stages in each formation. The LOWESS lines correspond to the y-axis on the right with 2σ confidence intervals from 100 random samples. The colored boxplots on the left are calculated from stages with no observed earthquakes in each formation with the same colors used for the LOWESS lines. The notch of each boxplot represents the 2σ confidence interval from 10,000 random bootstrap samples. (b) The composition of the SHAP sum value for seven representative stages (marked by numbers 1-7 in (a)). Red and blue colors correspond to a positive and negative value for the specific parameter. Table S2 in Supporting Information S1 provides all parameter values. treating pressure bears no significant importance relative to other parameters, including the vertical distance to the Precambrian basement. The differences found here might highlight a variable seismic response to operations in the NMP relative to Kiskatinaw in the south. In contrast to the studies by Wozniakowska and Eaton (2020) and , we do not use the maximum horizontal compression orientation, SHmax, due to a lack of resolution in the small study region. (Published values in Kiskatinaw are based limited observations of ∼2 wells (Bell & Grasby, 2012;Heidbach et al., 2018;Shen et al., 2019)). A recent induced earthquake focal mechanism study in Kiskatinaw (Roth et al., 2022) found predominantly left-lateral strike-slip focal mechanisms with nodal planes at shallow angles to the regional SHmax, suggesting low spatial variability in SHmax in our study area. We note that SHmax is a parameter worth investigating in cases where observations exist at sufficient spatial resolution (e.g., from borehole measurements) near injection wells and/or seismicity. The abundance of strike-slip faulting mechanisms within the local study area also makes it difficult to use structural elements with potentially vertical throw (e.g., Wozniakowska et al., 2021) to identify correlations.
In contrast to previous models noted above, our model includes additional parameters. It suggests thickness of the target formation (for both Lower and Upper Montney) bears significant importance for the likelihood of inducing earthquakes (Figure 2b), particularly in the Lower Montney (Figure 2a1)). The overall correlation is positive, meaning a thicker target formation is more likely to host an induced earthquake; the correlation becomes negative when formation thickness increases to values ≥200 m. The positive correlation is also consistent with the theoretical model of Galis et al. (2017), which assumes that most induced earthquakes are self-arrested ruptures caused by local pore-pressure perturbations propagating on pre-stressed faults. One notable aspect is that the model takes the area of intersection between the reservoir and a pre-existing fault into account. The positive observed correlation with the formation thickness is consistent with an increase in the probability of having a potential pre-existing fault intersect a thicker reservoir or target formation while also increasing the size of an intersection. The more likely intersection with a fault should subsequently increase the likelihood of inducing an earthquake. We discuss below why the reverse trend in thicknesses ≥200 m does not fit the above interpretation, including what other factors could influence the observed negative correlation.
The HF treatment date also exhibits an important role. For example, stages completed in late 2019 to early 2020 are more likely to induce an earthquake compared with stages in 2018 (Figure 2c). The likelihood reverses after the seismically quiescent COVID-19 lockdown period (e.g., Salvage & Eaton, 2021). The positive correlation of increasing likelihood with time could point to a potential cumulative effect of injection, leading to increasing pore pressure, poroelastic, and/or Coulomb static stress changes (Brown & Ge, 2018). The lockdown likely relaxed stress perturbations to some degree, decreasing the probability of inducing earthquakes afterward. Oklahoma and Kansas also exhibit similar time-dependent effects where forecasted seismicity rates decrease with injection rates (Langenbruch et al., 2018). However, the time-dependent effects in the central US are related to wastewater-disposal-induced earthquakes attributed to pore-pressure perturbations over km length scales (e.g., Peterie et al., 2018). In contrast, perturbations for HF-induced earthquakes are likely more localized due to lower permeability, and may drive stress changes by both poroelastic and pore-pressure diffusion (e.g., Deng et al., 2016). For example, temporal poroelastic stress changes generate cumulative static Coulomb stress changes, which increase the likelihood of an induced earthquake by subsequent HF stages. We note that the interpretation related to any cumulative effect is speculative, albeit consistent with the seismicity rates following the lockdown. Confirming the interpretation would require testing other factors, such as the dependence on changes in stimulation strategies or other geological parameters. The lack of additional information in this data set prohibits further testing here. However, it does show the merit of testing the validity of cumulative effects on data sets using a longer time window where HF operations are sustained over longer periods ( Figure S3a in Supporting Information S1). The time-dependent observation highlights the potential benefit of the machine learning approach to identify essential injection parameters that should be studied in greater detail to identify fault activation processes.
The positive correlation with total injected volume and negative correlation with the basement distance has also been observed using similar approaches Wang et al., 2022;Wozniakowska & Eaton, 2020). Studies also correlate the volume to maximum magnitude (e.g., Galis et al., 2017;McGarr, 2014) and seismicity rate (e.g., Schultz et al., 2020). Figure 2d further highlights the detailed and variable response of the target formation, where correlation is positive for all formations. Where fluid volumes ∼≤500 m 3 appear to be less important for the shallower formations (Lower Montney, undifferentiated sediments), they increase 9 of 13 drastically for volumes ∼≥500 m 3 . The relationship could imply a minimum fluid volume is needed to reach deeper layers, presumably with more pre-existing faults. Once a minimum fluid volume is exceeded, it becomes important for estimating the likelihood of inducing an earthquake. Another explanation could be related to M C . For example, Schultz et al. (2018) show that a minimum fluid volume is necessary to increase probability of causing an event above M C . The observed increased likelihood associated with larger fluid volumes in shallower layers combined with larger basement distances (colorbar in Figures 2d and 2e) suggests the necessary combination of conditions to generate earthquake inducing stress perturbations in or near the basement. Similar dependences on the distance to the basement have also been observed for HF-induced earthquakes in the United States (Skoumal et al., 2018). Again, the combined importance of fluid volume and distance to the basement needs to be studied in more detail in future work, and goes beyond the scope of this paper.
One potential shortcoming of machine learning models is that the selection and availability of parameters may fail to represent the truth correctly or may introduce non-negligible bias. For example, the inverse correlation of the SHAP values with formation thicknesses ≥200 m 3 (Figure 2b) contradicts the theoretical model of Galis et al. (2017), which could suggest that the spatio-temporal distribution of stages biases the data. A possible explanation is that the observed trend of decreasing SHAP value with increasing formation thickness coincidentally reflects a higher number of HF stages targeting thicker formations completed after the COVID-19 lockdown. The temporal coincidence could result in a decreased likelihood of an induced earthquake after the lockdown period that overshadows the likelihood associated with thicker formations (Figure 2c). However, we do not observe a clear spatio-temporal distribution supporting the above hypothesis, as stage completions occur at arbitrary times in any spatial region (Figures S7,S8 in Supporting Information S1). In contrast, stages in the eastern part of the study area show significantly lower SHAP sum values compared to the center ( Figure S9 in Supporting Information S1). The lower eastern SHAP values also correspond to higher formation thickness of the Lower Montney and lower distances to the basement ( Figure S7 in Supporting Information S1). The spatial distribution of stages could explain the observed inverse correlation with formation thickness, and it is possible that other stage parameters that were not included in the model change in the same area. We infer from the distribution of mapped faults that the parameter changes might be related to changing geological factors, such as proximity to pre-existing faults ( Figure S9 in Supporting Information S1). Another potential bias could arise from the possibility that HF operations target thicker formations less frequently, resulting in less pronounced cumulative effects that can not be ruled out with the limited time window of our data set. Furthermore, our data set consists exclusively of stages that are a priori associated with seismicity. The prohibitive number of PDF files that would need to be processed manually to include non-seismogenic stages make the latter potential bias inevitable (see Section 2).
The unavailability of other potentially important parameters is another potential source of bias. As discussed above, the local SHmax orientation is likely influential, but lacks sufficient spatial accuracy in the study area. Another influential (unavailable) indicator could be lithium concentration in formation waters, which could be used as a proxy for vertical, fault-related fluid flow (Pawley et al., 2018). Formation overpressure has also been found to correlate with induced seismicity, but is also unavailable here (e.g., . Many potentially influential parameters could explore similar false positive and false negative cases shown as in Figure 3 to identify important stage parameters to be included in future models. In addition, false predictions offer the potential to identify methodological limitations in cases where more detailed information about the interaction between injection and seismicity is available. An additional potential bias might be unbalanced variable distributions, for example, the abundance of events, or lack thereof, within a formation unit. Since hypocentral resolution of our enhanced catalog does not permit association to formation units, we are unable to quantify the effect with this data set, but it demonstrates an avenue of future use. We note that we do not calculate seismic risk based on our machine learning models or attempt to predict the MOM, as done by previous studies. The fitting error using a regression model corresponds to ∼≥0.5 and includes significant outliers, as the larger magnitude range is not quantitively well constrained (where 8 stages events have M Lmax ≥3 and 544 stages with M Lmax ≥2). Therefore, the time frame considered will not provide a regression model for a robust risk assessment. The MOM in Kiskatinaw may be more strongly influenced by the local distribution of pre-existing faults or fluid conduits (Peña Castro et al., 2020), which such models are unable to take into account. Such limitations would also apply to using the total number of earthquakes or total moment release as the prediction target. In addition, the observed varying M C of different seismicity clusters would further bias the input data.

Conclusion
We investigate the relationship between an enhanced catalog of 40,046 earthquakes and detailed injection parameters for 12,903 HF stages using a supervised machine learning approach in the Kiskatinaw area, British Columbia, Canada. We use the gradient boosted decision tree classifier XGBoost (ROCAUC score 0.89) combined with the tree-explainer SHAP to assess the importance of specific injection parameters on the likelihood of an HF stage to induce a detectable earthquake. We find that: 1. Geological parameters show higher importance than operational parameters. 2. The thickness of the target formation, the completion date, the total injected fluid volume, the vertical distance to the Precambrian basement, and the target formation are the most influential parameters in order of relative importance. 3. Stages in the Lower Montney Formation have a higher likelihood of inducing an earthquake compared to the Upper Montney or shallower undifferentiated sediments. 4. Thicker target formations are generally more likely to host a seismogenic stage. 5. The increasing probability of HF stages over time may suggest the significance of cumulative effects of pore pressure, poroelastic, and/or Coulomb static stress changes in areas with HF. Cumulative effects seem to partially abate in Kiskatinaw during the COVID-19 lockdown in 2020. 6. Larger injected fluid volumes are more likely to induce an earthquake. Stages in shallower layers (Lower Montney, undifferentiated sediments) show a significant increase in importance for volumes ∼≥500 m 3 . The latter observation could point to a minimum fluid volume needed to reach deeper layers with more pre-existing faults. 7. Stages closer to the Precambrian basement are generally more likely to host a seismogenic stage. 8. The accuracy of our model is highest in the Lower Montney compared to the Upper Montney and undifferentiated sediments, whereas the false positive rate is also higher in the Lower Montney.
This study highlights a fundamental application of machine learning models toward induced seismicity: identifying essential parameters to be studied in greater detail in efforts to reduce the overall risk of industrial operations. Machine learning models also have a general potential to provide new information about the mechanisms driving induced seismicity. Although these models do not consider the geological context, they exhibit a potential to guide detailed studies focusing on the most relevant parameters. Furthermore, false-positive and false-negative predictions could be used in detailed studies of specific stages to provide further insights on parameters that might drive earthquake generation.
We emphasize that the open-access availability of HF stage information through BCOGC enabled this study. Open data availability is crucial to enabling scientific study to decipher the correlation between induced seismicity and HF (or other alternative energy production) operations, as the risk associated with induced seismicity is of great public concern. Publicly available data, as well as dense regional seismic networks to study potential relationships using small earthquakes will help, in ideal cases, to identify the processes leading to (and avoiding) damaging earthquakes. All the above steps will only help increase public acceptance and tolerance of the manageable hazard and risks associated with alternative energy production.

Data Availability Statement
Waveform data used in this study are archived at IRIS under network codes XL, 1E, PQ, and EO (e.g., https:// ds.iris.edu/gmap/XL). Well data are provided by British Columbia Oil and Gas Commission (BCOGC; https:// files.bcogc.ca/thinclient/, last accessed May 2022). Topographic information in Figure 1