Source Properties of Hydraulic‐Fracturing‐Induced Earthquakes in the Kiskatinaw Area, British Columbia, Canada

This work presents a high resolution source property study of hydraulic fracturing induced earthquakes in the Montney Formation, a low‐permeability tight shale reservoir in the Kiskatinaw area, northeast British Columbia, Canada. We estimate source parameters, including focal mechanism solutions (FMSs), seismic moment, spectral corner frequency, and static stress drop values of earthquakes recorded between July 2017 and July 2020. Waveform‐similarity‐based event classification of 8,283 earthquakes yields 52 event families (clusters) and 1,014 isolated events (individuals). We calculate a total of 64 FMSs of events ML > 2.5 with high‐quality waveforms. Of the 64 solutions, 54 come from events within families, and are used to infer an additional 3,500 focal mechanisms of smaller‐magnitude events with similar waveforms. The other 10 are isolated events. The dominant faulting style inferred from FMSs highlights multiple cascading, shallow, strike‐slip events and generally isolated, larger‐magnitude reverse‐style events in close proximity to the Fort St. John Graben system. Inferred nodal planes of strike‐slip events are at low‐angles to regional SHmax, suggesting optimally oriented, left‐lateral faults. Reverse faulting nodal planes are roughly perpendicular to SHmax and have an orientation consistent with the reactivation of pre‐existing normal basement faults. Source spectral analysis using three approaches, including spectral‐ratio fitting, suggests a constant stress drop of 1–10 MPa and self‐similarity for induced events. Constant stress drop scaling breaks down at magnitudes smaller than ∼ML 2.0, likely due to observational bandwidth limitations.

