Deep learning detects uncataloged low-frequency earthquakes across regions

Documenting the interplay between slow deformation and seismic ruptures is essential to understand the physics of earthquakes nucleation. However, slow deformation is often difficult to detect and characterize. The most pervasive seismic markers of slow slip are low-frequency earthquakes (LFEs) that allow resolving deformations at minute-scale. Detecting LFEs is hard, due to their emergent onsets and low signal-to-noise ratios, usually requiring region-specific template matching approaches. These approaches suffer from low flexibility and might miss LFEs as they are constrained to sources identified a priori. Here, we develop a deep learning-based workflow for LFE detection, modeled after classical earthquake detection with phase picking, phase association, and location. Across three regions with known LFE activity, we detect LFEs from both previously cataloged sources and newly identified sources. Furthermore, the approach is transferable across regions, enabling systematic studies of LFEs in regions without known LFE activity.


Introduction
Stress release on tectonic faults can happen in two ways: fast and slow.Fast deformation happens in the form of earthquakes; slow relaxation is observed as creep or episodes of accelerated slip, so-called slow slip events (Dragert et al., 2001;Ozawa et al., 2002;Lowry et al., 2001;Ide et al., 2007a).The complex interactions between fast and slow deformation might be at play during the initiation of large earthquakes (Radiguet et al., 2016;Socquet et al., 2017;Cruz-Atienza et al., 2021).However, studying these interactions requires detailed catalogs of both deformation types.While the impulsive nature of earthquakes causes clear signatures on seismic recordings, detecting slow slip is substantially more challenging.Its detection commonly uses geodetic observations with a limited spatial and temporal resolution (Michel et al., 2019;Okada et al., 2022;Costantino et al., 2023).
An alternative way to map slow deformation is by detecting and characterising its seismic markers.One such type of markers are low-frequency earthquakes (LFEs), weak seismic signals with a duration on the scale of seconds.Recent research shows that the rate and magnitude of LFEs track the slow deformation (Frank and Brodsky, 2019;Mouchon et al., 2023).LFEs are similar to regular earthquakes in some characteristics, e.g., distinct phase arrivals or predominantly double-couple sources, but also have clear differences (Shelly et al., 2007;Ide et al., 2007b;Royer and Bostock, 2014;Imanishi et al., 2016;Supino et al., 2020;Wang et al., 2023).First, they have an eponymous depletion of energy in the high-frequency band (above a few Hz).Second, in consequence of missing high frequencies, they do not exhibit impulsive arrivals but are emergent, making them hard to detect.Third, they often occur in intense bursts with inter-event times of only seconds, leading to superimposed waveforms commonly referred to as tremors (Shelly et al., 2007).
To illustrate the challenges these characteristics cause for LFE detection, it is worth contrasting LFE detection with the identification of regular earthquakes.
Detecting regular earthquakes traditionally relies on a two-step procedure: (i) phase picking, i.e., identifying P and S waves arrival times at seismic stations; (ii) phase association, i.e., selecting sets of picks across stations that are consistent with a common source location and origin time.Downstream analysis can then determine the event location and additional source parameters.In this workflow, a side benefit of the phase association step is that it acts as a quality control to remove spurious phase picks.At the moment, such a workflow is usually not applicable to LFE detection, as the low signal-to-noise (SNR) ratio and the emergent onsets make it impossible for classical algorithms to pick phase arrivals.There are exceptions, notably the JMA catalog (Japan Meteorological Agency, 2023), but these rely on high-quality and high-density data, manual intervention, and high SNR LFEs.Instead, LFE detection usually relies on manual identification (Shelly, 2010), beamforming (Frank and Shapiro, 2014), or phase coherence (Gombert and Hawthorne, 2023).These approaches often suffer from high computational demand or requirement for manual labour.However, they can be used to generate LFE template waveforms forming an initial catalog for a subsequent matched filtering search on long running recordings.As the initial approaches often fail to identify all existing LFE sources, such catalogs will be biased towards certain sources.
Due to its high sensitivity, matched filtering, also known as template matching, has become the de facto standard for LFE detection (Shelly, 2017;Bostock et al., 2015;Frank et al., 2014).Once initial templates are identified, the method identifies repeat occurrences of the template events by correlating these with the continuous waveforms.In addition to detecting occurrences, this procedure groups the LFEs into families according to their matching templates.This allows to stack waveforms and accurately locate the families.While highly sensitive, matched-filtering presents several disadvantages: templates are always region and station specific, matched filtering does not provide locations for individual events, and the model can not detect LFEs outside the initially detected families.Especially the last limitation shapes our understanding of LFEs, as template matching can only recover repeating events, potentially skewing our view of overall LFE behavior by the most repetitive sources.Furthermore, the grouping into families is partially artificial, as template matches often overlap, i.e., many detections can not be uniquely assigned to one family.
A closely linked task to the detection of LFEs is the detection of tectonic tremors.Typical methods leverage either coherency across station, through source scanning (Kao et al., 2005), waveform coherency (Armbruster et al., 2014), envelope correlations (Bombardier et al., 2023), or repetitiveness of waveform motives within or across stations (Rubin and Armbruster, 2013).However, while the underlying processes are closely related, the tasks pose distinct challenges.Tremors are usually several minutes long, making them easier to detect than LFEs.In addition, these longer waveforms make it easier to locate them, as more characteristics, e.g., envelopes, can be used for location.In contrast, LFEs are short signals, with waveforms lasting at most a few seconds, making detection and location more difficult.However, the short duration of LFEs also implies that they can monitor underlying processes at a higher resolution than tremors, thereby providing additional insights into slow deformation.In some cases, waveform coherency methods similar to tremor detection can be applied for LFE detection, but the results are usually restricted to high signal-to-noise ratio examples (Savard and Bostock, 2015).
To build a flexible LFE detector addressing the disadvantages of template matching, it would be appealing to make a more traditional earthquake detection workflow applicable to LFEs.The critical point for this is a viable automatic phase picker for LFE arrivals.We borrow from the recent breakthroughs in seismic phase picking with deep learning, where recent neural network models have substantially improved earthquake detection (Zhu and Beroza, 2019;Ross et al., 2018;Münchmeyer et al., 2022).These neural network models are trained on millions of manually labelled phase arrivals and thereby learn to accurately discern seismic phase arrivals from noise and accurately determine arrival times.The application of these models to continuous data has allowed to substantially increase the completeness of earthquake catalogs (Tan et al., 2021;González-Vidal et al., 2023;Moutote et al., 2023).
For tremor and LFE detection with deep learning, only few studies exist.Rouet-Leduc et al. (2020) identify tremor episodes in single-station records, but do not attempt to detect or locate individual LFEs.Thomas et al. (2021) focus on LFEs on the San Andreas fault and test model configurations on cataloged events.The preprint of Lin et al. (2023) presents an LFE detection workflow similar to the one presented here but focus exclusively on Southern Vancouver Island.Here, we develop a deep learning based LFE picker and show its applicability to three independent study regions: Cascadia, Guerrero and Nankai.To train our picker, we develop a novel strategy for synthetic data generation that allows for fine-grained control of the training process.Using this method, we set up a classical earthquake detection workflow and demonstrate how to automatically create LFE catalogs across different world regions.Our model successfully identifies and locates individual LFEs, even without using any training examples from the target region.The resulting catalogs are coherent with classical catalogs but have been obtained in a fully automated and region-agnostic manner.Furthermore, the catalogs identify LFEs missing from the reference catalogs, showing that our approach can uncover sources missed in the template matching procedures.We make the trained picker available with a userfriendly interface through SeisBench (Woollam et al., 2022).

