Machine Learning Aids Rapid Assessment of Aftershocks: Application to the 2022–2023 Peace River Earthquake Sequence, Alberta, Canada

Detecting events and picking seismic phases are fundamental to seismological workflows. Accurately identifying and characterizing seismic events and precise phase picking are critical to discern earthquake dynamics, subsurface structures, and seismic activity patterns (e.g., Cattaneo et al., 1999; Waldhauser and Schaff, 2008; Bormann, 2012; Vasyura‐Bathke et al., 2023). However, traditionally, these tasks can be time and computationally intensive using methods such as the short‐term to long‐term amplitude method (STA/LTA), manual picking, and template matching (e.g., Trnkoczy, 1999; Skoumal et al., 2015; Salvage and Eaton, 2022). The adoption of machine learning (ML) models for seismic event detection and phase picking represents a paradigm shift in seismology, offering a promising alternative to traditional approaches, saving time, and minimizing human bias (e.g., Mousavi and Beroza, 2022). The degree of automation significantly reduces routine workload, and enhances the precision and consistency of seismic event detection and phase picking, thus contributing to more robust earthquake catalogs and improving seismicity analysis and seismic hazard assessments (e.g., Tan et al., 2021; Jiang et al., 2022). In addition, the ability of ML models to analyze large‐scale and continuous seismic data streams rapidly makes them invaluable tools for characterizing and monitoring induced seismicity and developing natural sequences (e.g., Wang et al., 2020; Anikiev et al., 2023).

In particular, deep learning models, which include convolutional neural networks (CNNs) and recurrent neural networks (RNNs), have become a popular tool in seismology (e.g., Perol et al., 2018; Ross et al., 2018). Zhu and Beroza (2018) developed a deep‐neural‐network‐based arrival‐picking model (PhaseNet), which analyzes three‐component seismic waveforms to provide accurate arrival times of P and S waves through the production of probability density distributions. Woollam et al. (2019) presented a CNN for classifying seismic phase onsets for local seismic networks. With a small training dataset (411 events), this CNN‐based approach outperforms one of the classical event detection approaches based on changes in amplitude (STA/LTA). Mousavi, Weiqiang, et al. (2019) and Mousavi et al. (2020) introduced two ML‐based models, namely the CNN‐RNN Earthquake Detector (CRED) and the Earthquake Transformer (EQTransformer). The former was successfully applied to a continuous dataset recorded in Central Arkansas and tested to be efficient and promising in earthquake detection, and the latter has shown the potential for detecting and characterizing more and smaller events.

In this article, we evaluate the effectiveness of two popular ML models when applied to a temporary array deployed for the rapid assessment of aftershock sequences. We find that ML substantially improves the detection of aftershocks in terms of both number of detected events and analysis speed, compared with manual detection using a regional network of seismometers. The dataset was acquired through the Peace River Induced Seismic Monitoring (PRISM) project, jointly conducted by the University of Alberta and Alberta Geological Survey in the Peace River (PR) region of north‐central Alberta, Canada (Fig. 1), within one week of the $ML$ 5.6 earthquake on 30 November 2022 (Salvage et al., 2023; Schultz et al., 2023; Vasyura‐Bathke et al., 2023). Eight temporary three‐component short‐period stations with the sampling frequency of 250 Hz (TGS ZLAND nodes) were deployed within 15 km of the mainshock location. Continuous seismic data from 7 December 2022 to 13 January 2023, inclusive, used in this study were generously provided by the Alberta Geologic Survey (Fig. S1). This sequence represents one of the largest recorded in Alberta and is controversial in nature, because it was originally determined to be natural (Alberta Energy Regulator, 2022), but subsequent analysis has suggested that it may have been induced by nearby disposal operations (Salvage et al., 2023; Schultz et al., 2023; Vasyura‐Bathke et al., 2023).

