Bayesian inference elucidates fault-system anatomy and resurgent earthquakes induced by continuing saltwater disposal

,

The sequence occurred near the eastern end of the Peace River Arch (PRA), an early Paleozoic emergent landmass surrounded by epicontinental seas and fringed by a Late Devonian reef system (Leduc Formation).The PRA collapsed in the Late Carboniferous, forming extensive, ancient grabens bounded by basement-rooted normal faults 16;17 .Several faults are intersected by conjugate graben structures.The M W 5.4 Dawson Creek earthquake in 2001 (green star) was a natural event with a reverse source mechanism that reactivated a normal fault in the Precambrian basement 4 .(B) Simplified Paleozoic stratigraphy of the study region 18 .
that are similar to natural events 3 , although they typically nucleate at shallower depths 4 and may exhibit spatiotemporal activation patterns consistent with an expanding pore-pressure triggering front 5 .Statistical correlation between injection rate and seismicity rate, often accompanied by a time lag of several months, has been invoked to infer causality 1;6 .Although saltwater disposal (SWD) causing injection induced seismicity (IIS) is well established in the US midcontinent 2;7;8 , IIS in western Canada has been previously attributed mainly to hydraulic fracturing 9;10;11 .Nonetheless, induced earthquakes caused by SWD have been increasingly documented in the Western Canada Sedimentary Basin (WCSB) in Alberta, including seismicity clusters near Cordel Field 12 and Musreau Lake 13;14 .Evidence that long-term injection can accelerate fault activation has also been presented 15 .
An enigmatic earthquake sequence 19 commenced in November 2022 at a previously seismically quiescent location near the town of Peace River in north-central Alberta, Canada (Fig. 1).The moment magnitude (M W ) 5.2 mainshock on November 30 was preceded by four M W > 4.0 foreshocks and was followed by a prolonged aftershock sequence, including two large events on March 16, 2023 (M W 4.5 and 4.6) (Fig. 2).
Industrial activity in the epicentral region includes widespread heavy oil production at ∼1-km depth and resultant saltwater disposal at various depths.Based on an original focal-depth estimate (6 km) greater than expected for an induced event, and the absence of precursory changes in rates of fluid disposal, Alberta's Energy Regulator initially assessed the earthquake to be a natural event 20 ; consequently, injection at nearby wells continued unabated.However the earthquakes have since been interpreted as induced primarily by long-term disposal of saltwater into a well, located about 3.5 km east of the mainshock epicentre 19 (well I in Fig. 2).If this interpretation is correct, it makes this the largest induced earthquake to date in western Canada, and implies an exceptional time lag of more than a decade between the onset of injection and fault activation (Fig. 2).
The M W 5.4 earthquake near Dawson Creek (Fig. 1) on April 14, 2001, located about 200 km west of the Peace River earthquake sequence, provides evidence of prior natural tectonic activity in the same geological setting.The 2001 earthquake is considered to be natural because it predates unconventional hydrocarbon development such as hydraulic fracturing and associated SWD, and it occurred at ∼9-km depth 4 , placing it well within the Precambrian basement.The 2001 event had a reverse mechanism 4 and likely reactivated a normal fault associated with the Late Carboniferous Fort St. John graben (Fig. 1).Crustal-scale normal faults of this age 21 are comprised of a series of secondary NW-trending normal faults (e.g., Tangent Fault, Fig. 1) that are intersected, and possibly truncated, by conjugate graben structures (e.g., Whitelaw graben, Fig. 1).The grabens are filled with syn-tectonic sediments (Stoddart Group) and the normal fault systems are step-like, characterized by a set of parallel faults that delineate structural corridors 21;22 .Collectively known as the Dawson Creek Graben Complex (DCGC), this regional fault system marks the collapse of the early Paleozoic Peace River Arch (PRA), a paleo-emergent landmass that was fringed by a massive carbonate reef complex (Leduc Formation).The fringing reef is >300 m thick and is characterized by exceptionally high permeability (up to 10 −12 m 2 ), making it well suited for SWD.Throughout most of the WCSB the Leduc Formation is underlain by older Paleozoic rocks, but in the epicentral region of the 2022 earthquake sequence it directly overlies Precambrian basement, with a discontinuous intervening unit (Majeau Lake Formation).
This study characterizes and interprets the Peace River earthquake sequence, focusing on the effects of continued SWD after the fault system was activated.The resurgence of this sequence in March 2023 correlates with an increase in injection (Fig. 2).A prolonged sequence of > 130 aftershocks between November 23, 2022 and March 23, 2023 (Fig. 3) is investigated by precise relocation analysis 23 .Using Bayesian joint inference 24 , mainshock rupture characteristics are inferred using interferometric synthetic aperture radar (InSAR) and seismic data.Source mechanisms of seven fore-and aftershocks are studied using Bayesian inference with seismic data alone.The earthquake distribution reveals a fault-system anatomy consisting of four nearly parallel NW-trending faults that are favourably oriented for slip in the present-day stress regime, truncated by a conjugate fault that is misoriented for slip.Seismic and geodetic rupture imaging of the M W 5.2 mainshock provides absolute location with low uncertainty and resolves deep nucleation, with arrest at the top of the Precambrian basement and peak slip of 0.24 m.Some areas of the imaged fault system have not yet been ruptured, pointing to possible future seismic risks.