Training and validation of a deep learning LFE phase picker
For detecting LFEs and determining their phase arrival times, we build a deep learning network.Our network is closely modeled after PhaseNet (Zhu and Beroza, 2019) due to the model's simplistic architecture and the excellent performance on earthquake data (Münchmeyer et al., 2022).PhaseNet is a 1D U-Net, i.e, a neural network consisting of convolutional encoder and decoder branches and skip connections (Ronneberger et al., 2015).We provide the model with 60 s of 3-component waveforms sampled at 20 Hz, bandpassfiltered between 1 and 8 Hz, the band in which LFEs are typically observed.The model outputs probability curves for P and S phase arrivals.We provide full details on the model and training procedure in the supplement.
In contrast to traditional earthquake pickers, training the model on cataloged LFE waveforms is suboptimal.First, LFEs occur in bursts, i.e., around one LFE arrival there are often further arrivals many of which have not been labelled.This leads to incorrect labelling and in addition makes a quantitative analysis of the model performance difficult.Second, most LFE catalogs are based on template matching, i.e., individual arrivals need to be inferred from arrival times on templates.Due to the low SNR, these times are often highly inaccurate, leading to high model uncertainties.Instead, we train our model on synthetics.For this, we combine LFE stacks with real seismic noise, allowing us to control the number and timing of LFEs and the SNR (Figure 1).We use up to three LFEs per trace to train the model to recognise events with low inter-event times.We use seismic noise from the INSTANCE dataset for Italy (Michelini et al., 2021) as it contains no known LFEs.
We train our model using four regions: Southern Vancouver Island in Cascadia (Canada/USA) (Bostock et al., 2015), the central section of the San Andreas fault (USA) (Shelly, 2017), Guerrero (Mexico) (Frank et al., 2014), and Nankai (Japan) (Japan Meteorological Agency, 2023).Figure S1 shows the distribution of events and stations in the reference catalogs.For Cascadia, San Andreas and Guerrero, we use template matching catalogs and the previously described strategy for generating examples.For Nankai, we apply the classical training scheme as used for earthquakes as individual picks are available.Further details on the datasets can be found in the supplement.
We evaluate our trained models quantitatively on synthetic examples generated with the previously described noise plus stack strategy.The performance on synthetic data can serve as a proxy for the expected performance on real data.We exclude the Nankai catalog from the analysis, because the catalog incompleteness precludes the extraction of challenging yet guaranteed LFE-free time windows.As this study focuses on the generated LFE catalogs, we only provide a synopsis of the analysis on synthetics here and refer to the supplement for further details.
Overall, the models show excellent detection performance for both P and S waves, with an area under the curve (AUC) of the receiver operator characteristics of 0.97 to 1.00 in all regions for positive SNR in dB scale (Figure 1).The performance degrades mildly at -2.5 dB SNR and more sharply after, but all AUC values stay above 0.88 even at -10 dB SNR.Models transfer well across regions with the worst results for a model trained exclusively on Cascadia (Figure S2).The best performing model is the one trained jointly on all four regions.Therefore, we use the model trained jointly on all regions in the subsequent analysis unless explicitly stated otherwise.
Analysing the pick time residuals, clear regional differences are visible, with lowest residuals in Cascadia (Figure S3).In all regions, average residuals are about 0.3 s larger for P arrivals than S arrivals, indicating that these are more difficult to pick.With standard deviations between 0.3 and 1.3 s (at 0 dB), residuals are substantially higher than for traditional earthquake pickers (Münchmeyer et al., 2022).Nonetheless, the residuals expose only low or no bias across all regions.For the regional differences in performance, we think that they can primarily be attributed to the heterogeneity in data quality and SNR.For example, the Cascadia stacks show the highest SNR, leading to the lowest pick residuals.In turn, this implies that no conclusions about inherent regional differences in difficulty for picking LFEs can be inferred.