In the realm of ML models designed to automate seismic event detection and phase picking, two widely used methods are the Earthquake Transformer (EQTransformer; Mousavi et al., 2020) and PhaseNet (Zhu and Beroza, 2018). These methods have emerged as leading tools in the field, offering significant advancements in analysis of seismicity. In a comprehensive study conducted by Munchmeyer et al. (2022), an evaluation of various ML approaches demonstrated that EQTransformer and PhaseNet outperformed most other existing models (e.g., CNN‐RNN Earthquake Detector Mousavi, Weiqiang, et al., 2019; Generalized Phase Detection Ross et al., 2018) in event detection, phase identification, and in the determination of event onset timing.

EQTransformer has demonstrated proficiency in automating seismic event detection and precise phase picking from very large datasets. Its ability to discern intricate temporal relationships within seismic waveforms enables it to detect subtle seismic signals that may prove challenging for conventional event detection methods (Mousavi et al., 2020; Munchmeyer et al., 2022). Similarly, PhaseNet has emerged as a powerful model for phase picking, exhibiting its capability to extract phase arrival times from seismic waveforms. Leveraging U‐Net architecture, PhaseNet achieves high performance in accurately identifying P and S phases (Zhu and Beroza, 2018). To associate phase picks with events, the Gaussian Mixture Model Associator (GaMMA) is often used, which assumes a hyperbolic moveout of arrival times of different phases and amplitudes (Zhu et al., 2022).

In the evaluation of both EQTransformer and PhaseNet using the Yangbi and Maduo earthquake datasets, Jiang et al. (2021) found that neural networks with deeper layers and complex structures may not necessarily enhance earthquake detection performance, suggesting that PhaseNet may perform better, in particular, in local rather than global studies. However, detailed comparison studies using both ML models on local datasets (i.e., geographically restricted) is lacking, and so utilization and analysis of more local datasets are needed for a more comprehensive understanding and further development of new models. This study employs EQTransformer and PhaseNet in direct comparison to analyze aftershocks following the $ML$ 5.6 earthquake on 30 November 2022 in Alberta.

Hyperparameter selection significantly impacts model performance, but finding the optimal configuration can vary depending on the specific model and is often computationally demanding (Wu et al., 2019). To ensure meaningful comparisons between models, we adopt fixed model architectures to focus solely on model performance by selecting (to the extent possible) identical hyperparameters (e.g., batch size, phase pick thresholds, central processing units [CPUs] etc.). This strategy enables us to concentrate on evaluating and directly comparing the effectiveness of the EQTransformer and PhaseNet models for rapid seismic event detection, providing insights into their respective capabilities. As such, we adopt the Seisbench interface (Woollam et al., 2022) for its role in standardizing the workflow and granting access to a diverse array of cutting‐edge seismological machine‐learning models and datasets. The seisbench application programming interface contains several pretrained models, including both EQTransformer and PhaseNet. In this study, both the models have undergone pretraining using the well‐established STanford EArthquake Dataset (STEAD; Mousavi, Sheng, et al., 2019)—a dataset acknowledged for its suitability in evaluating these models’ performances (Munchmeyer et al., 2022), which is one of the main aims of this study. STEAD contains 1.2 million traces, including 450,000 earthquakes (ranging from magnitude 0.5 to 8) and ∼1 million P and S picks.

As a result of pretesting trials, the hyperparameters employed in this study were chosen to balance speed and performance. The parameters used in this study are: batch size of 512; overlap of 256; P threshold of 0.55; S threshold of 0.55; and the number of CPUs as 10. After the evaluation with CPUs, we subsequently incorporated compute unified device architecture acceleration into both the models. For EQTransformer, we used a detection threshold of 0.7; all other parameters remain identical irrespective of the model architectures. After undergoing preprocessing steps, including linear detrending, down‐sampling to 100 Hz to fit the pretrained models, and high‐pass filtering at 1 Hz, the two ML models were applied to determine phase picks (Figs. S2 and S3, available in the supplemental material to this article).