Results
We present earthquake hypocenter relocations for 132 events 25 and detailed rupture characteristics for eight of the largest events.The mainshock is studied by joint inference of InSAR and seismic data obtained with the Bayesian Earthquake Analysis Tool (BEAT      27 .The Leduc fringing reef is >300m in thickness and has a near-vertical margin.A basement ridge is evident from the missing Majeau Lake section.(B and C) Interpretation of earthquake locations (circles, radius proportional to magnitude) in terms of four clusters (F1: orange, F2: black, F3: yellow, F4: purple) that delineate four dipping faults.The extent of the Leduc Formation (blue), Majeau Lake Formation (green), and Precambrian basement (transparent grey) are shown 27 .In (B), well locations for three wells (white hexagon), reef edge (dashed line) and an arrow indicating the view direction in (C) are shown.Specific events considered in the discussion are numbered chronologically, starting with two foreshocks (negative numbers) that are not considered with Bayesian inference (see also Supplement Table S2).
is important when combining data types to ensure that the information contributed by the various types is weighted appropriately without subjective input 26 (see Methods for details).

Reactivation of four subparallel faults
The relocation was carried out using a double-difference method 23 to refine relative hypocenter locations.
The improvement in fit is demonstrated by significant reduction of residual errors and removal of depth artifacts (see Fig. S8 and Methods for details).In map view (Fig. 4 B), aftershocks reveal an apparent NEtrending cluster of seismicity; however, a fault with this strike direction is inconsistent with inferred nodal planes (Fig. 3).When visualized in three dimensions, it is evident that the seismicity occurs along a series of NE-dipping features.By careful visual inspection, the event distribution was separated into four clusters.
The clustering was then evaluated by applying least-squares regression to fit planes to each cluster (Fig. S10).
The planar fit of events is significantly improved when considering four planes, justifying the choice of four clusters.The four fault planes interpreted from regression are labelled F1, F2, F3, and F4 (Fig. 4).These have strike values of 290 • , 289 • , 295 • , and 297 • (ordered from F1 to F4), and dip angles of 52 • , 55 • , 54 • , and 50 • .Therefore, this analysis in 3D reveals a basement fault system consisting of four subparallel faults, with seismicity mainly confined to the area covered by the Leduc fringing reef and where a basement ridge exists (Fig. 4).

Mainshock nucleates deep and arrests near top of Precambrian basement
Using functionality available in BEAT, we applied a systematic analysis to investigate the mainshock source characteristics.The mainshock is first modelled with a deviatoric moment tensor (Figs.5-A, and S1), which resolves a reverse fault and suggests that the rupture can be reasonably explained by slip on a planar fault (i.e., a double-couple, DC, mechanism).The modest compensated linear vector dipole (CLVD) component obtained by this analysis could be due to spatial extent of the ruptured asperity which is not fully explained by the point-source assumption.The predicted InSAR scenes provide a general match to the observed line-of-sight (LOS) deformation of ∼3.5 cm (towards satellite) (Fig. S2).A region of ∼-0.8 cm deformation (away from satellite) is also predicted.Variance reductions (VR) for the CMT source are (>81, >72, and >54%) indicating that the point-source approximation of CMTs cannot match all InSAR features; hence, fault models that consider spatial extent are justified.
Next, the results for a kinematic rectangular source model are considered (Figs.5-B and 6).Notable for an M∼5 event is the resolution of kinematic model parameters with low uncertainties, which is enabled by joint inference using three InSAR displacement maps and seismic data.Fault location is resolved with uncertainties of < 200 m in the lateral direction and < 100 m in the vertical direction.The model resolves strike at 310 • , dip at 45 • and rake at 118 • , with angular uncertainties of <2 • .The fault is estimated to be 2.7 km long (along strike) with 0.5-km uncertainty (95% credibility-interval width, see Methods for details) Figure 6: Inversion results for the mainshock with kinematic rectangular source constrained by InSAR and seismic data in terms of marginal and joint-marginal posterior probability densities.Note the small uncertainties and significant reduction in parameter correlations compared to inversion with only seismic data (Fig. S5). and 4.75 km wide (along dip) with 0.3-km uncertainty (Fig. 7).The average slip on this fault area is ∼0.2 m with an uncertainty of 0.06 m.The duration of moment release is 6.7 s with 0.3-s uncertainty and the rupture propagates across the fault in 1.3 s, resulting in a rupture velocity of 3 km/s which is similar to the shear-wave velocity at the rupture depths.The rupture nucleated at the bottom of the fault, displaced 0.5 km from the fault center (Fig. 5-B) at ∼5-km depth, and propagated upwards with reverse sense of motion to ∼1.9-km depth, where basement rocks are in direct contact with the disposal zone (Leduc Formation).While the kinematic solution constrains the evolution of rupture with time, it does not resolve spatial distribution of slip.
Figure 7-A shows observed and predicted InSAR data for the mainshock obtained with a kinematic rectangular fault rupture model 28;29 .For the three scenes, VRs are >82, >78, and >86%, a significant improvement over fits achieved with the CMT model and which is consistent with the notion that the more complex model explains more data features.The seismic waveforms include extensive coda (Figs.7-B, and S3), likely due to the layering structure of the WCSB.The predictions generally fit these features (Fig. S3).
For the mainshock with rectangular-fault parametrization, VR mean values of >58% for Z components and >35% for T components (Fig. S4) are achieved.Note that some stations (e.g., RV.PECRA, Fig.

Discussion Combining InSAR and seismic data reduces uncertainty
The mainshock was well recorded by three InSAR maps and seismic data, which is remarkable considering the relatively small magnitude of M W 5.2.Bayesian inference can be employed to combine these data sets in joint inference without the requirement of choosing subjective data weights (see Methods).As a result, the rupture location and characteristics exhibit low uncertainty (Fig. 6).In comparison, uncertainties for results obtained exclusively from seismic data (Fig. S5) are much larger.In addition, the results from seismic data show strong parameter correlations and multiple modes in the marginal probabilities.These effects vanish for the joint inference (Fig. 5), which furnishes evidence that complimentary information is contained in the various data sets.For example, uncertainties for fault location (including depth), fault extent, and rupture nucleation are significantly reduced.Kinematic parameters (e.g., nucleation point and rupture duration) are improved as a result of the constraints provided by the static InSAR data (e.g., fault orientation and slip amount).Since the absolute location of the mainshock has the lowest uncertainty, we applied a spatial shift to the relative hypocentre locations derived from the double-difference algorithm to achieve consistent absolute locations.

Fault-system anatomy
At first glance, the map pattern of epicentral locations (Fig. 3) suggests activation of conjugate NE-and SE-trending faults.However, clustering in 3D leads to our interpretation of four steeply dipping, subparallel faults.The seismicity on these faults concentrates along a NE-trending linear feature that we interpret as a conjugate fault that was not activated.Hence, seismicity appears to concentrate along the intersection with a misoriented fault, possibly due to stress concentrations at fault terminations 33 .This geometry is consistent with the structural style of Carboniferous faulting in the PRA region (e.g., Tangent fault and Whitelaw graben) as well as step-wise parallel sets of normal faults 34 .Both the rectangular-fault and finitefault models show that the rupture arrests at the same location where aftershocks are observed to cluster, providing further evidence of an intersecting fault at this location.
The pattern of seismicity is also consistent with fault orientations relative to the present-day stress, since the conjugate fault that truncates the NW-trending faults is misoriented for slip.Stress inversion 35 results (Fig. S11) from the eight events with M W ≥ 3.8 show that maximal compressive stress (σ 1 ) is nearly horizontal (25.3 • ±4.4 • ).The orientations of σ 2 and σ 3 exhibit low uncertainty in their orientations, but with uncertain plunges.The ratio of σ 2 to σ 3 is 0.89.Therefore, the four activated faults are favourably oriented to the current stress field 35 , and the conjugate fault is not activated because of unfavourable orientation in the current stress field.According to our model, seismicity concentrates near fault intersections and activates slip on a system of four favourably oriented, sub-parallel faults.While, carrying out stress inversion with eight similar focal mechanisms provides only limited constraints, the inferred directions are plausible and broadly consistent with regional stress.

Implications for IIS risk
The sequence activated the four faults with resurgent characteristics, likely due to the continued saltwater disposal after the mainshock.The sequence initiated on F3 with foreshocks at 4.4-and 5.3-km depth on November 23 and 24, respectively.The next two foreshocks on November 29 occurred on F2 at 3.7 and 5.6-km depth.The mainshock (No. 3 in Fig. 4) occurs on F1, rupturing from 5-km to 1.9-km depth.It is followed by two aftershocks on F3 at 3.7-(November 30) and 5.8-km (December 1) depth.We also observe that most aftershocks occur down-dip of the main rupture on F1 and on faults F2-F4.However, the M W 4.5, March 16, 2023 aftershock occurred on F1, at the SE edge of the asperity ruptured by the mainshock with nearly 4 months delay (No. 7 in Fig. 4).The resurgence was preceded by a 2-month period of quiescence and coincides with an increase in disposed volume (Fig. 2).During this resurgence, significant seismicity was only observed in this SE area with ruptures occurring on F1, F2, and F3.
The deepest fault (F4 in Fig. 4) hosts more aftershocks than F1-3, but no foreshocks.There are three events of M L > 4, all on November 30.Note that only M L estimates are available for these events which show a consistent bias of 0.86 higher than M W estimates for the catalog (Fig. S7).The centre region of F4, located beneath the disposal well, does not contain notable seismicity.This lack of fault activity suggests the possibility of future seismogenic potential.Ruptures due to natural earthquakes on similar faults to the west were shown to extend to much greater depths 4 than observed here.It is therefore possible that the faults inferred by our study may extend much deeper than the asperities activated to date.
The inferred upward migration of foreshock hypocenters, as well as the up-dip propagation of the mainshock rupture in our geometric fault model, favours a cascade nucleation model 36 that progresses from deeper to shallower levels.This 'retrograde' earthquake nucleation pattern is opposite to rupture direction that may be surmised based on naïve expectations for models of outward propagating triggering fronts 37 .
We infer that the observed updip rupture direction likely reflects trade-offs that result from superposition of the effects of downward diffusion of pore-pressure with depth-dependence of instability on a rate-state friction controlled fault 38 ; however, further work is needed to unravel the nucleation process associated with this sequence.
As noted previously 19 , the exceptional latency period of 11 years between start of injection and onset of seismicity likely reflects the prolonged time required to sufficiently pressurize faults that intersect the enormous regional saline aquifer system hosted within the Leduc fringing reef.More generally, this consideration counters possible arguments that lack of observed seismicity during multiyear operation of a disposal assures that the induced seismicity is not expected to occur in the future.In this case, direct stratigraphic contact between the disposal zone and Precambrian basement also represents a risk factor that should be avoided.

Methods Data
Three interferograms (InSAR) derived from synthetic aperture radar data acquired by the Sentinel-1 and Radarsat Constellation Mission (RCM) satellites were processed by using the GAMMA software 39 (Supplement Table S1).The topographic phase was corrected using the advanced spaceborne thermal emission and reflection radiometer digital elevation model 40 and interferograms were filtered with a Goldstein adaptive phase filter 41 .Finally, interferograms were unwrapped with a minimum cost-flow algorithm 42;43 to obtain the final deformation map employed in the inference (Fig. 7).
Seismic waveforms were rotated to radial (R), transverse (T) and vertical (Z) components and T and Z components were used in the analysis (Figs. 7 and S3).For the smallest event, 6 Z and 7 T were retained after quality control and for the largest event 17 Z and 23 T components.The waveforms were corrected for instrument response, filtered between 40-and 10-s periods, and cropped to 130 s after the predicted P-wave arrival.This processing focusses predominantly on P and S waves.Seismic amplitude spectra are based on waveforms filtered between 40 and 7 s periods and provide additional information for the same 130-s time window.Note that amplitude spectra do not consider phase information and can be considered for stations where waveform fits are poor.

Modelling of data
Data predictions for InSAR data are based on Green's functions (GFs) computed with the PSGRN/PSCMP software 44 which computes the co-and post-seismic deformation, geoid and gravity changes for a viscoelastic, layered halfspace (Table S3).Data predictions of seismic waveforms and amplitude spectra are based on GFs computed with QSEIS 30;45 and include near-field effects.To resolve source location, a 1-km spacing for GFs is computed in the source region with a 2-Hz sampling rate.All predictions assume a one-dimensional layered elastic half-space 46;47 (Supplement Table S3).

Bayesian Inference
To study the earthquake sequence, the Bayesian Earthquake Analysis Tool (BEAT) is applied 24 which utilizes nonlinear Bayesian inference with sequential Monte Carlo sampling 48 .The BEAT software is versatile and can be applied with various source models (full centroid moment tensors with unknown source time function, kinematic rectangular faults, static and kinematic finite faults) and utilize a variety of geodetic and seismic data.
In Bayesian inference, the unknown model parameters are considered to be random variables and probabilities express degrees of belief 49;50 .Therefore, uncertainty quantification is intrinsically addressed.Bayes' theorem expresses the posterior probability density for a vector of N observed data d and M model parameters m as where l(m) is the likelihood function, p(m) is the prior information, and p(d) a normalizing constant that does not need to be considered for this work.BEAT estimates the posterior p(m|d) with a sequential Monte Carlo (SMC) sampler 48 .

Joint inference and noise covariance estimation
Data errors, including theory errors stemming from the choice of Earth velocity model and source parametrization, play a crucial role for parameter and uncertainty estimates in Bayesian inference 51;52 .Some assumptions about the noise statistics are required to employ the method while others are commonly made to further simplify the application.For example, the type of statistical distribution must be assumed while the parameters for the assumed distribution may be assumed or estimated.We employ data covariance estimation that has been shown to address both errors due to the measurement process and due to theory assumptions (e.g., velocity model), it does not rely on assumptions about uncertainties that may exist in the assumed velocity model, and does not have significant computational cost 26 , as explained below.
This work considers joint inference with three data types: Seismic waveforms, seismic amplitude spectra, and deformation maps derived from synthetic aperture radar data.We assume that the noise observed on each data type is independent from the noise observed on the other types.Gaussian-distributed noise is assumed with zero mean and the noise covariance matrices are estimated from the data.Each data type consists of observations taken at various locations in space and time.For simplicity, the following considers these observations are combined into one data vector for each type.However, we will later explain specific details about the data.The mathematical treatment is sufficiently general to address these details.Under these assumptions, the likelihood function is where | • | denotes the matrix determinant, k indexes the various data types, and d k (m) are predicted data with N k being the number of data points.Estimating the covariance matrices C k is particularly important for Bayesian joint inference, where the contributions of various data sets to the solution is governed by these matrices.Therefore, each data type requires the estimation of a specific data covariance matrix.By considering eq. ( 2), the importance of the covariance matrix becomes clear since the likelihood function depends on the number of observations in the various d obs,k .For instance, waveform data typically contain many more data points than amplitude spectra, since the sampling rate is typically high.Similarly, InSAR data contain pixels at a resolution far greater than the spatial wavelengths of features under investigation.
The covariance matrices provide the appropriate statistical information that governs how much each data type contributes to the solution in eq (1).Including full covariance matrices in joint inference provides the means to quantify uncertainty objectively, without the need for human-assigned weights.
Of particular interest is the covariance estimation for InSAR data.For the mainshock, the three available scenes exhibit different resolutions.To include these scenes in joint inference, the resolution must be accounted for by estimating a full covariance matrix for each scene including variance and covariance terms.We estimate non-Toeplitz covariance matrices 53 for the spatial InSAR scenes.The method is a combination of empirical iterative covariance estimates originally developed for frequency-domain reflection coefficient data 53 and nearest neighbour search 54 via the k-nearest neighbour algorithm (KNN) 55;56 , following methods developed for potential-field data 57 .To estimate the covariance matrix, the SMC sampler initially assumes data noise to be uncorrelated with unknown standard deviation, similar to the approach implemented in BEAT for waveforms and amplitude spectra 26 .Subsequently, spatial data are binned via KNN clustering to estimate covariance terms.The approach also includes hierarchical scaling 58;59 to lower the dependence of covariance matrix estimates on the iterative scheme.
A challenge arises from SAR scenes being acquired with significant time between subsequent satellite passes (i.e.multiples of ∼12 days, Tab.S1, Fig. 7), resulting in multiple earthquakes contributing to the observed displacement fields.Acquisition times of SAR pairs are such, that the interferogram of track 020 includes only the first four inverted events, whereas both other tracks (122 and RCM) include all events except those in March 2023 (Tab.S2).Consequently, it is likely, that the residual errors (Fig. 7) include contributions from fore-and aftershocks.This is also in agreement with track 020 having the highest VRs as it includes the least events.However, the difference in M W between the mainshock and the next largest event is 0.5.In addition, the hypocenter depth of the mainshock is 1.25 km shallower than the next largest event.Therefore, such contributions are likely sufficiently small to be accounted for in the non-Toeplitz noise covariance estimate.

Earthquake source models
Source parameters are inferred by assuming point-source moment tensors 60 for the fore-and aftershocks of M>3.8.Note that the moment tensors include unknown centroid locations and source-time function parameters, resulting in highly non-linear inverse problems.In addition, moment tensors are c onsidered in Lune space 61 , which provides the significant advantage that parameter units are in degrees, making prior specification simple compared to working with fource-couples.The combination of high data quality and the availability of three InSAR scenes enables inferring parameters of spatial fault models for the mainshock.
Bayesian inference for both a kinematic rectangular fault 28;29 and a finite fault model 62;63 are considered.

Uncertainty quantification and marginalization
All parameter inferences are based on marginalization of the Bayesian posterior probability density, where the high-dimensional posterior density is simplified by integrating over all parameters except those to be interpreted by the marginal probability.For example, the marginal density for centroid depth quantifies the probability of centroid depth while integrating over all other parameter values in the posterior.This process accounts for uncertainty due to the parameters that are integrated over (Fig. S6 a) and provides more general results than sensitivity studies.It is also useful to consider the marginal probability for two parameters (joint marginals) which illuminates the correlation of two parameters (Fig. S6 a).Finally, we employ the concept of marginalization to lower-hemisphere stereographic projections of the focal mechanism ("beach-balls").We refer to these marginals of projections as fuzzy beach-balls (Fig. S6 b).Uncertainty estimates are taken to be the widths of 95% highest probability density credibility intervals estimated from marginal probabilities.

Stress inversion
Inferring the direction of principal stresses (stress inversion) is often hindered by insufficient knowledge of which nodal plane is the fault plane 64;65 .We assigned the fault plane using an iterative process 35 by selecting the plane with the highest Mohr-Coulomb failure stress within a given stress field 66;67 .We use 1000 random noise realizations in addition to the input nodal planes from the focal mechanisms for events with M W ≥ 3.8 (Fig. 1), over 100 bootstrap iterations for each inversion.This technique yields the orientation of the three principal compressive stresses (σ 1 > σ 2 > σ 3 ), as well as a relative measure of the magnitude of the intermediate principal stress (R) 68;64 with 0 ≥ R ≤ 1.Small values of R suggest that σ 1 and σ 2 are close in magnitude, whereas large values of R suggest that σ 2 and σ 3 are similar in magnitude 69 .

Earthquake relocation
An earthquake catalog 25 containing 132 events between November 23, 2022 and March 23, 2023, with the smallest event estimated to be M L 1.93, is employed to study the spatial and temporal distribution of the seismicity.The catalog locations include many solutions with fixed depths.The accuracy of the few free depths is poor because the closest two stations are just over 70 km away and the third-closest station is at 130km distance.Therefore, relocation based on picked catalog arrival times and up-sampled cross-correlations for earthquake pairs is applied.This method carries out an iterative, linearized inversion with re-weighting of data and automated discarding of poor data 23 .The final results are based on 150,000 data, of which 90.4% are retained to constrain the final result.
The relocation process reduces earthquake-pair residuals from 0.42 s to 0.23 s, cross-correlation residuals from 0.38 s to 0.05s, and absolute travel-time residuals remain similar (8.65 s and 8.63 s, Fig. S9).The residuals convey that relative locations of hypo-centres are constrained well but absolute locations for the sequence are less well constrained.For the interpretation presented here, relative locations are most important, and absolute locations are well-constrained for the largest event via the rectangular-fault solution that is constrained by InSAR and seismic data.Therefore, we apply a static, horizontal shift for all relocation epicenters based on the difference between mainshock centroid and mainshock relocation epicentre.The shift is 1.5 km to the north and 1.8 km to the east.

Data Availability
The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used to access seismic data, related metadata, and/or derived products used in this study under https://ds.iris.

Figure 1 :
Figure 1: (A) Geologic setting of the Peace River earthquake sequence in western Canada (red box) and location of the M W 5.2 mainshock (red star).The sequence occurred near the eastern end of the Peace River Arch (PRA), an early Paleozoic emergent landmass surrounded by epicontinental seas and fringed by a Late Devonian reef system (Leduc Formation).The PRA collapsed in the Late Carboniferous, forming extensive, ancient grabens bounded by basement-rooted normal faults 16;17 .Several faults are intersected by conjugate graben structures.The M W 5.4 Dawson Creek earthquake in 2001 (green star) was a natural event with a reverse source mechanism that reactivated a normal fault in the Precambrian basement 4 .(B) Simplified Paleozoic stratigraphy of the study region 18 .

Figure 2 :
Figure 2: (A) Stacked monthly injected volume of saltwater at three injection wells located near the epicenters of the Peace River earthquake sequence.Well I has operated since December 2012 and is injecting into the Leduc Formation, while wells II and III are injecting into shallower formations.The dashed box indicates the extent of data in (B) and (C).(B) Stacked monthly injected volume of saltwater at three wells near the epicenters from November 2022 to March 2023, where M w is estimated from M L for the smaller events using the relation in Fig. S7.(C) Time series of moment magnitudes (M W ) for the Peace River earthquake sequence, between November 23, 2022 and March 23, 2023.Events analysed using BEAT are shown in orange.Saltwater injection continued after the mainshock, with reported increase in March 2023.

Figure 3 :
Figure 3: Relocated seismicity (white circles) between November 23, 2022 and March 23, 2023, plotted on a ground displacement image from InSAR (background colour).Circle sizes are scaled by magnitude, and waste-water injection wells (blue hexagons) are shown.The focal mechanisms are our solutions for eight events of M W ≥ 3.8 and are numbered chronologically.The rectangle is the surface projection of a north-east dipping rectangular rupture area inferred from seismic and InSAR data, where the thick line represents the top of the fault.Displacements in the satellite line of sight (LOS) were calculated from InSAR data between November 13, 2022 and December 12, 2022.

Figure 4 :
Figure 4: (A) Surfaces showing formation tops with disposal well I (black), top of Leduc Formation (blue), top of Majeau Lake Formation (green), and top of Precambrian basement (grey)27 .The Leduc fringing reef is >300m in thickness and has a near-vertical margin.A basement ridge is evident from the missing Majeau Lake section.(B and C) Interpretation of earthquake locations (circles, radius proportional to magnitude) in terms of four clusters (F1: orange, F2: black, F3: yellow, F4: purple) that delineate four dipping faults.The extent of the Leduc Formation (blue), Majeau Lake Formation (green), and Precambrian basement (transparent grey) are shown27 .In (B), well locations for three wells (white hexagon), reef edge (dashed line) and an arrow indicating the view direction in (C) are shown.Specific events considered in the discussion are numbered chronologically, starting with two foreshocks (negative numbers) that are not considered with Bayesian inference (see also Supplement TableS2).

Figure 5 :
Figure 5: (A) Deviatoric moment-tensor solution and its decomposition into CLVD and DC for the mainshock.(B) Planar view on the rectangular fault with rupture fronts (solid black, increments in 0.1 s after nucleation time), nucleation point (black star), and the rupture extent (thick black rectangle).(C) Slip distribution for the mainshock inference from static finite-fault model.The rake direction for each patch (black arrows), 95% confidence intervals (black ellipses), and the size of the rupture inferred from the rectangular source (dashed) are shown.
7) show near-field effects in the form of long-period signals.The modelling of seismic waves with QSEIS included nearfield effects30 which match the observed near-field effects to a limited extent.Seismic amplitude spectra are generally fit more closely than waveforms (Figs.7 and S4), due to the lesser information content, specifically the lack of information about phase in the spectrum.In summary, the data fits indicate that inferred parameter values reflect the observations, which raises confidence in our results.Finally, we consider Bayesian finite fault inversion31;24  to study the spatial distribution of slip based on only InSAR data.An extended fault which employs strike and dip from the rectangular solution (fixed at the maximum posterior probability values) with square patches of 1 km 2 .Patch size was determined by sensitivity analysis32 to avoid excessive regularization.The extended size is 7 km along strike and 6 km in dip direction, with a total of 42 patches.The top of the fault is located at 1-km depth.The results (Fig.5C) provide constraints on the peak slip which is not possible with the rectangular fault.The maximum slip of 25 cm occurs at 2.5-km depth.The rupture is approximately ellipsoidal with the long axis along the rake direction (115 • ), long radius of 3 km, short radius of 2 km, and most slip occurring between 2-and 4.5-km depth.Taken together, the rectangular-and finite-fault results suggest relatively little slip where the rupture initiated at depth, with slip increasing to a maximum value at shallower depths.This finite-fault solution overall resulted in the highest VRs for the three InSAR displacement maps (>90, >89, and >82%).

Sentinel- 1 22 Sentinel- 1 22 Figure 7 :Figure 8 :
Figure 7: (A) Data (left column), MAP-parameter predictions (middle column) for the finite-fault model, and residuals (right column) for three InSAR scenes (top, middle, and bottom) used in the joint inference for the mainshock.Variance reductions and standardized residuals are shown (right-column inserts).Surface projections of rectangular fault (small) and finite fault (large) are shown (middle column) with a thick edge for the top of the faults.(B) Observed seismic waveforms and amplitude spectra for three representative examples (black).The posterior predictive density (red), standardized residual probability densities (gray), VR probabilities (yellow), station names, epicentral distances, and backazimuths are shown.Details on interferometric pairs are given in TableS1.

Figure 8 summarizes
Figure 8 summarizes Bayesian CMT results for seven earthquakes of M W > 3.8 with two 3D views.The agreement between observed and predicted data for these events is comparable to those achieved for the mainshock.Focal mechanisms for all events are highly similar oblique-reverse with rake angles varying between 100 • and 120 • (where 90 • is pure reverse).The inferred faults strike at ∼310 • and dip at angles between 37 • and 55 • .Such steep dip is atypical for reverse faults, but consistent with reactivation of normal faults within the DCGC.The 3D-view in Fig. 8-B illustrates that the moment-tensor solutions activate the complex fault architecture in several stages, with ruptures are highly consistent between faults and highest magnitudes at shallow depth.Most aftershocks are located on the deepest fault below the M5.2 earthquakes.Notably, the M W 4.5 aftershock on March 16, 2023 (No 7) appears to rupture the same asperity as the mainshock, after a 4-month delay.
edu.IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience (SAGE) Award of the National Science Foundation under Cooperative Support Agreement EAR-1851048.The earthquake catalog was obtained from Natural Resources Canada 25 .Monthly injected volume data is available through Petrinex, under Conventional Volumetric Data Download with IDs ABWI100141808217W500, 24