Building deep learning LFE catalogs
Using our phase picking model, we set up an LFE catalog workflow similar to the classical earthquake detection workflow.Here we provide an overview of the workflow with further details in the supplement.First, we pick P and S phases by applying the trained deep learning model to continuous waveforms using SeisBench (Woollam et al., 2022).Second, we use the PyOcto phase associator (Münchmeyer, 2024) to identify coherent arrivals across stations.Third, we use NonLinLoc (Lomax et al., 2000) with a 1D velocity model to perform absolute location of the events.To avoid false detections, we filter the events based on the number of phase picks and the location residuals.For comparison, we create earthquake catalogs using the same waveforms and workflow but using a PhaseNet model trained on INSTANCE implemented in SeisBench as the picker (Zhu and Beroza, 2019;Michelini et al., 2021;Woollam et al., 2022).
As we observed a certain number of events detected as both earthquakes (EQs) and LFEs, we remove these events from the LFE catalogs (Figure S4).We note that it is not clear whether these events should be classified as LFEs or EQs.The level of overlap is dependent on the region, with almost no overlap in Cascadia and Guerrero, a 5% overlap on the San Andreas fault, and a 40% overlap in Nankai.While we are not certain what causes this different behavior, we speculate that in Nankai, LFEs and earthquakes show a wide range of apparent spectra, due to the diverse event distribution and frequency dependent attenuation.This might lead to a higher number of overlapping detections.
We apply our workflow to compile LFE catalogs for the four study regions.As we focus on studying the performance of the model and its resulting catalogs, we restrict ourselves to short study periods: Figure 2 shows the spatial event distributions.While the overall event locations are scattered, we notice strong similarities with the reference catalogs.For Cascadia (10211 events detected), LFE activity is distributed along a band underneath South Vancouver island.For Guerrero (876 events detected), LFEs occur mostly in a band between 100 • and 99 • West and around 18.25 • North.For Nankai (2525 events detected), a clear band of LFEs is visible in Southwestern Nankai.Further LFEs around 135.5 • E match the second band of LFEs com- In the reference catalogs for Cascadia (Bostock et al., 2015) and Guerrero (Frank et al., 2014)  monly observed in Nankai.For San Andreas (975 events detected), the new catalog deviates from the previous observations, with the detection broadly distributed in space instead of along the fault (Figure S7).This is likely caused by very poor locations due to the station geom-etry.As many events are only detected by the Parkfield borehole array with very dense station spacing, the aperture is small.Together with high pick uncertainties, this makes determining accurate locations challenging.Therefore, we will exclude San Andreas from the following analysis.
In all catalogs, the event depth exposes high scatter.Nonetheless, the largest density of events is around the previous estimates of LFE source depths.For Cascadia, events within the network show less scatter on depth than outside the network.We suggest that this is caused by the high timing uncertainty of the picks.In particular, a high P pick uncertainty will cause poor depth constraints as the P to S time is indicative of depth.Template matching catalogs alleviate this problem by locating LFE families instead of single events.In Guerrero, we observe an arc-shaped depth distribution.This is most likely related to the station distribution that traces out almost a straight line, leading to poor location constraints perpendicular to the station line (Frank and Shapiro, 2014).
Figure S5 shows the event density of the detected events.This visualisation further highlights the match with the reference catalogs both in terms of latitude and longitude and in terms of depth.For Japan, the highest event density occurs in the South-Western band of LFE activity.For Cascadia, the fine structure of detected LFEs is compatible with earlier publications, e.g., clear overlap is visible with patches A and C in the visualisation of LFE density of Figure 7 by Savard and Bostock (2015).
Even though the overall shape of the catalogs is consistent with the previous catalogs, this alone does not confirm that the identified events are indeed LFEs.We therefore conduct additional analysis into the newly obtained catalogs.Figure 3 shows spatial and temporal patterns in the catalogs.In Cascadia, LFEs in all three observed sequences show a clear North-Westward migration.This is consistent with the slow slip and tremor migration patterns in the area during these episodes (Wech, 2021).We conducted an additional analysis, comparing the PNSN tremor catalog and LFEs detected using our method for a 31 day period in 2021 (Figure S6).This analysis shows a high agreement between LFE and tremor locations, density, and temporal development.This holds even though for this time period we used a less dense station coverage and considered a larger part of Vancouver Island than in the 2003 to 2005 episodes.
For Nankai, we observe a migration in the North-Eastward direction.Notably, the migration is not continuous but rather has a gap and an additional, earlier cluster at the far North-West.This pattern matches exactly the migration pattern in the JMA catalog.We do not observe clear spatial migration patterns in Guerrero, however, such patterns have previously only been identified with very precise location estimates (Frank et al., 2014).In all regions, the evolution of daily event rate between the deep learning and reference catalogs is highly similar with Pearson correlation coefficients between 0.74 and 0.93.In absolute numbers, the deep learning method detects substantially fewer events than the template matching, but more events than the manual detection procedure of the JMA.We note that the number of events from deep learning is highly dependent on the chosen quality control parameters, which we set rather conservatively to avoid false positive detections.
Figure 4 shows a comparison of the velocity spectra of LFEs detected by our method, earthquakes and noise in the three regions.The spectra clearly show the characteristics of the different event types.The earthquake spectra show increasing or at least constant energy up to about 10 Hz.In contrast, the LFEs show a continuous decrease or at most constant levels of energy from low frequencies onward.The LFEs only show substantially higher energy than the noise in a small frequency band, while the EQs show substantially higher SNR at high frequencies.This depletion in energy at higher frequencies is the key property of LFEs.