For phase association we used the Gaussian Mixture Model Association (GaMMA) package (Zhu et al., 2022). For this method, each earthquake is represented as a cluster encompassing P and S phases that exhibit an approximately hyperbolic moveout of arrival times and a decline in amplitude as a function of distance. An event’s underlying distribution of phase selections is characterized using a multivariate Gaussian distribution, in which the mean values are dictated by the anticipated arrival time and amplitude derived from the seismic event responsible for the observations. GaMMA necessitates the input of an average velocity model. In this study, we adopted regional velocity averages of 4868.5 m/s for P‐wave velocity and 2863.8 m/s for S‐wave velocity (Schultz et al., 2023). To process the data effectively, we apply the DBSCAN (Ester et al., 1996) clustering option within GaMMA. This method excels at grouping data points in close proximity while efficiently identifying noise and outliers. The minimum sample parameter for DBSCAN is set to 3. In addition, the following hyperparameters were selected for filtering lower quality associations: minimum number of picks = 6; minimum P‐picks = 3; minimum S‐picks = 3; phase‐time residual (⁠$σ11$⁠) maximum value of 2.0 s; phase‐amplitude (⁠$σ22$⁠) maximum value of 1.0 $log10m/s$⁠; and covariance (⁠$σ12$⁠) maximum value of 1.0.

After event detection and phase association using the ML‐based models, we used NonLinLoc (Lomax et al., 2000, 2009) to compute hypocenters. NonLinLoc is a global search method that employs a probabilistic framework to determine event location using estimated posterior probability density functions. We use the regional velocity model by Schultz et al. (2023) for hypocenter determination. Magnitudes were calculated using the formula of Babaie‐Mahani and Kao (2020), which was developed for induced seismicity in the western Canada Sedimentary basin. Unfortunately, errors within the instrument response files meant that magnitude calculations for events detected using the local short‐period array were implausibly low (⁠$Mw$ < 0). As such, we performed a magnitude calibration (Fig. S4) using events detected in the Peace River region by Salvage et al. (2023) on a regional network of seismometers during the same time period.

EQTransformer determined 12,578 picks and detected 1,241 events, whereas PhaseNet determined 14,729 picks and detected 1,078 events (Fig. 2). Event pick quality varied across stations, and between picks from PhaseNet and EQTransformer (Figs. S2 and S3). The average processing time per station with EQTransformer was 7.63 min with 10 CPU cores, which is ∼1.8 times longer than with PhaseNet (4.07 min). When utilizing graphic processing units (GPU) acceleration, the processing time for EQTransformer and PhaseNet was 2.63 min and 2.22 min, respectively.

Both PhaseNet and EQTransformer yield hypocenters in close proximity to the mainshock on 30 November 2022, as well as to the boundary of the Leduc fringing reef and several active disposal wells in the area (Fig. 2). Such proximity was also identified by aftershocks determined using the regional monitoring network (Salvage et al., 2023; Vasyura‐Bathke et al., 2023). However, the event locations for PhaseNet appear more tightly clustered, potentially revealing subtle fault‐like structures at depth (Fig. 2c). Average errors associated with event locations detected using PhaseNet were 0.9 km (x, y direction) and 1.1 km in depth. In comparison, the hypocenters identified using EQTransformer‐detected events were significantly worse, with an average error in the x, y direction of 6.2 km and an average error in Z of 5.5 km. Artifacts relating to the implementation of NonLinLoc’s grid‐search algorithm are evident in depth, for both PhaseNet (Fig. 2c) and EQTransfomer (Fig. 2d), as revealed by several unlikely linear features (e.g., at ∼2 km depth). These artefacts appear more prominent when using EQTransformer for event detection, potentially related to larger associated errors in event location using this method. Temporally, events detected by PhaseNet and EQTransformer (Fig. 2e,f, respectively) appear moderately concurrent, with the maximum daily event counts occurring between 10 and 15 December 2022, although EQTransformer did detect events earlier in the recording period (from 7 December onward, Fig. 2f) compared with PhaseNet (Fig. 2e).