Source parameter estimates are an important aspect in determining differences in the rupture process of induced and tectonic earthquakes. Some studies estimate that stress drop values of induced and tectonic earthquakes in the central and eastern United States are on the same order of magnitude (Ellsworth, 2013;Huang et al., 2017). Boyd et al. (2017) observe differences correlated with varying hypocentral depth and faulting style. A number of studies in the WCSB observe similar stress drop values between induced and tectonic events over a broad region in Alberta and northeast British Columbia (Clerc et al., 2016;Zhang et al., 2016). However, local-scale studies reveal differences that are consistent with resulting from heterogeneous rock mechanical properties. For example, Yu et al. (2019Yu et al. ( , 2020 interpret lower stress drop values within ∼1 km of HF wells to result from rupture in the damaged rock volume near the wellbore. Distinct differences in rupture behavior may be resolvable from data sets with dense azimuthal coverage. But, comparing source parameters of small earthquakes from different studies is challenging, due to differences in instrumental coverage and varying methodological approaches. For example, Holmgren et al. (2019) emphasize the importance of the underlying models and constants in estimating corner frequencies and stress drop values of 87 induced earthquakes in the WCSB. Their results show variations with station azimuth by more than a factor of four in some cases.
Nevertheless, some of the variability between seismic source parameter estimates, particularly within individual studies, may be attributable to differences in source behavior. The Montney Formation, the largest unconventional shale gas play within the WCSB, with over 6,400 production wells (Wozniakowska & Eaton, 2020), has hosted several events of magnitudes 4+, including an M w 4.6 event on 17 August 2015 near Fort St. John (Babaie Mahani et al., 2017;, and an M w 4.2 (M L 4.5) event on 30 November 2018 near Dawson Creek (Babaie Mahani et al., 2019;Peña Castro et al., 2020). Both events likely resulted from reverse-faulting in basement rocks, whereas the relocated aftershock sequence of the 2018 event in the shallower sedimentary layers provides evidence for strike-slip faulting as well. The above studies are limited to individual sequences surrounding the M4+ mainshocks, and document only a few hundred events in their respective catalogs. This study uses a dense local seismic array to build a comprehensive data set of stress drop values and focal mechanism solutions (FMSs) of HF-related earthquakes over a ∼3-year time period in the Kiskatinaw area (Figure 1), which covers parts of the Montney Formation (BC Oil and Gas Commission, 2012), extending between Fort St. John and Dawson Creek (Davies et al., 2018, purple area in Figure 1b). With the dense seismic instrumentation in this region (Figure 1b), the aim of this study is to obtain a more comprehensive view of the predominant faulting styles and stress drop variations, as well as to gain insight into the interaction between injection activity and earthquake generation.
Only a small number of HF wells in the WCSB are associated with a vigorous seismic response, and even fewer with M3+ earthquakes . Areas that do experience significant levels of HF induced seismicity share a number of common features: waveform similarity between swarm-like sequences of earthquakes, proximity of earthquake hypocenters to the wellbore and/or basement, and temporal correlation with well pad completion, particularly in the Kiskatinaw area . Previous studies estimated FMSs of moderate-sized events in this area that suggest the activation of strike-slip faults at low angles to S Hmax , as well as events of reverse faulting perpendicular to S Hmax (Babaie Mahani et al., 2017, 2019, 2020;

Earthquake Catalog and Event Classification
The earthquake catalog used in this study contains a total of 8,283 events in the Kiskatinaw area (circles in Figure 1) in the time period from 12 July 2017 to 31 July 2020 (updated from Roth et al. (2020)). Catalog detections result from an automated short-term average/long-term average (STA/LTA) trigger and analyst-reviewed phase arrivals. We refer the reader to Roth et al. (2020) for further details of the catalog development. Here we use a total of 25 broadband surface stations with a sampling rate of 100 Hz (triangles with solid edges in Figure 1) operated by the Ruhr University Bochum, McGill University, and Natural Resources Canada for event location and source parameter analysis. The stress drop estimation uses 15 additional broadband stations deployed in 2020 by University of Calgary and Geoscience BC (dashed triangles in Figure 1; Salvage et al., 2021). Of the 15 additional stations, 13 have a sampling rate of 200 Hz and 2 have sampling rates of 100 Hz. As we consider the signal up to the Nyquist frequency, signal processing does not require downsampling of the 200 Hz stations for comparison. Table S1 shows the (non-uniform) operation periods of the individual stations during the study period.
We first cluster the 8,283 catalog events based on waveform similarity using the same clustering approach detailed in Roth et al. (2020). We also cluster events in the time domain by identifying time windows with at least four events on consecutive days, leading to a total of 29 windows. We perform waveform-similarity clustering within each time window and then delve into secondary spatial clustering. The two clustering steps result in 51 groups with individual cross-correlation coefficients (CCC) ≥ 0.6, and values ranging up to 0.875. One additional (52nd) group of events exhibits a clear spatial and temporal relation, in spite of having CCC < 0.3. This particular group consists of earthquakes associated with the 30 November 2018 M L 4.5 earthquake and its 3 aftershocks with M > 3.7. The number of isolated events with CCC < 0.6 totals 1,014. We refer the reader to Roth et al. (2020) for further details of the clustering algorithm and Peña  for details on the M L 4.5 sequence.
Defining event groups as described above enables designating the largest magnitude event within each group as a representative template event to be used for moment tensor inversion. The 51 event groups provide template events with a magnitude of M L 1.5+ available for focal mechanism estimation. We also identify 10 additional events with M L 2.5+ out of the 1,014 isolated events. As we use the same methodological approach for FMS as Peña , we take the solutions of the 2018 M L 4.5 and its two largest aftershocks from Peña  for the 52nd "isolated cluster" and do not recalculate them here.

Source Parameter Inversion: Focal Mechanism Estimation
Using the 61 template events together with the three reported solutions from Peña  provides a total of 64 events for which we perform (and report) moment tensor inversion estimates. Together, the 64 events comprise the FMS catalog detailed within the following sections, which is listed in Table S2. The FMSs comprise the output of a full moment tensor inversion using broadband data from a range of 6-22 stations ( Figure 1) and the 1D homogeneous layered velocity model (Table S3 in the Supporting Information S1) used in Roth et al. (2020).
The following sections first describe the inversion procedure and the FMS catalog, and then present the inferred faulting styles from similar events.

Full Moment Tensor Inversion
We estimate the FMSs for the 61 representative earthquakes determined from clustering using the probabilistic earthquake source inversion framework Grond (Dahm et al., 2018;Heimann et al., 2018). The relative advantage of Grond for this data set is that the Bayesian framework is well-suited to small events and provides a good overview of the parameter tradeoff and uncertainties. Dahm et al. (2018) demonstrate the utility of Grond to estimate FMSs for two earthquakes of comparable magnitude with similar station coverage (M w 3.2 and M w 2.8 events with 15-19 stations distances between 34 and 94 km). Grond inverts for the optimal hypocentral location and moment tensor solutions by applying a bootstrap-technique on fully inverted waveforms and comparing them to the modeled data set. It calculates the optimal FMS based on both time and frequency domain information, as well as waveform envelopes. It can also incorporate geodetic data, when available.
The algorithm procedure first calculates a Green's function database from a regional 1D velocity profile (Table  S3 in the Supporting Information S1) and a 3D grid of potential hypocenters. It uses additional programs, Qseis (R. Wang, 1999) and Fomosto (Heimann et al., 2017), to do the computation within the source region (Dahm et al., 2018;Heimann et al., 2018). We calculate the initial Green's functions database at a grid spacing of 200 m (horizontally and vertically) for source station distances ranging from <5 up to ∼50 km. The frequency content of the synthetic Green's functions is modeled for stations with a sampling rate of 10 Hz, which permits fitting of the full waveforms at frequencies of up to 5 Hz. We use time-windows ranging from 4.2 to 5.1 s that encompasses both P and S-phases for comparison between data and synthetic waveforms. The range is dictated primarily by the source-station travel time. We review the picks to ensure phase arrivals are unambiguous and distinguishable from the noise level and exclude stations with low data quality. We then perform the waveform modeling and 10.1029/2021JB022750 5 of 29 compare with data from between 6 and 22 stations meeting the azimuthal coverage at source distances between roughly 5 and 50 km.
We perform 1-5 inversions for each earthquake, and vary inversion parameter settings in order to find the parameter set which yields the lowest global misfit. Input parameters that influence the fit quality include the filtering frequency band and the number of iterations. We test frequency bands between 0.6 and 1.6 Hz for P-waves and 0.5-1.5 Hz for S-waves, where we estimate the specific ranges by trial and error, and run a total of 26,000-28,000 iterations for each frequency band. In general, the fundamental assessment of FMS quality rests on the comparison between the observed and synthetic waveforms, which is quantitatively measured by the global misfit value. The global misfit and its mean value are typically set at a given threshold for an individual data set determined by empirical testing. Values lower than 1.0 describe a reduction of the misfit and an improvement of the model. We refer the reader to the Supporting Information for further details of the various parameters tests and fit assessments ( Figure S1 in the Supporting Information S1). In our study, all 64 FMSs, each representative of either an individual cluster or an isolated event, have minimum mean bootstrap global misfits between 0.76 and 0.26 and a total mean of 0.57. As such, the catalog better captures the full spectrum of earthquake source behavior in the Kiskatinaw area, and the solution quality can be weighted according to the misfit value. Figure S1 in the Supporting Information S1 shows the data-model agreement, the convergence curves, the global misfit value, and the FMSs solution probability for four representative earthquakes (with global fit values on the upper and lower quality ends).
The focal mechanism catalog (Table S2, Figure S2 in the Supporting Information S1) provides a graphic and numeric description of the 64 FMSs (including the three solutions from Peña ). Given that 51 of the 64 solutions stem from event families with similar waveforms, we can deduce that many of the smaller events must share similar types of deformation as the larger events for which FMSs can be inverted. We therefore use waveform similarity to identify additional events in the catalog with similar FMSs to template events. We identify all events that occur within a 5 km radius of each of the 64 catalog FMSs and have a mean CCC > 0.8 across three components of a set of three stations common to the event pairs (i.e., we calculate the CCC on nine channels). A high CCC on multiple stations is a necessary condition for a similar fault orientation, as well as a large portion of overlapping rock volume along the travel path, and therefore a strong line of evidence for similar FMSs (i.e., the hypocenters lie roughly within a maximum distance of one quarter of the wavelength of the smaller event; Got et al., 1994). We further detail the waveform similarity determination in the Supporting Information S1, and group the earthquakes that are classified as similar to the FMSs for the purposes of the discussion and interpretation of faulting processes in the Kiskatinaw area.

Results: Focal Mechanism Catalog
The 64 FMS catalog is documented in Figure S2 in the Supporting Information S1 and in Table S2. We classify the solutions following Frohlich (1996Frohlich ( , 2001 using the azimuth and plunge angles of the three principle axes, P, T, and B. A vertical orientation of each respective axis corresponds to the following deformation style: B-axis to strike-slip, P-axis to normal, and T-axis to reverse faulting (Frohlich, 1992).
Using the 64 template earthquakes in the FMS catalog, we infer an additional 3,500 FMSs based on waveform similarity between events within a family. We define waveform similarity as having a CCC-threshold of 0.8 at three representative stations as described above (see Figure S3 in the Supporting Information S1 for a 10.1029/2021JB022750 6 of 29 representative example of template and matching waveforms). We provide additional details on the requirements for similarity in the Supporting Information S1.
After classifying events based on waveform similarity, we find that the classification of the entire catalog is largely similar to the classification of the inital 64 FMSs. For example, from a total of 3,564 events with calculated (64) and inferred (3,500) FMSs, roughly 93% (3,306 events) are classified as strike-slip events (StSl), 2% (87 events) are predominantly strike-slip events with a reverse-faulting component (StSl-R), 5% (161 events) are predominantly reverse-faulting events with a strike-slip component (R-StSl), and <1% (10 events) are reversestyle events (R).  (Table S2; Álvarez-Gómez, 2014; Kagan, 2005). Each corner denotes one major faulting class ( templates that do not bear similarity to other families. Figure 2a suggests a disproportionately large number of reverse-faulting mechanisms originate from the isolated templates (yellow dots). Conversely, the majority of detections derived from the templates with multiple associated events in the same family indicate strike-slip faulting.
The panels in Figures 2b-2d show the distribution of events with respect to moment tensor components, where in contrast to (a), the circle radii correspond to the earthquake magnitude. We use the MoPaD software (Krieger & Heimann, 2012) to quantify and illustrate the isotropic (ISO; Figure 2b), pure DC; Figure 2c), and compensated linear vector dipole components (CLVD; Figure 2d; Knopoff & Randall, 1970). The panel (c) demonstrates the dominantly DC component of a majority of the solutions, particularly for the strike-slip events. There are a number of solutions with a significant proportion of CLVD motion as shown in (d), which also have varying degrees of reverse motion. The magnitude distribution reveals no obvious correlation between faulting style and local magnitude M L , nor local magnitude and non-double-couple components. We note that the largest magnitude in the data set, a M L 4.5 on 30 November 2018, has a predominantly reverse-faulting mechanism and a dominantly DC moment tensor component. Figure 3 shows the spatial distribution of estimated FMSs. Figure 3 shows the 14 of 64 solutions with a non-negligible reverse faulting component (Figure 3a), and the remaining 50 exhibit a dominant strike-slip faulting mechanism ( Figure 3b). The earthquakes in Figure 1a exhibiting varying degrees of strike-slip-reverse to dominantly reverse motion occur in close proximity to the Fort St. John Graben and the Peace River (e.g., Barclay et al., 1990;Davies, 1997;Eaton et al., 1999). In contrast, the strike-slip events to the southwest (Figure 3b) appear to occur on left-lateral faults striking at roughly 30° to S Hmax (as inferred by relocations in Roth et al. (2020)). The orientation of the assumed fault planes in both (a and b) are consistent with both the mapped faults and inferred faults based on lineations in seismicity, as well as the maximum horizontal stress S Hmax direction (Bell & Grasby, 2012). We would like to note that we do not emphasize non-DC components of the solutions, as velocity model uncertainties are documented to reduce the robustness of FMS solutions as well as create spurious non-DC solutions (Bean et al., 2014;Hardebeck & Shearer, 2002). Given that the earthquake relocations point to the reactivation of optimally oriented faults as noted above, the most straightforward explanation is of shear rupture on a 2D surface . Figure 4 shows the moment tensor decomposition of the 64 FMSs in a ternary diagram with faulting style represented by distinct symbols. We observe a trend of larger magnitude events with smaller global misfit as would be expected ( Figure S4 in the Supporting Information S1). But, there is no obvious correlation between the global    Figures 2c and 4 suggest that 76% of strike-slip events are dominated by shear DC motion of 50% or greater and a relatively low volumetric component of motion. The data show no significant indication of volumetric deformation for dip-slip nor strike-slip faulting events. Out of the 7 events with a CLVD component larger than 60%, 5 exhibit some proportion of reverse faulting motion (Figures 2d and 4).
A study by R. Wang et al. (2018) quantifying the amount of DC and non-DC components of 16 FMS of HF induced events within the WCSB show an overall similar trend. They report an overall small contribution of ISO components (<5%), a dominant proportion of DC components corresponding to ∼76%, and a CLVD component contribution of ∼20%. Other studies observe large ISO components in injection environments, suggesting fault opening and closure for normal-and reverse-faulting events, respectively, and strike-slip events consistent with shear movement (Martínez-Garzón et al., 2017). While the non-negligible CLVD components within our results may represent either some volumetric component of deformation, or successive shear ruptures on parallel structures, we can not unambiguously distinguish between the two processes. We note that the proportion of moment tensors with a high CLVD and high global misfit constitute a small proportion of events. Furthermore, the solutions with low overall misfit are generally consistent with the fault orientations implied by earthquake locations, and mapped geological structure (Berger et al., 2009;Cui et al., 2017;Davies et al., 2018;Hayes et al., 2021;Norgard, 1997;Roth et al., 2020). We therefore use inferences from earthquake locations, mapped structures, and the reasons noted above related to the uncertainties in velocity structure to interpret the events as resulting predominantly from shear failure.

Source Parameter Inversion: Spectral Analysis
We consider the 8,283 events from the catalog for spectral fitting, taking the subset of events surpassing the quality control criteria detailed in the next subsection. The objective of spectral fitting is to constrain the long period spectral amplitude and the spectral corner frequency. We then use the long-period spectral amplitude to estimate seismic moment, and together with the corner frequency, we use both parameters to estimate the static stress drop value for each event. As commonly noted in source parameter inversion studies, non-source related factors such as travel path attenuation and site effects can bias spectral corner frequency estimates. For example, single-spectrum estimates have been shown to consistently underestimate corner frequency values for small earthquakes and to bias stress drop estimates to lower values, thereby introducing an artificial scaling with magnitude (Abercrombie & Rice, 2005; Ide & Beroza, 2001;Ide et al., 2003). The difficulty in isolating the source information from the non-source-related terms in the earthquake spectrum lies in the trade off between attenuation and corner frequency during the fitting process. Using alternative approaches to spectral estimation that correct the spectral signal for attenuation in different ways can help alleviate some of the bias. For example, approaches such as the "Clustered Event Method" (CEM) correct for attenuation by using multiple stations and events to put tighter constraints on seismic quality factor, Q, estimates (Ko et al., 2012). So-called spectral-ratio methods work by exploiting common travel paths to cancel the effects of attenuation altogether (e.g., Abercrombie & Rice, 2005;Ide et al., 2003). The trade-off between fitting the attenuation and corner frequency term still exists in methods similar to the CEM approach, and the spectral-ratio approach is arguably the most robust way to isolate the source signal within a spectrum. However, the number of events for which one can estimate spectral parameters decreases between the single-spectrum and CEM approaches, and decreases drastically for spectral-ratio fitting. The relative merits between comparing the results using different approaches is that it can help identify biases in individual methods while retaining a larger number of events for which parameter estimates can be interpreted (Pennington et al., 2021).
Single-spectrum fitting produces reliable seismic moment estimates for earthquakes, as the low-frequency portion of the earthquake spectrum on which the moment estimate is based is not as hampered by attenuation effects. In order to robustly estimate seismic moment and to obtain first order estimates of the spectral corner frequency to identify the approximate frequency band of interest, we perform single-spectrum fitting (detailed below). We then refine corner frequency estimates using both the CEM and spectral-ratio fitting. As noted above, the strict criteria required to use the spectral ratio approach typically render only few event (pairs) available for analysis (Baltay et al., 2011;Prieto et al., 2006;Shearer et al., 2019). As a consistency check on both the single spectrum estimates and the spectral-ratio results, we perform the CEM estimates of moment and corner frequency. We refer the reader to the Supporting Information for the details of the CEM method. As the spectral-ratio estimates provide the best constraints on spectral corner frequency, we report and interpret the corresponding stress drop estimates and potential scaling based on spectral ratios in the Results section.

Spectral Estimation and Quality Control
We begin by screening the initial 8,283 earthquakes in the catalog to select candidates for the spectral analysis that meet a series of quality control and spectral fitting criteria. The first quality control criterion requires all waveforms to have a signal-to-noise ratio (SNR) > 3 on a minimum of 5 stations in a specific frequency band (discussed below). The minimum station number requirement helps insure a variation of less than ∼65% from the mean value for 95% of the event corner frequency and stress drop estimates (Kemna et al., 2020). The initial frequency band limits for which the SNR > 3 requirement is imposed works to cull poor-quality waveform recordings of the initial 8,383 events for which we impose additional fitting quality criteria.
We first determine the time windows to be used in the phase and noise spectra calculation. The time window determination starts with a magnitude-based estimate to ensure a sufficient window length to contain source information for the lowest stress drop events while minimizing contamination from secondary phase arrivals. The procedure determines the time window based on (the inverse of) a theoretical corner frequency estimate, f c,initial , calculated using a low-end stress drop value of 0.01 MPa. The next step refines the time window length using the f c,initial calculations to insure that it covers the relevant part of the phase energy for individual stations. Considering both the station and magnitude information helps account for dispersive effects and variability in energy recorded at variable station distances. It also helps remove potential systematic bias in stress drop estimates of small magnitude events between M 2.3 and 1.8 (Bindi et al., 2020). We follow the approach outlined by Bindi et al. (2020) and calculate the radiated phase energy as the integral of the squared velocity (Izutani & Kanamori, 2001) in order to dynamically estimate the per-station time window that encompasses a 90% (for hypocentral station-event-distance <25 km), 80% (hypocentral distance between 25 and 50 km), or 70% (for larger hypocentral distances) fraction of the cumulative energy. We refer the reader to the Text S2 and Figure S5 in the Supporting Information S1 for the specific details of the magnitude and energy dependent time window determination.
Once the appropriate time windows have been determined for the spectral estimation, we calculate the P and S-phase spectra as well as the noise spectra for determining the SNR. We use a multi-taper approach (Prieto et al., 2009) on the individual vertical components for P-waves and on the vector sum of both horizontal components for S-waves. In cases where one of the horizontal components contains data gaps, we calculate S-spectra using the remaining horizontal component with no gaps. To avoid potential bias of the curve fitting procedure toward higher frequencies, we resample each spectrum using an equally spaced, logarithmic sampling interval. We remove the instrument response to displacement using a low-frequency corner between 80% and 100% of the expected corner frequency for a given event magnitude assuming a stress drop of 0.005 MPa. We use upper corners for the instrument response removal at 40 and 45 Hz for stations in the 1E, XL, PQ networks and stations FSJ1 and FSJ2 in the EO network (with sampling rates of 100 Hz). We use upper corners of 90 and 95 Hz for remaining EO stations (with sampling rate of 200 Hz).
We calculate corresponding noise windows with the same window length, spaced one-time window before the P-arrival. If no P-arrival is detected, we place the noise window 5 s before the origin time. To be considered for the spectral analysis, we require an individual earthquake have a SNR > 3 on a minimum of 5 stations over a specified frequency band. The lower frequency limit for the SNR requirement corresponds to the lower frequency corner in which the instrument response is removed. The upper frequency limit is set to the frequency corresponding to a 25 dB-drop in the spectral amplitude Ω(f) relative to the flat part of the spectrum (Ω 0 ). We choose the upper frequency limit to ensure a sufficient number of data points in the high-frequency fall-off portion of the individual spectrum to obtain a numerically robust fit. Figure S6 in the Supporting Information S1 shows an example of a single spectrum fit with the corresponding waveforms plotted on the right. Of the 8,283 events in the initial catalog, 2,474 events meet the SNR > 3 requirement on a minimum of 5 stations for at least one phase, and are therefore considered in the spectral analysis. Although the SNR > 3 and 5-station-minimum requirements remove nearly 3/4 of the events in the catalog, they insure reliable data quality with sufficient frequency bandwidth for the spectral analysis.
In addition to SNR and station requirements, we also visually inspect event spectra to insure data quality. Manual inspection of the P-wave spectra and corner frequency fitting on the station recordings within 50 km of the 9 largest (M3+) events reveal a numerical bias in the form of apparent high-frequency oscillations. The observation suggests some of the time windows (and source-station distances) may be too short to encompass all of the isolated P-phase energy without contamination of S-wave energy. We therefore remove the 36 affected P-spectra, leaving 117 P-wave individual station spectral estimates. The remaining total of 17,763 individual P-wave spectra, and total of 16,338 S-spectra (34,101 phase spectra combined) are available for spectral fitting. We then fit the event phase spectra to obtain estimates of seismic moment and spectral corner frequency, which we in turn use to estimate the static stress drop, as described below. We note that the velocity recordings for the M L 4.5 event on 30 November 2018 are clipped on the closest station, MG01. We therefore use the integrated ground acceleration data from a co-located accelerometer in place of the velocity data MG01 ( Figure S7 in the Supporting Information S1).

Single Spectrum Fitting and Seismic Moment Estimation
The objective of fitting the phase spectra with a given spectral model is to recover estimates for the long-period spectral amplitude, Ω 0 , which is proportional to the scalar seismic moment, M 0 , and for the spectral corner frequency, f c . As noted above, while single spectrum estimates are disproportionately hampered by attenuation effects at high frequencies, the Ω 0 estimates generally provide robust M 0 estimates. We therefore use the single spectrum fitting to obtain M 0 estimates for all the three fitting methods (single-spectrum, CEM, spectral-ratio), and the initial estimates of f c . We then refine the more sensitive f c estimates with the CEM and spectral-ratio methods (while the latter will be the focus of the interpretation).
We estimate the Ω 0 and f c values for individual event spectra using a Boatwright source model (Boatwright, 1978): where the spectrum Ω(f) depends on the long period spectral amplitude Ω 0 , frequency f, source-receiver travel time t, spectral corner frequency f c , high frequency falloff rate, n, and the seismic attenuation, Q. The parameter γ controls the shape of the corner, where γ = 2 corresponds to the Boatwright source model. The value of t comes from the difference between the phase arrival time and catalog origin time. We hold Ω 0 fixed at the maximum value of the calculated spectrum. An inherent tradeoff occurs if all three parameters Q, n, and f c remain unconstrained during fitting, which may not represent real physical variations between attenuation and corner frequency. We therefore perform a grid search to find the best combination of Q and n that produce the lowest fitting residual, and then hold them constant (n = 2.5 and Q = 850) during the spectra fitting of all the events so that the f c remains the only free parameter (Please see Supporting Information Text S2 and Figure S8 in the Supporting Information S1 for details on determining the Q and n values used for fitting as well as the details of the fitting procedure, including uncertainty estimation). We compute the fit over the same lower and upper limit frequency band limits used for the SNR calculation. We use a least squares minimization on the vertical component of ground displacement for P-phases and the vector sum of the horizontal components for S-phases, and determine the fitting uncertainties using a jackknife approach (Tukey, 1958). We retain the range of Ω 0 value estimates from individual stations to determine a station averaged M 0 estimate, and determine station averaged initial f c estimates.

Cluster Event Method (CEM) Spectral Fitting
Single-spectrum fitting retains the highest number of source parameter estimates and provides stable M 0 estimates for a large number of events. However, it commonly leads to underestimation of f c values (particularly for earthquakes with M ∼ < 3) due to insufficient attenuation correction through the assumption of a constant Q.
Recent studies demonstrate lower uncertainties of f c estimates when solving for source parameters of clusters of events with an assumed constant Q within clusters using a so-called CEM (Ko et al., 2012). One reason for the reduced uncertainty might be that allowing Q to vary between clusters of events with similar source volumes may be more physically representative of crustal conditions. For example, Yu et al. (2020) point out that HF injection conditions can affect the mechanical properties and seismic Q in rock volumes close to the wellbore, as well as lead to high spatial variation. As the events in this data set can be easily classified into a number of clusters with similar waveforms (Section 2), the CEM approach is well suited to provide a consistency check on f c estimates based on spectral-ratio estimates, for which significantly fewer events qualify for analysis. Because the CEM approach is used as a consistency check on f c refinement, we have included the technical details of the approach in the Text S3 in the Supporting Information S1 for completeness, and focus on the spectral-ratio approach to corner frequency refinement. The latter approach forms the basis of the interpretation of our results.

Spectral Ratio Fitting
The spectral ratio approach is restricted to the earthquake source region and relies on co-located event pairs recorded at a common station to isolate source spectra (Hartzell, 1978;Hough, 1997). It employs the principle that the ratio between the single spectra of two events within approximately one source dimension (defined by the larger, "target", event) will result in cancellation of all non-source related terms that are convoluted with the source time function in the time domain. The smaller event is typically referred to as the "empirical Green's function" (eGf, Hartzell, 1978). We start by identifying viable candidate event pairs using the initial 2,474 earthquakes meeting the SNR and minimum station requirements for which we calculate single spectrum fits. Viable event pairs must meet criteria described below for co-location, magnitude differences, and bandwidth requirements based on the initial single-spectrum f c estimates.
The first quality requirement of event co-location insures that the travel path attenuation, site, and other non-source related terms are common to both events in the pair so that they cancel in the ratio. Hypocenter location estimates can be larger than the earthquake source dimensions, and therefore alone may not be helpful for identifying pairs co-located within one source dimension. However, events with similar source locations and fault plane solutions will have similar waveforms on stations recording at a range of azimuths. Identifying events from similar source volumes typically requires event pairs to have a high CCC (i.e., waveform similarity) in a frequency range up to the expected f c of the larger event (Abercrombie, 2015). We therefore impose both a hypocenter location and waveform similarity requirement in order to identify co-located event pairs. The eGf must have a maximum catalog hypocentral offset of 5 km relative to the to the target event, and a CCC ≥ 0.8 on at least one channel of the closest station. We calculate the CCC using the same window length as the one used for the template matching to enhance the FMS catalog, starting 0.5 s before the P-wave onset and ending 1.8 × (t S − t P ) after it. Preliminary processing steps revealed that the CCC ≥ 0.8 criterion at the closest station (which has the highest SNR) helped identify viable eGf-target pairs, particularly in the beginning of the study period (July 2017-October 2018) where fewer stations were operating simultaneously.
The second quality requirement relates to the magnitude difference of the event pair, which is used to insure the difference in f c values between the two events. A minimum magnitude difference ranging from 0.3 (e.g., Clerc et al., 2016) to 1.0 (e.g., Hartzell, 1978) helps insure resolution of the two spectral corner frequencies, and fitting of the target f c . In cases where the frequency band of resolution is wide enough, the eGf f c fit can also be constrained. In this study, we use a minimum magnitude difference of 0.5. Although a magnitude difference greater than 1.0 was required by some studies (e.g., Abercrombie et al., 2016;Hartzell, 1978;Hatch et al., 2018), our choice of 0.5 avoids preemptively removing a large number of event pairs that may have resolvable differences in corner frequencies. A number of other studies (including within the WCSB) use a similar magnitude difference in order to retain a larger number of event pairs from catalogs with large numbers of small events (Holmgren et al., 2019;. In addition, considering event pairs with a magnitude difference of 0.5 has the added benefit of using multiple eGfs, which results in multiple individual corner frequency estimates for the target event and minimizes potential bias from using a single eGf (Abercrombie, 2015). A number of additional studies have also used smaller magnitude ranges to yield significant results (Abercrombie, 2013(Abercrombie, , 2014Clerc et al., 2016;Harrington et al., 2015;Holmgren et al., 2019;Kwiatek et al., 2014).
The third quality requirement for candidate event pairs relates to additional SNR bandwidth requirements. In order to be considered for fitting, individual event SNRs must exceed 3 for the frequency band ranging from 0.75% of the single-spectrum target corner frequency estimate, f c,t , to 300% of the initial eGf corner frequency estimate f c,eGf . We set the upper fitting frequency limit at 30 Hz, as the SNR of the surface stations deteriorates at higher frequencies. The maximum 30 Hz upper boundary ensures a corner frequency resolution of 10-15 Hz, based on conservative estimates (Abercrombie, 2014). Abercrombie et al. (2017) and Ruhl et al. (2017) report an even higher cutoff criteria of 1/2 to 2/3 of the upper frequency limit of high SNR resolution, corresponding to 15-20 Hz for this data set. We therefore interpret our absolute values of stress drop for values up to 15 Hz, but scaling up to f c values of 10 Hz in alignment with more conservative estimates of bandwidth. As part of the bandwidth requirements for analysis, we use a time window length for spectral estimation of both waveforms that corresponds to the magnitude-estimated window length for the target event. The requirement implies that the quality control criteria described above must also be fulfilled for the eGf under the constraints of the target event window length.
Once we identify all possible target-eGf-pairs, we calculate the spectral ratios on individual stations meeting the above criteria (co-location, magnitude difference, and bandwidth criteria) and stack them into a single target-eGfratio for fitting. As a final quality requirement, we consider only target-eGf pairs for which at least three eGfs meet the above criteria. The reason for the last requirement is that visual inspection of spectral ratios reveals that target events with multiple eGfs generally produce smoother, more stable spectral ratios compared with target events with single eGfs.
The analytical expression for the stacked target-eGf spectral ratio (Equation 1) is as follows (Abercrombie, 2014;Abercrombie & Rice, 2005;Viegas et al., 2010): where Ω 0 and Ω 0 are the constant low frequency displacement amplitude of the target and eGf events, respectively, f is the frequency, and are the respective corner frequencies, n is the high frequency fall-off rate, and γ is a constant that controls the shape of the corner. As in the case for the single-spectrum and CEM estimates, we set γ = 2 for the Boatwright source model. The analytical expression demonstrates the effect of dividing the two individual spectra, which cancels the instrument and non-source related effects, such as site effects and travel path, leaving the ratio of the source spectra. We fit the stacked spectral ratios for the and values using Equation 2, where n is fixed at 2.5 (as in the single spectrum estimates), and the long-period ratio is held fixed at the maximum observed value of Ω 0 ∕Ω 0 within the reliable frequency band. We require a minimum of five stations to compute and fit a spectral ratio stack. We calculate the best fit using a least-squares fitting procedure for P-phases on the vertical component, and S-phases on the vector sum of the available horizontal components, respectively. As multiple eGfs will result in multiple estimates for , we calculate the mean out of all estimates and the respective 95% confidence interval with the Jackknife method, with one corner frequency estimate being randomly removed from the remaining estimates each time. Figure 5a shows a representative example of a spectral ratio estimate and corresponding fit for the same target event shown in Figure S5 in the Supporting Information S1. The target event has a M L 3.1 and the eGf event M L 1.7. The event-pair CCC value is 0.91 on the closest station and has an average of 0.84 on three selected stations. The colored lines in (b and c) show nine individual station single spectra for the target (b) and eGf (c), respectively. The individual ratios in Figure 5a shown in translucent colors are calculated from the individual spectra in (b and c), and their stack is shown as a solid black line. We perform numerical fitting of the stacked ratio in (a) using Equation 2 with and as free parameters, (and the long period ratio Ω 0 ∕Ω 0 , n = 2.5, and γ = 2 held fixed), in the frequency band highlighted in orange. The dashed black line shows the optimal fit. The solid-line circle shows the estimate and the dashed-line circle highlights the estimate. As noted above, a number of studies show that stable corner frequency estimates require a high SNR at frequencies that are 2-3 times the f c value (Abercrombie, 2014;Abercrombie et al., 2017;Ruhl et al., 2017). Many of the values in this study typically fall beyond the upper frequency resolution limit of 15 Hz. Even for some cases where f c,eGf < 15 Hz, the amplitude ratio of the low-frequency-level (LFL) to high-frequency-level (HFL) (at frequencies above the f c,eGf ), LFL/HFL, is too low to constrain f c,eGf robustly. Ruhl et al. (2017) suggest a minimum amplitude ratio of 3, while Abercrombie (2013) suggests higher amplitude ratios of 5 to ensure the shape of the spectral ratio is well fit. We therefore follow the accepted convention of reporting f c values of the target events only here, and omit reporting f c,eGf estimates that border the frequency limit of resolution (e.g., see Figure 6; Abercrombie, 2014Abercrombie, , 2021Ide et al., 2003;Shearer et al., 2019).

Estimation of f c , M 0 , and Δσ
Once seismic moment (M 0 ) and corner frequency (f c ) values have been estimated from single spectra and corner frequency values have been refined by spectral-ratio fitting, we remove estimates from stations where estimated values diverge more than two standard deviations from the linear mean. We then compute static stress drop (Δσ) values for the three types of f c fitting approaches: Single-spectra, CEM, and spectral ratio fitting. The following steps describe the seismic moment (M 0 ) and static stress drop (Δσ) calculations, where the former uses the long-period spectral amplitude from single spectra fitting (Ω 0 ) only and the latter is calculated using both Ω 0 and f c . We first estimate seismic moment for individual events and stations using the long period spectral amplitude constrained on single spectra following the relation defined in Brune (1970Brune ( , 1971):  where β is the depth-dependent shear wave velocity, ρ is the depth-dependent rock density with values between 2.46 g/cm 3 and 2.86 g/cm 3 (as reported in CRUST1.0 (Laske et al., 2013, Table S3 in the Supporting Information S1), R is the source-receiver distance, and U ΦΘ is the radiation pattern for P-and S waves, assumed to be 0.53 and 0.63, respectively (Aki & Richards, 2002). We first calculate the station mean M 0 and then estimate the eventwise confidence intervals based on a Jackknife replacement method. We calculate confidence intervals with the Jackknife method by repeating the average calculation N−1 times, where N is the number of station estimates. We then remove individual M 0 station estimates for each iteration (Prieto et al., 2007). Single spectrum M 0 estimates are then used in all stress drop estimates.
Similarly to the mean and confidence interval estimates for M 0 , we also estimate the mean f c values and confidence intervals of single-spectrum fitting. We calculate the mean values using individual values for events with a minimum of five single-spectrum fits. We determine mean f c and confidence interval values for the CEM approach by choosing 20 random model parameters of f c , Q, and n for which the misfit is individually minimized. We refer the reader to the Text S3 in the Supporting Information S1 for further information about the model parameter, optimization process, and confidence interval estimation of the CEM approach.
Because the f c values are fit from stacked spectral ratios from a minimum of 5 stations, the f c itself is already a mean value. We estimate uncertainties using the deviation in f c estimates for individual target-eGf stacks, and discard f c -estimates for which the confidence intervals exceed 20% of the mean of the f c values. We further ensure the quality of the spectral ratio fits through visual inspection, which enables manual removal of stations from the stack where the individual station ratio quality is poor (Roth et al., 2021). Visual inspection also allows narrowing the frequency fitting band if the initial corners allow the inclusion of spectral noise. For example, visual inspection reveals that estimated stress drop values of <0.1 MPa typically arise from erroneous corner frequency fits or anomalous spectral ratio shapes relative to the analytical function in Equation 2. As stated in Section 4.4, we report f c,t only, based on two factors: (a) the magnitude range in the catalog does not allow for large differences between target and eGf magnitudes, and/or (b) the f c,eGf may exceed the frequency bandwidth limitation of 10-15 Hz. The former may render the f c,eGf as being less reliably resolved from the f c,t , while the latter is a typical SNR frequency bandwidth limitation inherent to surface station data (Onwuemeka et al., 2018). We note that by removing qualitatively poor fits, we potentially neglect that some events may exhibit real source complexity. However, differentiating between noisy spectral features and real complexity effects is outside the resolution of the data set of surface stations and therefore the scope of this study.
With M 0 values from single spectra fits and f c values from all approaches in hand, we then calculate the static stress drop values of individual earthquakes. For the earthquakes with CEM and/or spectral ratio refined f c estimates, we calculate static stress drop values for each corresponding frequency estimate. We use a circular crack model (Eshelby, 1957), where r describes the radius of the rupture calculated from the radiated S-wave spectrum (Brune, 1970(Brune, , 1971). The expression for r is given by: where k is a constant that relates the spherical average of corner frequencies applied to a specific theoretical rupture model. Here, we use k = 0.38 and k = 0.26 for P-and S-wave estimates, respectively, from Kaneko and Shearer (2014), who adjusted the initial model of Madariaga (1976) to variations in dynamic models. We calculate confidence intervals for Δσ estimates for each method based on the 95% confidence interval of all station estimations of Δσ and the Jackknife remove-one algorithm.

Results: Stress Drop Estimates
In total, we estimate the spectral corner frequencies, seismic moment, and stress drop values of 1927 P-phases and 1882 S-phases by single-spectrum fitting (gray circles in Figures 6a, 6b, S9a and S9b in the Supporting Information S1), 1406 P-phases and 1117 S-phases using the CEM method (blue squares in Figures 6c, 6d, S9c and S9d in the Supporting Information S1), and 242 P-phases and 247 S-phases by spectral-ratio fitting (green diamonds in Figures 6e, 6f, S9e and S9f in the Supporting Information S1). While the frequency bandwidth resolution limits the magnitude range available for analysis, the results show consistency with the CEM and spectral ratio estimates of smaller-magnitude events in a nearby study area documented in Yu et al. (2020) (plotted in black diamonds for comparison in Figures 6e, 6f, S9c and S9f in the Supporting Information S1). Figure 6 shows the relationship between seismic moment (left) and moment magnitude (right) with corner frequency, and Figure  S9 in the Supporting Information S1 shows the scaling of seismic moment (or moment magnitude) with stress drop for P-and S-phases, respectively. Isolines of constant stress drop in Figure 6 assume a shear wave velocity of 3.2 km/s, which corresponds to the velocity of the most common hypocentral depth between 2 and 4 km.
As discussed above, a number of studies caution that corner frequency value may be underestimated if it exceeds the rough limit of ∼1/3-1/2 of the highest frequency for which the SNR exceeds the threshold of 3 (corresponding to ∼30 Hz for surface stations and a limit of 10-15 Hz). Figures 6 and S9 in the Supporting Information S1 therefore show the frequency range above 10 Hz as gray shaded areas in order to emphasize that the scaling of estimates above 10 Hz should be interpreted with caution. (An additional caveat is that the gray shaded areas are based on the isolines that assume a constant shear velocity, where the event estimates use a depth-dependent shear velocity). In order to demonstrate the robustness of the spectral ratio estimates, we scale the green diamonds by the number of eGfs per target event. Table 1 gives an overview of the mean stress drop estimates for all events, calculated by the respective method and their Jackknife confidence intervals. Entries marked with an asterisk include events with a f c < 10 Hz, and are consistent with all events meeting the quality criteria guidelines for analysis, to within estimated errors.
The dashed isolines of stress drop in Figure 6 demonstrate the theoretical slope of −3 for constant stress drop scaling, that is, they show 0 ∼ Δ −3 for a given shear velocity value (Aki, 1967). In addition, we have added dashed maroon linear regression lines showing the scaling between M 0 and f c for the population of events with f c < 10 Hz for reference.
Out of 1,927 f c estimates for P-phases, 787 fall below the 10 Hz line, as well as 1,765 out of 1,882 estimates for S-phases. The respective single-spectrum stress drop values range between 0.1 and 10 MPa for both P-and S-phases. The frequency bandwidth of resolution appears to affect the estimates at frequencies lower than 10 Hz for S-phases (Figure 6b). Figure S9 in the Supporting Information S1 shows a similar decrease of stress drop with seismic moment for magnitudes of roughly M w < 3 for both P-and S-phases. The gray area suggests that the smaller-magnitude events are primarily affected by bandwidth limitations, and the apparent scaling disappears for events with M w > 3 where SNR considerations allow f c to be well resolved. The number of estimates with higher magnitude are less numerous however, which makes broad inferences about scaling difficult. Nevertheless, the decrease in stress drop values with magnitude for events with M < ∼2.3 is commonly observed in single-spectrum estimates on surface stations (e.g., magnitudes ∼2; Lanza et al., 1999).
The same general trend of decreasing stress drop with moment for events with f c > 10 Hz persists in the CEM estimates (Figures 6c and 6d). However, the CEM seems to provide a better attenuation correction for frequencies up to 10 Hz, as demonstrated by the constant stress drop trend for both P-and S-waves. The spectral ratio estimates are less numerous but show consistent results with the CEM, and follow the stress drop isolines even above the 10 Hz value, suggesting a more effective attenuation correction up to frequencies of ∼15 Hz. The spectral ratio estimates are more limited in their magnitude range, which is the likely reason the linear regression line deviates Note. Entries with * are mean values using estimates with f c < 10 Hz, only.

Table 1 Mean Static Stress Drop Estimates With Confidence Intervals for All Respective Methods
more from the stress drop isoline than the CEM. Nevertheless, the CEM and spectral-ratio estimates provide consistent results up to the 10 Hz range, especially for S-waves.
The most reliable estimates from spectral-ratio fitting indicate 37% of events have stress drop values between 1 and 3 MPa ( Figure S9e in the Supporting Information S1). Roughly 48% of P-phase stress drop values estimated using CEM are between 0.3 and 1 MPa ( Figure S9c in the Supporting Information S1). Overall, estimates from spectral ratio fitting for P-phases are higher compared to the mean CEM values (Table 1). However, average stress drop values derived from S-phases ( Figures S9b, S9d, and S9f in the Supporting Information S1) using CEM are lower compared to spectral-ratio estimates (1.18 MPa ± 0.07 MPa using CEM and 7.07 MPa ± 2.04 MPa using spectral-ratio fitting). The overall range of stress drop values are similar between the two methods, as shown by the sidebar histograms in Figure S9 in the Supporting Information S1. Roughly 41% of spectral ratio estimates range between 1 and 3 MPa and 29% range between 3 and 10 MPa, while 49% of the CEM estimates range between 0.3 and 1 MPa, and 36% range between 1 and 3 MPa. The single-spectrum estimates do not exceed the 10 Hz resolution boundary, particularly for S-phases. Table 1 shows a similar range in values (to within error) for populations both within and outside of the conservative resolution bandwidth limits. Contrary to single spectrum estimates, both spectral ratio and CEM estimates show now obvious visible trend of increasing stress drop with seismic moment for events within the frequency bandwidth of resolution. The spectral ratio results highlight a less-pronounced scaling between seismic moment and corner frequency, implying a constant stress drop of ∼10 MPa, and a range from 1 to 100 MPa.

Focal Mechanism Solutions and Faulting Style
The FMSs constrained in this study agree with the general trend estimated in previous regional studies (Babaie Mahani et al., 2017Onwuemeka et al., 2019;Peña Castro et al., 2020), as well as with the overall observations in the WCSB (Eaton & Babaie Mahani, 2015;Schultz et al., 2017;Schultz, Mei, et al., 2015;R. Wang et al., 2018;Zhang et al., 2016). The assumed fault plane orientations are consistent with the direction of the maximum horizontal stress S Hmax (Bell & Grasby, 2012): strike-slip faults at roughly 30° to S Hmax that are consistent with left-lateral deformation, and reverse faults with nodal planes roughly perpendicular to S Hmax . The additional events considered in the extended catalog are also similar to Roth et al. (2020).
One way to provide further constraints on the ambiguity between the fault and auxiliary planes is using evidence from rupture directivity measurements. We estimate rupture directivity based on the per-station f c -variance with azimuth for 57 strike-slip events, of which 20 events have a FMS according to the FMS catalog (shown on the map in Figure 7). The FMSs suggest either right-lateral deformation on faults striking roughly N-S, or left-lateral deformation on faults striking ENE-WSW. We refer the reader to Text S4 and Figure S10 in the Supporting Information S1 for detailed information on the directivity measurements. The inset diagram in Figure 7 shows a rose diagram of the individual station estimates, with a dominant trend between 185° and 245°. The overall trend is also consistent with left-lateral motion on strike-slip faults that are optimally oriented in the maximum horizontal stress direction. The results differ from Holmgren et al. (2019), who calculated rupture directivity for 37 events within the WCSB and found rupture oriented primarily in the N-S direction. However, the events they analyzed are located further to the east in Alberta, and are also consistent with the regional stress direction there.
The relation between faulting style and hypocentral depth does not provide any obvious correlation between event type and depth distribution (Figure 8). While the low number of reverse-style events relative to strike-slip-style events hinders establishing any dominant trend, we can still examine if the distribution suggests any systematic differences. Figure 8a shows the distribution of catalog hypocentral depths and Figure 8b the depths estimated from the Grond FMSs as a function of the T-axis plunge (Figure 2). The lower plunge angles correspond to strike-slip deformation and values closer to 90° correspond to reverse-fault deformation. As Figure 2 suggests, the amount of normal faulting is not significant, therefore we consider only the two dominant fault styles in the figure.
The catalog depths do not show any distinctive systematic differences between hypocenter depths of strike-slip and reverse-style events. Conversely, the FMS-inferred hypocentral depths suggest that larger magnitude events might occur at greater depth in the basement units. While the catalog hypocentral depth of the largest event, the M L 4.5 reverse faulting event on 30 November 2018), is located at 1.5 km, relocation results from Peña Castro While the depth of the largest M L 4.5 is recorded as 1.5 km in the initial catalog, comparison of S-P arrival times to events with common epicenters suggests the depth may be up to ∼4 km (see text). The lack of viable eGfs for the M L 4.5 suggest an isolated event.   are based on a multi-station matched filter catalog and relative relocations. Because relative relocations can not provide constraints on absolute depths, we also use S-P times of earthquakes with common epicentral locations to place better constraints on the true depth of the M L 4.5. A comparison of the S-P differential times of the M L 4.5, its largest aftershock (M L 4.2, StSl), and one M L 2.1 strike-slip event that occurred in May 2018 at the same epicenter shows a greater S-P arrival time suggests a relatively deeper hypocenter of the M L 4.5 compared to the strike-slip events (Text S5 and Figure S11 in the Supporting Information S1). Another line of evidence consistent with the deeper hypocenter is that template matching does not produce any viable eGf events with similar waveforms to the M L 4.5. A lack of a viable eGf suggests the mainshock likely occurred as an isolated event. While the M L 4.5 is an isolated case, it is the largest in the data set and therefore noteworthy, particularly because the above lines of evidence suggest it ruptured a pre-existing basement fault. Its properties suggest that it resulted from a fault activation process that may be distinct from the shallower strike-slip faults.
For the M L 4.5, we estimate Δσ values of 2.86 MPa, and 5.93 MPa for P-and S-phases, respectively, using single spectra with an assumed depth of 1. The spatial distribution of the 64 FMSs highlights that a larger number of events with a reverse-faulting component are located in the northwestern part of the study area. In contrast, the strike-slip style events concentrate mainly in the central to south part of the study area and have nodal planes oblique to the east-west-striking mapped graben bounding fault (Figure 3). The orientation of mapped faults (Berger et al., 2009;Cui et al., 2017;Davies et al., 2018;Hayes et al., 2021;Norgard, 1997) are generally consistent with the observed spatial trend of FMSs and the orientation of nodal planes. The north to northwest portion of the study area is in close proximity to the Fort St. John Graben, which hosts normal faults formed during graben formation and basin infill (Barclay et al., 1990;Berger, 1994;Davies, 1997). The prevalence of reverse faulting events in proximity to the graben is consistent with reverse-sense reactivation of pre-existing normal faults under the present-day stress regime. We note that large-magnitude HF induced events have also been linked to the margins of near-basement carbonate reef structures in the neighboring Duvernay Formation in Alberta . Schultz et al. (2016) interpret the seismicity as reactivation of pre-existing basement faults associated reef growth during the Devonian. In contrast, the southeast part of the study area in this work may host more predominantly shallow strike-slip events on more recently formed faults, possibly due to a lower number of potential deep faults away from the graben boundary. Reverse faults may also be present under current stress conditions due to the thrust-faulting belt in the Rocky Mountains foreland that formed during the Laramide orogeny (Kao et al., 2018;Pană & van der Pluijm, 2015). If reverse faults are indeed present in the southern part of the study area, they do not appear to have been active during the time period of this catalog.
Another factor that can explain the activation of two different fault types in the same region may be the stress conditions. Stress constraints taken from borehole data and focal mechanism inversions in Fox Creek, Alberta suggest that S v and S hmin are similar in magnitude (Shen et al., 2019). If close in magnitude, it would be consistent with the activation of both reverse and strike-slip faulting. For example, work by Shen et al. (2019) suggests S v may be larger than S hmin at shallower depths in the WCSB, while S hmin could be larger than S v at greater depth.
If the S v transitions from being the intermediate principle stress at shallow depths to being the least principle stress at greater depth, it would favor strike-slip deformation that transitions to reverse deformation at the point where the critical stresses change orientation. Having similar values of S v and S hmin would also be consistent with overpressure conditions at depth, which tend to reduce the differential stress (Peška & Zoback, 1995;Zoback & Townend, 2001). According to Shen et al. (2019), the stress regime transition occurs between 2.8 and 17.3 km, with the optimal depth at 5.9 km near Fox Creek, Alberta ( Figure S12 in the Supporting Information S1, 95% significance). The hypocentral distribution of the catalog used here does not point to any obvious correlation of a transition at a specific depth in the Kiskatinaw area (Figure 8), at least with the resolution of hypocentral depth error alone. Examining the correlation between faulting style and relative relocation depths from enhanced catalogs does suggest that the orientation of the principle stress directions transitions below depths of ∼2 km, that is, below the concentration of strike-slip earthquakes . It is worth noting that reverse faulting events in this data set generally do not exhibit waveform similarity with many other events (Figure 2a), and are located in close proximity to the Fort St. John Graben system in the northwest portion of the study area ( Figure 3a). If the hypocenter estimates of the STA/LTA catalog are correct and reverse faulting events do occur at similar depths to strike-slip events, it would be consistent with potentially laterally heterogeneous overpressure conditions, and/or conditions where S v ≈ S hmin (Shen et al., 2019, Figure S12 in the Supporting Information S1). However, focal mechanisms alone make a transition depth difficult to pinpoint, and suggest a similar stress study to Shen et al. (2019) may be warranted.
We also perform a consistency check on the inference of similar magnitudes of S v ≈ S hmin . We perform inversions of the ambient stress field using FMSs with a global misfit <0.6 derived in this study and the software Stressinverse (Vavryčuk, 2014, Figure S13 in the Supporting Information S1). Our inversion results in a stress ratio of R = 0.9, and support the inferences outlined above from the two types of fault observed and from studies in Alberta indicating σ 2 and σ 3 as being on the same order of magnitude. The regional distribution of both strike-slip deformation (simple shear) and shortening (contraction) suggest that the area hosts a transpressive stress regime.
It is also worth considering the seismotectonic setting and geological context to further consider the generation of strike-slip style earthquakes. Figures 2 and 3 demonstrate that the majority of reverse faulting events occur on shallow to moderately dipping surfaces perpendicular to S Hmax , and strike-slip events occur on near-vertical, roughly parallel discrete slip surfaces at low angles to S Hmax Salvage et al., 2021). The strike-slip events also have a depth distribution that is predominantly uniform with the injection depths (Eaton et al., 2018;Roth et al., 2020;Riazi & Eaton, 2020;Schultz et al., 2017;Schultz & Wang, 2020, Figure 8). Results from Roth et al. (2020) suggest multiple optimally oriented parallel faults hosting left-lateral shear deformation (e.g., see the large number of strike-slip faults with parallel fault planes highlighted in detail in Figure 3b). One explanation for the observed seismicity pattern might be ruptures along Riedel shear structures, which commonly develop in the growth of strike-slip fault systems (Riedel, 1929). Riedel structures typically consist of deformation bands that are both synthetic and antithetic (Katz et al., 2004) to a major slip zone. The inset map Figure 3b highlights the scheme of a potential brittle shear zone which could host such structures.
The FMSs shown in Figure 3 would be consistent with shear on synthetic Riedel structures at low angles to S Hmax (Figure 3b), corresponding to optimally oriented strike-slip faults. Rock records of young, newly created Riedel structures suggest that Riedel shearing can also occur aseismically by micron-scale displacement in the damage zone and/or granular flow in gouge within the fault core (Pallister et al., 2012). Aseismic stress transfer has already been postulated to occur in areas of induced seismicity by both modeling work and observations, and provides a plausible explanation for why we do not observe seismicity along antithetic shears in the study region (Bhattacharya & Viesca, 2019;Eyre et al., 2019;Yu et al., 2021). In addition, Ahlgren (2001) points out that antithetic shears might also develop in further stages of faulting system growth, therefore it could be that the fault zone is still too nascent. Alternatively, the observation period may have been too short to capture any strikeslip events on larger north-northeast oriented structures for which FMSs could be computed at the low tectonic loading rates observed in the region. An additional point to note is that the scale of injection activity in the region means that virtually all of the earthquakes in our data set can be inferred to be induced. Therefore, pore pressure and stress perturbations likely affect individual faults, and/or entire fault systems, including secondary, unmapped structures (Kao et al., 2018). Potential unmapped structures might further be activated by an increase of Coulomb stress, after individual faults are activated following HF operations (Deng et al., 2016).
In general, the moment tensor decomposition suggests a high proportion of DC motion for all faulting styles and magnitudes (Figures 2-4). Independent of faulting style, the ISO component is <20% for 60 out of 64 events (Figure 4). The majority of strike-slip events have the largest DC portion of the full moment release, describing pure shear motion (Figure 4). R. Wang et al. (2018) found similar proportions for the ISO, DC and CLVD components for moment tensor solutions of induced seismicity in Fox Creek, Alberta (part of the WCSB and roughly 300 km southeast to our study area), that is, a dominant DC component ∼76%, a minor CLVD component of ∼20%, and a relatively small ISO component <4%. The authors interpret the minimal ISO component combined with the CLVD orientations along the faulting plane (northeast-southwest), as a consequence of tensile-crack opening and/or contemporaneous slip on en echelon faults. Our results generally agree with their findings (R. Wang et al., 2018). For comparison of other regions with induced seismicity, Martínez-Garzón et al. (2017) report ∼65% of analyzed injection induced earthquakes at The Geysers geothermal reservoir with non-DC components of above 25% and find the amount of non-DC motion to be correlated with the injection volume.
Out of the 64 FMSs estimated here, 9 events have a CLVD component >50% (Figure 4). We note that the percentage of CLVD component is not inversely correlated with event size. Large (vertical) CLVD components have commonly been observed in volcanic settings, and are displayed by FMSs representing vertical tension or vertical pressure resulting from slip along inward or outward-dipping ring-shaped faults around the caldera (Shuler et al., 2013). Similar "fried egg" FMSs can also be indicative of induced dip-slip/horizontal-slip earthquakes.
For example, some studies observe horizontal fracture opening along vertical structures leading to slip along the bedding-plane (Tan & Engelder, 2016). However, the CLVD component is typically lower than as observed at volcanoes (Nettles & Ekström, 1998). A third geological interpretation of a high CLVD component follows from the rupture on multiple, partially overlapping faults or fractures in time and space. When multiple ruptures occur close in time and space, it can lead to an inherent ambiguity in distinguishing multiple processes and result in higher CLVD components for inverted FMSs (Zhang et al., 2016;Lee et al., 2019). In particular, rupture on parallel oriented (Riedel) shear structures resulting from regional scale stress loading and/or larger areas of increased pore pressure could potentially cause multi-segment strike slip faulting in the southwestern portion of the study area. Figure 4 shows two additional reverse faulting events that occurred in November/December 2018 with high CLVD components. The events were part of the earthquake sequence of the M L 4.5 on 30 November 2018. As discussed in Peña , the events occurred within a complex fracture network with nearly contemporaneous rupture of multiple fault strands. Nearly simultaneous rupture of (sub-)parallel multifault strands may produce FMSs with significant CLVD components. While multi-segment rupture on a fault network is consistent with the known geological conditions, we note that interpreting FMSs, particularly of small magnitude events, should be performed with caution. Bean et al. (2008) show using synthetic tests that poor resolution of near surface velocity structure and topography can result in spurious volumetric components of moment tensor solutions, and many studies show the robustness of solutions depend on the degree to which the velocity model is known (Hardebeck & Shearer, 2002;X. Wang & Zhan, 2020). It is indeed probable that the 3D velocity structure in the region differs from 1D velocity model used here, and would affect the stability of FMS inversions (Hardebeck & Shearer, 2002). We therefore provide the caveat that the moment tensor decomposition with high non-DC and/or significant misfit values should be interpreted with caution.  Hanks, 1977). The values from our study also agree with the more broad average stress drop for small to moderate size events in the WCSB, which is 7.5 ± 0.5 MPa (Holmgren et al., 2019). However, comparing estimates for M w (or M 0 ) of the largest event in the data set, the M L 4.5, reveals deviations in the moment magnitude following P-wave and S-wave estimates:

Stress Drop and Self Similarity
We estimate a M w 4.3 using only P-phases (from 9 stations), and M w 4.6 using only S-phases (from 14 stations). Fletcher (1982) also observed higher seismic moment estimates from S-wave estimates compared to P-wave estimates in southern California. Boatwright and Fletcher (1984) provide a possible explanation for the discrepancy between P-and S-wave estimates. It rests on the fact that nodal source/receiver geometry for the P-waves and/or a set of site effects could potentially systematically amplify the S-wave motion relative to the P-wave motion. While Peña  report an overall M w 4.2, we note that other studies indicate a higher moment magnitude of M w 4.6 Babaie Mahani et al., 2019).
Although stress drop appears to decrease in magnitude below values of ∼M w 2.3, it exhibits no obvious scaling at larger magnitudes (Figures 6 and S8 in the Supporting Information S1). We interpret the difference in scaling to result from the frequency bandwidth of observation of the surface stations. Frequency bandwidth limitations disproportionately affect small-magnitude events, and the magnitude range of events analyzed here constitute mainly M w < 4 (Abercrombie, 1995;Viegas et al., 2010). A similar breakdown in constant stress drop scaling at M w ∼ 2 has been reported for a number of other studies that estimate source parameters with surface stations and single spectrum fitting (e.g., Dysart et al., 1988;Kemna et al., 2020;Onwuemeka et al., 2018;Shi et al., 1998;Viegas et al., 2010). We prioritize spectral ratio estimates in stress drop interpretations, as they provide more robust f c estimates, particularly for smaller earthquakes (e.g., Ide & Beroza, 2001;Ide et al., 2003). CEM estimates likely provide a better estimate of f c compared to single spectrum fitting, as well as spatial and temporal changes in the seismic attenuation (Figures S14 and S15 in the Supporting Information S1). The CEM results are broadly consistent with the spectral-ratio results. However, Ko et al. (2012) caution that they are strongly dependent on the pre-processing clustering because directivity effects and a sparse station dispersion can contribute to errors in both f c and Q. Again, such limitations would disproportionately affect small-magnitude (M < 3) events with higher SNR restrictions. The advantage of the CEM estimation is that by resolving all clusters individually, it is potentially able to account for 3D variation in attenuation from fluid accumulation following the HF operations.
In addition to Figures 6 and S8 in the Supporting Information S1, we provide an additional comparison of the f c -distribution for each method, as well as the f c -ratio of different methods, as a function of magnitude (Figures S16-S18 in the Supporting Information S1). The figures highlight the spectral ratio method produces higher estimates of f c , particularly for larger magnitude events within the bandwidth of resolution. We infer the higher values to result from a superior attenuation correction.
Taking the limited magnitude range into account, the stress drop values of the induced events in the Kiskatinaw are consistent with a constant stress drop scaling. The observed scaling agrees with scaling laws for self-similar tectonic events (Aki, 1967) as well as with other types of induced seismicity (Goertz- Allmann et al., 2011;Kwiatek et al., 2010). However, the general repetitive character of induced events (Skoumal et al., 2015), which we also observe in the WCSB Roth et al., 2020), may be an indicator of the rupture of co-located events on fault patches with different fault radii that generate similar waveforms. If stress drop values remain constant with magnitude and rupture velocity can be assumed roughly constant over time, it suggests that the relative proportion of average slip to fault dimension is independent of earthquake magnitude within repeating event families. Figure 9 shows the relationship between faulting style and stress drop values (the latter estimated using all three approaches). The plot contains events from families included in the FMS template-matching catalog, that is, all events with waveform similarity of CCC > 0.8 to one of the 64 representative events. We plot the stress drop values of events with M w ≥ 2.5 versus the T-axis plunge, color coded by faulting style. While stress drop estimates using the single spectrum approach contain all faulting styles, CEM is limited by clustering methods and contains exclusively strike-slip events. The spectral ratio events consist of primarily strike-slip events and few reverse faulting events. The few reverse faulting events that are included suggest reverse-faulting events tend to have similar or higher stress drop values. However, the number of non-strike-slip events is low in all estimates, therefore we can only make consistency arguments about why stress drop may correlate with faulting style.
Studies of stress drop values as a function of faulting style in North America (Boyd et al., 2017) suggest normal faulting events generally have lower stress drop values compared to strike-slip events, and that reverse (thrust) faulting events have higher stress drop values compared to strike-slip events. In contrast to the broad depth range of tectonic earthquakes in the Boyd et al. (2017) study, our source parameter data set is limited to induced seismicity within a limited depth range surrounding the horizontal injection wells. Our catalog contains 93% strike-slip events and no known normal faulting events, hindering a study of systematic differences between faulting style. Stress drop values of induced earthquakes in the WCSB seem to be homogeneous, although the number of non-strikeslip events is low if one assumes that the distribution in this study is representative of the region. While the theoretical shear strength S of reverse faults is generally higher than normal faults, with S reverse > S strike-slip > S normal (e.g., Huang et al., 2017), the shear strength is not the only factor controlling stress drop. For example, laboratory experiments correlating the fault roughness to stress drop suggest that larger stress drops occur on smooth faults compared to relatively lower values observed on rough faults (Blanke et al., 2021). A higher fault roughness on the relatively younger strike-slip fault surfaces would be consistent with both laboratory and field observations (Harrington & Brodsky, 2009;Sagy & Brodsky, 2009), and could result in lower stress drop values compared to the older and potentially smoother reactivated graben faults. Regional field scale studies provide similar observations. For example, results from a 3D seismic reflection analysis of the Osa-Peninsula megathrust subduction plate boundary in Costa Rica show heterogenities in fault surface roughness with dimensions of a few kilometers (Kirkpatrick et al., 2020). The authors conclude that spatial variation in strength and in stress drops follow the distribution of rheological asperities. In addition, the observations of Boyd et al. (2017) show lower stress drop values for induced events compared to tectonic events. Our results are also consistent with stress drop values at the mid-to-lower range of values typically observed for tectonic earthquakes. However, the values observed in the Kiskatinaw area appear to be generally lower compared to values observed in the Duvernay Play in the WCSB (Clerc et al., 2016;Zhang et al., 2016). Hence, the lack of obvious observable variations in stress drop with fault- ing style or spatial variation may be a result of seismicity being generated on both potentially smoother, mature reactivated reverse faults and younger, rougher strike-slip faults.

Conclusion
We present earthquake source parameter estimates of ∼8000 HF induced events in the Kiskatinaw area in the Montney Formation, British Columbia, between July 2017 and July 2020. A total of 64 FMSs calculated with Grond using M L > 1.5 events allow us to infer fault-planes of 3,564 out of a total of 8,283 cataloged events. A majority (93%) of the 3,564 solutions are consistent with strike-slip faulting and <1% are consistent with reverse faulting. The remainder of events suggest intermediate faulting styles. In general, FMSs imply predominantly shear motion with a low ISO component (<30%). Out of the 50 strike-slip template FMSs, 27 events have a DC component >60%. Large CLVD components may be an indicator for rupture on a system of parallel faults and/ or in some cases correspond to FMSs with a relatively large misfit. A number of factors, including the relation of FMSs to mapped faults and varying S-P arrival times, suggest that reverse faulting occurs on deeper, pre-existing normal faults reactivated in a reverse sense under the current stress regime. We interpret strike-slip faulting to occur as the result of the activation of younger, shallower faults or zones of weakness that are optimally oriented in the current ambient stress field. We note that the orientation of the regional stress field, the spatial distribution of FMSs and mapped fault distribution are consistent with strike-slip faulting resulting from slip on Riedel shear structures in a broad fault zone.
The earthquake spectral analysis suggests earthquakes exhibit constant stress drop scaling with respect to size for induced events in the Kiskatinaw area within the bandwidth of resolution, with values ranging between 1 and 10 MPa. Stress drop values are on the medium to low end of the range observed for tectonic events, and consistent with observations from other studies in the WCSB. Observation of constant stress drop within families of repeating events implies that the ratio of slip to fault dimension is constant, assuming that the rupture velocity remains constant in the same place over the time of one event family. Finally, we observe no obvious systematic variations of stress drop values with faulting style, although the small population of reverse faulting events may not permit drawing strong conclusions regarding any correlation.