Increased diversity of LFE sources through deep learning
We compare the detected events to the reference catalogs (Figure S8).In Cascadia, for 64% of the LFEs from our workflow, we find an LFE in the reference catalog within 10 s; for Guerrero for 81% of the events.For Nankai, only 8% of our LFEs are in the reference catalog, however our catalog also substantially exceeds the original catalog in the total number of events.Conversely, we recover 39% of the cataloged events.Note that a loose threshold for matching is justified due to uncertainties in the origin times due to inaccurate locations.On one hand, these results are another confirmation that the method correctly identifies LFEs.On the other hand, the substantial fraction of uncataloged events suggests that our method identifies previously unidentified LFEs.In the following, we verify and analyse these detections.First, we rule out spurious detection.To this end, we scramble the picked arrival times of each station by applying small random shifts.We choose constant shifts for each station for every hour.This destroys the exact times, while keeping the pick distribution, the P to S times within a station, and the higher number of picks within tremor bursts intact.We then build "catalogs" by associating these scrambled pick times, using the same associator settings as for the actual catalogs.The scrambled "catalogs" only contain about 5% to 10% of the number of events contained in the original catalogs and show no spatial coherence (Figure S10).Even these numbers are still likely an overestimation of the false positive rate, as events recorded at many stations are likely to be unperturbed by our scrambling procedure.Therefore, at least 20% additional detections can not be attributed to spurious associations.
Mapping the events without matches in the reference catalog (Figure S11) reveals that they follow the same spatial extent and migration pattern as the full catalog.Notably, for Cascadia and Mexico there are changes in the temporal patterns.For Cascadia, the newly detected events concentrate early and late in the LFE sequence, coinciding with a spatial location around the southeastern tip of Vancouver island and towards the northwestern end of the LFE cluster.Nonetheless, there are additional detections dispersed throughout the whole region including the central region with good coverage in the reference catalog.For Mexico, the largest frac- tion of new detections clusters in time between days 30 and 50 of the analyses sequence.Visualising the interevent time (Figure 4) confirms these observations.Both the full deep learning catalogs and the catalog of events without a match in the reference catalogs show clear burst behaviour.In particular for Mexico, certain LFE bursts are contained virtually completely in the template matching catalog, while others have not been identified at all.This highlights that the newly detected LFEs do not only uncover new sources but even new LFE bursts.
To further validate this finding, we correlate the uncataloged detections with the family stacks from the reference catalogs (Figure S12).For Cascadia, the distribution of correlation values for these uncataloged detections are identical to the noise distribution, i.e., the new detections do not match known sources.For Guererro, some events show systematically higher correla-tions than the noise.Nonetheless, the vast majority of our additional detections do not match the known sources any better than the noise.This verifies that these new detections are systematically different and that these events can not be found with template matching without identifying further, novel templates.
Extending upon the finding that the model can generalise from known families to LFEs outside these families, we investigate the ability to detect LFEs in regions the model has not been trained on.For this analysis, we trained leave-one-out models, i.e., models trained on all but one region, and applied them on the left-out region.Figure S13 visualises the spatial and temporal migration patterns.Again, the clear migration patterns in Nankai and Cascadia are retrieved.Furthermore, the number of events correlates highly (Pearson correlation between 0.69 and 0.82) with the reference catalogs.The total size of the catalogs varies, with a substantially Figure 4 Spectra and interevent times in the different regions.The top part shows velocity power spectral density for noise, earthquakes (EQs), and LFEs detected using deep learning.For each region, all traces stem from one reference station (Cascadia -MGCB, Guerrero -MAXE, Nankai -IIDH).Noise example have been extracted outside tremor episodes.Spectra have been calculated from the horizontal components (20 s windows for noise, 11 s windows starting 1 s before the S arrival for events).Thin lines show individual spectra, bold lines median spectra.EQs were selected at a similar distance and depth range to the LFEs.Network averaged spectra are shown in Figure S9.The bottom part shows the development of interevent times during the LFE sequences in the reference catalogs, the deep learning catalogs, and for the unmatched events, i.e., all events from the deep learning catalogs that are not in the reference catalogs.Vertical stripes in the events indicate the occurrence of LFE bursts.For Cascadia, we only visualise the 2005 sequence for simplicity.We visualise all events from the reference catalogs without further declustering, leading to very low interevent times.
smaller catalog in Cascadia, a similarly-sized catalog in Nankai, and a far bigger catalog in Guerrero.However, these might be related to changes in the model confidence values rather than their actual quality as we produced all catalogs with fixed picking thresholds.The cross-regional analysis clearly illustrates that the mod-els can be transferred across regions and recover LFEs from families they have not been trained on.