Calibrated magnitude estimates (Fig. S4) suggest similar magnitude estimates for both PhaseNet and EQTransformer, and thus similar estimated b‐values (Fig. S5). The estimated magnitudes of completeness are lower for the ML‐derived catalogs (⁠$Mc=0.9$⁠) compared with that for events detected using the regional stations (⁠$Mc=1.6$⁠). This contributes to a lower b‐value for events detected on the regional network (b= 0.854, compared with b = 1.04 for PhaseNet and b = 1.08 for EQTransformer), compounded by the fact that larger magnitude events (⁠$ML$ > 4) were identified in this catalog (Salvage et al., 2023).

To directly compare events detected by both EQTransformer and PhaseNet, we identified “common” events in both the catalogs using a time tolerance factor. This method suggests that an event in the EQTransformer catalog is assumed to be the same event as detected in the PhaseNet catalog, as long as they fall within this tolerance. With a tolerance of 1 s, a total of 838 common events were identified. The spatial distribution of common events unveils a substantial overlap (Fig. S6). For PhaseNet, the hypocenters’ latitude, longitude, and depth ranges extend from 56.03° to 56.36°, −117.21° to −116.51°, and 0.82 to 7.78 km (excluding outliers), respectively. The hypocenters’ latitude, longitude, and depth ranges for EQTransformer encompass 55.97° to 56.34°, −117.2° to −116.41°, and 0.62 to 9.75 km (excluding outliers), respectively. The more dispersed hypocenter distribution (Fig. 2b,e) for EQTransformer‐detected events could be the result of this models’ apparent greater sensitivity to noisy (in particular low frequency) waveforms. Furthermore, an apparent systematic delay in phase picks using EQTransformer (Fig. 3) likely also significantly affects the location accuracy. Reducing the tolerance to 0.5 s resulted in 794 common events.

As an illustrative example, Figure 4 provides a segment of the recorded signal from station 3 corresponding to an event detected by both the models, occurring on 19 December 2022, at 05:06:46. Although the overall detection results using these methods display a high degree of similarity (e.g., similar event counts, similar hypocentral location ranges), subtle discrepancies do arise, primarily owing to the intricate nature of waveform patterns and the inherent capabilities of these models in capturing aftershock events. In general, the EQTransformer model used in this study registers consistently delayed arrival times for both P and S arrivals in comparison with PhaseNet, with a larger discrepancy observed in the identification of P arrivals (Fig. 3). This systematic bias could stem from various factors, including the selection of hyperparameters, but is primarily attributed to the distinct architectures of the two models. This suggests that human intervention for further evaluation following the initial rapid assessment is likely still necessary at the current level of development of these models.

In Figure 5, we present events that were independently detected either only by EQTransformer (9 December 2022 at 11:24:59) or by PhaseNet (26 December 2022 at 18:34:04). The apparent high signal‐to‐noise of these events suggests that they should have been detected by both the models. This highlights a limitation of ML phase‐picker models, which sometimes fail in unpredictable ways, emphasizing the need for human intervention following rapid evaluation via ML techniques.

Although the total counts of associated events from both the models are similar, it is essential to note that these results arise from fundamentally different methodologies. EQTransformer operates in a manner that encompasses both event detection and subsequent phase picking, with the event detection preceding the phase picking stage. Consequently, all phase picks are temporally confined to an event detection window. Conversely, PhaseNet revolves around the autonomous identification of P and S arrivals, decoupled from the event detection process. In terms of uncertainty quantification, EQTransformer typically generates probabilities for event detection and phase picking by the neural network model’s output layer and represents the likelihood of the presence of specific events or phases (Mousavi et al., 2020), whereas PhaseNet converts the waveforms into probability distributions with several spikes of P and S arrivals (Zhu and Beroza, 2018). Such distinction underscores the inherent differences in the outputs of these two pickers and the nuanced nature of their respective methodologies, potentially leading to the observed differences in event detections and locations (Fig. 2). EQTransformer typically engages a longer receptive field in its architecture, which has been argued to enhance its ability to detect teleseismic events that are dominated by low‐frequency signals (Munchmeyer et al., 2022). It is possible that the less well‐resolved event locations (and in particular errors) associated with EQTransformer detections are related to the models’ higher consideration for low‐frequency signals and/or events with lower signal‐to‐noise in general.