Conclusion
Our analysis shows that deep learning and template matching are complementary in the way they detect LFEs with specific advantages and disadvantages for either method.The biggest strength of deep learning is the flexibility.The model can be applied to additional stations, including temporary stations, shows higher diversity in terms of event families, and can be transferred across regions.In addition, our method directly allows to locate single LFEs, even though with substantial uncertainties in terms of depth.In contrast, template matching requires a predefined set of sources that is difficult to obtain and specific to each region and set of stations.While rigid, this leads to a more sensitive model, as evidenced by higher event counts.Furthermore, it allows template matching to identify LFEs with fewer stations than deep learning.For LFEs, commonly no individual location is performed after template matching, due to the difficulties caused by the low SNR ratio.However, it has been shown that relative locations of individual LFEs can be determined with classical methods as well (Shelly et al., 2009).A promising avenue might be the combination of deep learning and template matching, i.e., using deep learning to identify a diverse set of templates and afterwards use template matching to increase the completeness of the identified families.
Lastly, the deep learning method extends our view of LFEs by detecting previously unidentified sources.Building a comprehensive set of templates for template matching is challenging: the low bandwidth and SNR makes it difficult to distinguish between closely spaced sources, leading to a trade-off between missing sources and redundant templates.In contrast, the deep learning method is source-agnostic, i.e., no selection of sources needs to be performed for detecting individual events.Such a source-agnostic view is necessary to perform unbiased subsequent analysis that requires a complete view of LFE sources, such as estimates of slow slip.In addition, the fact that the model can be transferred across regions shows that LFEs have universal, regionindependent properties similar to earthquakes.Given our results, we expect that deep learning methods will allow to map LFEs across world regions with high consistency and diversity.