In the context of this study, we have purposefully chosen to maintain consistency by setting the hyperparameters for both the models in a comparable way, ensuring a relatively equitable basis for comparison. However, it is crucial to acknowledge that achieving an unbiased comparison is inherently challenging due to the differences in the mathematical formulations of the two models. Despite our efforts to align parameters, the underlying mathematical intricacies might introduce nuances that impact the fairness of the comparison. Furthermore, the utilization of the initial pretraining using global seismicity in STEAD, although intentional to provide a direct comparison of methods, may not be suitable for event detection at a local scale, which are likely to exhibit different characteristics, for example, in noise. For example, the STEAD database on which our models were pretrained contains events from many different tectonic settings and with events up to magnitude 8 (Mousavi, Sheng, et al., 2019), which are not applicable for our study. However, the use of this pretraining dataset allows us to effectively compare the performance of both ML models, with the most basic of inputs. Our preliminary findings suggest that these models exhibit a high degree of consistency and efficiency in their performance, even when applied to local, rather than global, data. This underscores the adaptability and generalization capabilities of these ML models, as well as the robustness of these models when confronted with variations in seismic data characteristics.

Based on our current parameter selection for rapid assessment of aftershocks, it appears that PhaseNet outperforms EQTransformer. PhaseNet offers the advantage of shorter processing times and demonstrates smaller inferred errors in event location (Fig. 2). Focal depth and location errors are significantly larger for events identified using EQtransformer, compared with those using PhaseNet, possibly related to the inherent architecture of the model system leading to inaccuracies in phase picks. Events determined using EQTransformer do not depict clear structures in map or depth view (Fig. 2b,d), and appear to have a systematic phase pick offset compared with those determined by PhaseNet (Fig. 3). In practice, optimizing model efficiency to enhance phase picking accuracy can be achieved through approaches such as grid‐search algorithms for fine‐tuning hyperparameters, considering the spatial coherence of seismic phases at different stations, or using transfer learning algorithms (e.g., Lapins et al., 2021; Chen and Li, 2022), although such efforts are beyond the scope of this study.

The consistency of results displayed in Figure 2 implies that the seismic events are concentrated close to active disposal wells on the edge of the Leduc fringing reef. When considering event focal depths, it is evident that the predominant cluster of earthquakes aligns with $ML≥4.0$ events detected from the regional network from November 2022 through to March 2023 (gray circles in Fig. 2a–d). Events depicted from the regional catalog in Figure 2 are only those of $ML>4.0$⁠, because these have the smallest associated errors. Additional events were detected by the regional network; for further details see Salvage et al. (2023) and Vasyura‐Bathke et al. (2023). Locations obtained for events detected using PhaseNet appear to depict several linear features extending from ∼2.5 to >6 km at depth (Fig. 2c), which may reflect several subparallel fault systems that have been identified by Vasyura‐Bathke et al. (2023) in this area using aftershocks from the same sequence using the regional network.

Our findings in this study demonstrate the potential of ML for rapid seismic event assessment. The primary advantage of employing ML in this context is the capacity to process large volumes of seismic data with high efficiency. Both PhaseNet and EQTransformer identified over 1000 events within the four‐week analysis period, with average horizontal and depth errors of 0.8 and 1.1 km, and 6.2 and 5.5 km, respectively. After calibration, the minimum relative computed magnitudes were $ML$ 0.18 (PhaseNet) and $ML$ 0.29 (EQTransformer). In comparison, a standard STA/LTA algorithm with manual phase picking identified only 109 events with average horizontal and depth errors of >30 km and a minimum magnitude of $ML$ 1.02 (Salvage et al., 2023). A significant challenge of using traditional detection and picking methods (e.g., STA/LTA) is that many small‐magnitude events are often lost in noise and therefore not detected. ML models offer a superior alternative, picking events of small magnitudes with sufficient accuracy (e.g., Walter et al., 2021). The additional detected events by ML allows a greater understanding of aftershock productivity in both time and space, such as potential fault structures at depth (Fig. 2c), which were only previously identified using regional networks after event relocation (Vasyura‐Bathke et al., 2023). Furthermore, the identification of greater numbers of events using ML methods allows a better understanding of the developing seismicity rate with time (decay characteristics), which is important for understanding the evolution of stress in the subsurface. This is particularly crucial when monitoring the evolution of felt aftershock sequences (as in this case).

We applied two popular ML models for event detection and phase picking, EQTransformer (Mousavi et al., 2020) and PhaseNet (Zhu and Beroza, 2018), to continuous waveform data from a temporary aftershock deployment following the $ML$ 5.6 Peace River earthquake on 30 November 2022. Using identical (to the extent possible) choices of hyperparameters, these two methods yielded a similar number of event detections: EQTransformer detected a total of 1241 events, whereas PhaseNet detected 1078 events. For both the methods, event detections peak in the first few days after station deployment (10–15 December 2022) and hypocenters occur in close proximity to both the edge of the Leduc fringing reef and several active water disposal wells. In depth, events are clustered close to the $ML$ 5.6 mainshock event and several significant (⁠$>ML$ 4.0) fore‐ and aftershocks detected by regional stations. However, we find that hypocenter locations are more accurate (smaller errors) utilizing PhaseNet, which is able to better depict subtle fault structures at depth, compared with EQTransformer‐detected events. In general, the application of ML models shows a high degree of promise for rapid analysis of aftershock sequences, based on its ability to accurately process large data volumes with high accuracy in shorter time periods than traditional methods.

Continuous waveform data and related metadata from 7 December 2022 to 13 January 2023 (inclusive) from a temporary seismometer array were generously provided by the Alberta Geological Survey, part of the Alberta Energy Regulator (AER). The dataset was acquired through the Peace River Induced Seismic Monitoring (PRISM) project, jointly conducted by the University of Alberta and the Alberta Geological Survey. Disposal well information is available from the Alberta Energy Regulator General Well Data database https://www.aer.ca/providing-information/data-and-reports/activity-and-data/general-well-data (last accessed November 2023) by Well Identifiers 14‐18‐082‐17W5/0 (A); 06‐14‐082‐18W5/0 (B); and 14‐25‐082‐18W5/0 (C). Data processing and analysis were conducted using PhaseNet (Zhu and Beroza, 2018) and Phase Association (Zhu et al., 2022), EQTransformer (Mousavi et al., 2020), Obspy (Beyreuther et al., 2010), and NonLinLoc version 7 (Lomax et al., 2000, 2009), all of which are open source. Figures 1 and 2 were produced using PyGMT (Uieda et al., 2023). The supplemental material contains several figures denoting temporal data availability per station, phase pick differences, and magnitude calculations using the two machine learning (ML) methods. Derived catalogs of seismicity using PhaseNet and EQTransformer from this study are available at doi: 10.5281/zenodo.10436534.

The authors acknowledge that there are no conflicts of interest recorded.

This work was supported by Consortium for Research in Elastic Wave Exploration Seismology (CREWES), Consortium for Distributed and Passive Sensing, formally the Microseismic Industry Consortium (C‐DaPS), the Natural Science and Engineering Research Council of Canada (NSERC) through Grant Numbers CRDPJ 543578‐19 and ALLRP‐548576‐2019, and the Canada First Research Excellence Fund through the Global Research Initiative in Sustainable Low Carbon Unconventional Resources (GRI) at the University of Calgary. The authors gratefully acknowledge the contributions of the Peace River Induced Seismic Monitoring (PRISM) project, jointly conducted by the University of Alberta and the Alberta Geological Survey, in the acquisition and sharing of this dataset. The authors would like to thank two anonymous reviewers and the Editor, Keith D. Koper for comments and suggestions that greatly improved this article, as well as Xiaojun (Jon) Liu at the Geological Survey of Canada for reviewing an early version of this article.