Figure 1
Figure 1 Data generation procedure and evaluation results for synthetic data.The top panels show (top to bottom): the combination of two LFE stacks from Cascadia, a 60 s noise segment from INSTANCE, the combination of signal and noise at 3 dB SNR, and the Gaussian pulse labels for the P and S arrivals.The bottom panels show the receiver operating characteristics (ROC) at different SNR.The numbers in the legend indicate the area under the ROC curve (AUC).For all plots, we use the joint model trained on all four datasets.

Figure 2
Figure 2 LFE catalogs obtained from deep learning (top row) and the reference catalogs (middle row).For deep learning, each dot represents one LFE.The bottom subpanels show depth cross-sections, showing longitude and depth of events.Color encodes event depth.The histograms on the left of the cross-section show the depth distribution of the detected events.In the reference catalogs for Cascadia(Bostock et al., 2015) and Guerrero(Frank et al., 2014)  each dot represents an LFE family.For Nankai (Japan Meteorological Agency, 2023) individual LFEs are plotted.The bottom panels show waveforms of associated LFE picks from deep learning in each region.Red lines indicate phase picks (dotted for P, solid for S).
Figure 2 LFE catalogs obtained from deep learning (top row) and the reference catalogs (middle row).For deep learning, each dot represents one LFE.The bottom subpanels show depth cross-sections, showing longitude and depth of events.Color encodes event depth.The histograms on the left of the cross-section show the depth distribution of the detected events.In the reference catalogs for Cascadia(Bostock et al., 2015) and Guerrero(Frank et al., 2014)  each dot represents an LFE family.For Nankai (Japan Meteorological Agency, 2023) individual LFEs are plotted.The bottom panels show waveforms of associated LFE picks from deep learning in each region.Red lines indicate phase picks (dotted for P, solid for S).

Figure 3
Figure 3 Spatial and temporal migration patterns of detected LFEs.Each dot represents one LFE.The time within the sequence is indicated by colour.The bottom panel shows the number of events per day for the LFEs from deep learning (blue) and reference catalogs (orange).The numbers in the upper left corners indicate the Pearson correlation coefficient between the daily number of events between the two catalogs.