A Silent Build-up in Seismic Energy Precedes Slow Slip Failure in the Cascadia Subduction Zone
AA Silent Build-up in Seismic Energy Precedes Slow Slip Failurein the Cascadia Subduction Zone
Claudia Hulbert , , Bertrand Rouet-Leduc , Paul A. Johnson Laboratoire de G´eologie, D´epartement de G´eosciences, ´Ecole Normale Sup´erieure, PSL Research University,CNRS UMR 8538, Paris, France Los Alamos National Laboratory, Geophysics Group, Los Alamos, New Mexico, USA ∗ C. Hulbert (email: [email protected]), B. Rouet-Leduc ([email protected]), and P. Johnson ([email protected])
Abstract
We report on slow earthquakes in Northern Cascadia, and show that continuous seismic energy in the subduc-tion zone follows specific patterns leading to failure. We rely on machine learning models to map characteristicenergy signals from low-amplitude seismic waves to the timing of slow slip events.We find that patterns in seismic energy follow the 14-month slow slip cycle. Our results point towards arecurrent build-up in seismic energy as the fault approaches failure. This behaviour shares a striking resemblancewith our previous observations from slow slips in the laboratory. a r X i v : . [ phy s i c s . g e o - ph ] S e p ince their discovery in Japan at the turn of the millennium [
1, 2 ], slow earthquakes and associated tectonictremor and low frequency earthquakes (LFEs) have been identified in most subduction zones as well as othertectonic environments [ ]. Slow slip is transitional between the fast rupture of regular earthquakes and stablesliding along a fault interface [ ], and occurs deep in faults, down-dip from the nucleation zone of damagingearthquakes at the transition from brittle to ductile zones [ ]. In subduction tectonics, slow slip and tremor arethought to take place where temperatures drive dehydration of subducting material that increase pore pressures,inhibiting brittle failure [ ]. Slow earthquakes are characterized by short-term slip with durations of days orlong-term slip with durations of months or years [ ]. The slowly slipping region is thought to denote the transitionfrom unstable (seismogenic) to stable (creep) friction and therefore may define the depth limit of megathrustrupture [ ]. How the slow slip may couple to the locked region is an open and fundamental question [
2, 9, 11 ].Slow slip has been shown to occur before some large seismogenic ruptures in some instances, for example beforethe M9.0 Tohoku earthquake [ ] and the M8.2 Iquique earthquake [ ], and has been observed preceding largeearthquakes in other tectonic environments including some transform faults [
8, 14 ]. These observations suggestthat coupling between slow slip and large earthquakes occur in some instances but has not been systematicallyobserved, perhaps due to limitations of geodetic measurements.In laboratory studies of slow slip [ ] generated by a bi-axial shear device [ ], we showed that the energyof continuous seismic waves follows characteristic patterns throughout the slip cycle, allowing us to estimate keyproperties of the laboratory fault (instantaneous fault displacement, friction, gouge layer thickness). In a first effortto generalize these results to Earth, the analysis of slow slip in Cascadia [ ] revealed that statistical characteristicsof continuous seismic signals encode the fault displacement rate, as measured by GPS. These characteristics arerelated to signal energy and resemble those identified in the laboratory.In the laboratory, continuous energy patterns also encode information regarding the future behavior of the fault(time remaining before the beginning and the end of the next slow slip event, duration and ‘magnitude’ of the nextslow slip event). Here, we seek to determine whether these same signatures can be used to infer the future behaviorof the slowly slipping zone in Cascadia beneath Vancouver Island.2 eismic data analysis and slow slip failure times in Cascadia We analyze slow slip beneath Vancouver Island, Canada, where the Juan the Fuca oceanic plate subducts north-westly beneath the North American plate (Figure 1). The Cascadia subduction zone exhibits slip events that arelong-term, on the order of weeks [
6, 20–22, 22–24 ]. Large slow slip events occur quasi-periodically (approxi-mately every 14 months), manifest by the North American plate lurching southwesterly over the Juan de Fucaplate. Smaller slow slip events occurring between these large periodic events have been identified recently, point-ing towards a large variability in the size and timing of slow earthquakes in the area [ ]. We focus on Cascadiain scaling from the laboratory because of the recurrent slow slip events observed, especially in the vicinity ofVancouver Island [ ], and because there have been continuous seismic recordings for more than a decade. Su-pervised machine learning (ML) used in the work described below, requires robust training and testing sets withmany slip events. The long history of slow slip observed in this region makes it an ideal case to determine if thereis information carried in the seismic signal characteristic of upcoming failure.Figure 1: (A) Map of the area analyzed in Cascadia and sketch of the subduction zone. The seismic stationsare represented by black triangles.The shaded region represents the area of interest, the southern Victoria Island.Our goal is to rely exclusively on continuous seismic waves (B) to estimate the time remaining before the next slowslip event. (C)
Slow slip timing from PNSN Tremor Logs (gray shaded areas), and smoothed catalogued tremorrates from the PNSN catalog for comparison.ML analysis of seismic data is an expanding field, with many studies over the recent years focusing in particular3n event detection [ ], phase identification [ ], phase association [
27, 28 ], or patterns in seismicity [ ]. Here,by relying on a supervised ML approach, we assess whether continuous seismic waves carry the signature of howclose the system is to failure. We use the PNSN Tremor Logs as a proxy for slow slip timing. These timings areestimated from bursts of tremor, and are represented by vertical gray bars in Figure 1 (C); tremor rates from thePNSN catalog [ ] are also plotted for comparison. Note that because the tremor logs take into account the wholeCascadia region and are not geographically limited to Vancouver Island, slip timings may not match perfectlyactual timings within the area of interest. We could use peaks in tremor rates within the area instead to identifyfailure times; the analysis below remains robust if we do so (see Supplementary). However, because the catalogonly starts in 2009, this leads us to discard about a third of the data. We can also use GPS data to try to identifyslow slip timing (see Supplementary).We rely on seven seismic stations from the Plate Boundary Observatory [ ], denoted by black triangles inFigure 1 (A). We find that borehole stations are much more robust than surface ones regarding contamination byseasonal signals. A continuous, clipped and de-sampled seismic waveform for one of the stations is shown inFigure 1 (B). Our goal is to assess whether a given time window of this continuous seismic data can be used toestimate the time remaining until the next slow slip event.We use the methodology initially developed in laboratory shear experiments and modified in previous work [ ]for application to Earth data. We first process the seismic data by correcting for the instrument gain, centering andde-trending the data, for each day during the period analyzed (2005-2018). The daily data intervals are band-passed between 8 and 13 Hz (within bands of 1Hz) and clipped, to limit the contribution of microseismic noiseand earthquakes respectively, and to focus on low-amplitude signals. Once the data are pre-processed, for each daywe compute a number of statistical features linked to signal energy (see Supplementary Information for details).These daily features are then averaged within a time window. Anomalous datapoints are detected within eachwindow and removed before averaging using Isolation Forests [ ]. The results shown below use a time windowof 3 months, but our methodology is robust to changes in the window size (see Supplementary). Each window isindexed by its latest day; successive time windows are offset by one day. The averaged features over these time4indows are used as input to the machine learning algorithm. In the following, ‘seismic features’ will refer tothose features averaged over a time window.As in any supervised machine learning problem, our analysis is divided into a training and a testing portion.In the training phase, the algorithm takes as input the seismic features calculated from the first 50% of the seismicdata (training set), and attempts to find the best model that maps these features to the time remaining before thenext slow slip event (label or target). The problem is posed in a regression setting. We select the model’s hyper-parameters through Bayesian optimization by 5-fold cross-validation. We use Pearson’s correlation coefficient(CC) as the evaluation metric.Once a model is finalized, it is evaluated on data the model has never seen–the features from the remaining50% of the seismic data (testing set). It is important to note that in the testing phase, the model only has access tothe seismic features calculated from the continuous seismic data, and has no information related to slow slip timing(the label). In the testing phase, the label is used exclusively to measure the quality of the model’s estimates, i.e. how close these estimates are compared to the true label values obtained from PNSN Tremor Logs. Estimating slip failure times from continuous seismic waves
We rely on gradient boosted trees, algorithms that are transparent in their analysis in contrast to many other meth-ods. These algorithms can be probed to identify which features are important in the model predictions, and why.Identifying the important statistical features allows us to compare with laboratory experiments, and gain insightinto the underlying physics.Estimations of the time remaining before the next slow slip event on the testing set are shown in Figure 2(A). This plot shows the ML slip timing estimations (in blue) from station B001 (3-month windows) and the timeremaining before the next slow slip event (ground truth, dashed red line). This ground truth can be understood asa countdown to the next slip, and is equal to zero whenever the PNSN reported an ongoing slow slip. The samemethod can also be applied successfully when using peaks in tremor rates or GPS to identify our proxies for failuretimes (see Supplementary). 5hen close to failure, the machine learning model is able to estimate the timing of most slow slip events - withthe exception of the 2018 slow slip, which is not estimated before it begins. Each point on the blue prediction curvein Figure 2 (A) is derived from a single time window of seismic data. Thus, the results demonstrate that at anytime during the slow slip cycle, a snapshot of continuous seismic waves is imprinted with fundamental informationregarding the time remaining before the next slow slip event. In particular, the seismic data long before failure andclose to failure appears very different to the trained model: although estimations far from failure are noisier, themodel easily distinguishes between these two extreme cases in all the examples in our testing set.Our results suggest that the system follows specific patterns leading to failure. Note that because the eventsare most often separated by 13 or 14 months, the model cannot rely on seasonal signals to make its estimations.To prove this is the case, we show that a seasonal sinusoidal cannot lead to good estimations of failure times (seeSupplementary). Therefore, the seismic data analyzed contains identifiable precursory patterns that match slowslip timing and are independent from the season. T i m e t o f a il u r e ( d a y s ) Date I Q - - H z , B ( m / s ) Experimental runtime (s) T i m e t o f a il u r e ( s ) A c o u s t i c p o w e r ( a . u . ) BC D
1e 10
1e 10
1e 10
Time to failure (days)
Date
CC=0.56 T i m e t o f a il u r e ( d a y s ) Time to failure (s)
ML estimatesTime to failure Time series of features A -1000100200300400 Figure 2:
Slip timing estimations, patterns identified by our model, and comparison to shear experiments inthe laboratory. (A)
ML estimates in testing of failure times (in blue), and time to failure (in red) from the PNSNtremor logs. (B)
The most important feature identified by our model plotted against time, for the best stations(B001), for time window intervals of 3 months (black curve) [The right-hand vertical axes show the interquantilerange (IQ 60-40) for the frequency band 8-9 Hz, and its evolution over a time window]. These energy-basedfeatures show clear patterns with respect to the time remaining before the next slow slip event (left axis). (C) shows the best statistical feature found in laboratory slow slip experiments for comparison (acoustic power, righthand vertical axis). The best features in the laboratory ( C ) and Cascadia ( B ) are related to the energy of the seismicwaves, and appear to follow very similar patterns, with a progressive increase in amplitude that peaks near the endof each slow slip event. (D) Distribution of these features with respect to time to failure, in between slip events -energy increases when failure approaches. Top: Cascadia slow earthquakes. Bottom: laboratory slow slips.6 atterns in seismic energy follow the slow slip cycle, with a silent build-uppreceding failure
Because we choose to rely on transparent machine learning algorithms, we can identify the most important featuresused by our model to make its estimations of failure times, and therefore make comparison with laboratory exper-iments. We find that the best features identified in Cascadia follow almost identical patterns as those identified inthe laboratory.The most important seismic features identified by our model for the best station are plotted in Figure 2 (B).The evolution of these features over time follows clear 14-month patterns, and is related to the timing of the slowslip cycles - which explains why our model is able to make its estimations of slow slip timing.We know from recent work that similar continuous seismic features can be used to accurately estimate thedisplacement rate of the fault in both the laboratory [
15, 32 ] and in Cascadia [ ]. In particular, these statisticalcharacteristics of the seismic data can be used to differentiate between the loading and slipping phases of the cycle.The estimations of failure times in Figure 2 (A) show that, just as in the laboratory, continuous seismic signalsare also encoded with information regarding the future state of the fault. This suggests that some of the frictionalphysics scales from the simple laboratory system to subduction in Cascadia, which argues in favor of self-similarityof slow slip rupture and nucleation.In the laboratory, cycles in seismic energy seem to be associated to displacement on the fault. Therefore ourresults suggest that slow slips often begin with a silent build-up in energy possibly associated to an acceleration onthe fault, that can be small enough to not be captured in the cataloged tremor - tremor at this early stage may be tooweak to be picked up systematically. We find that a model trained only on cataloged tremor is not able to estimateslip timings (see Supplementary). In contrast, the peak energy measured during slow slips are likely driven byidentified tremor events.In the case of laboratory slow slip, the most important feature by far for forecasting failure time is the seismicpower [ ], shown in Figure 2 (C). In the case of the Cascadia subduction zone we rely on inter-quantile ranges,closely related to the root of the seismic power but with outlier values removed (these outliers are more likely to7orrespond to anthropogenic noise and/or local and distant earthquakes).Figure 2 (B), top, shows that seismic energy follows similar patterns in both the laboratory and Cascadia: i)a progressive increase, with peaks in energy reached toward the end of the slow slip (as the slipping phase ofthe cycle emits larger seismic amplitudes than the loading phase, due to episodes of tremor); and ii) an oftenabrupt decrease within each slip cycle, towards the end of an event. Cycles in Cascadia are clearly apparent, butmuch noisier. This is likely due to background noise and to the fact that many fault patches located over a broadspatial range on the subducting interface may be contributing to the slip. Nonetheless, the subduction zone at largescale follows similar patterns in energy; the machine learning algorithm identified these patterns, enabling timingestimations.In the laboratory, these signals carry fundamental information regarding the frictional state of the system. Thefact that a similar behavior can be observed in Cascadia, with energy patterns matching the 14-month slow slipcycle, may also provide indirect information regarding the evolution of friction at much larger scale. The originof the seismic signal is presumed to be due to fault gouge grain-to-grain and grain-to-block interactions. Effort isongoing to explain its origin [
33, 34 ]. In the Cascadia subduction zone, we posit that part of this energy is emittedfrom a large number of asperities located on the fault interface. To refine the analysis and focus exclusively ontremor, we are training models to recognize tremor and estimate cycles in tremor energy [ ], which may ultimatelyhelp to improve the slip timings estimations. Conclusion
We find that seismic energy in the Cascadia subduction zone follows specific patterns throughout the slow slipcycle, with energy starting to build up steadily before failure. This initial increase in energy is not captured in thecataloged tremor - suggesting that slow slip events often begin with a silent phase. Slow slip has been observedprior to large earthquakes in some instances, but not systematically. If indeed slow earthquakes have identifiableprecursory signals, estimating their timing may ultimately help improve hazard assessment for destructive earth-quakes. 8 cknowledgements
We thank C. Marone, C. Bolton and J. Rivi`ere for the laboratory data. The seismic data used were obtained fromthe Plate Boundary Observatory [ ]. All the data is publicly available. This work was supported by DOE Officeof Science (Geoscience Program, grant 89233218CNA000001) and Institutional Support (LDRD) at Los Alamos.CH work was also done thanks to a joint research laboratory effort in the framework of the CEA-ENS Yves RocardLRC (France), and thanks to the European Research Council (ERC) under the European Union’s Horizon 2020research and innovation program (Geo-4D project, grant agreement 758210). Author contributions statement
C.H. and B.R.L. conducted the machine learning analysis. P.A.J. supervised the project. All authors contributed towriting the manuscript.
References
1. K. Obara, Nonvolcanic deep tremor associated with subduction in southwest japan.
Science , 1679–1681(2002).2. K. Obara, A. Kato, Connecting slow earthquakes to huge earthquakes.
Science , 253–257 (2016).3. D. M. Veedu, S. Barbot, The Parkfield tremors reveal slow and fast ruptures on the same asperity.
Nature ,361 (2016).4. Z. Peng, J. Gomberg, An integrated perspective of the continuum between earthquakes and slow-slip phenom-ena.
Nature Geoscience , 599 (2010).5. G. C. Beroza, S. Ide, Slow earthquakes and nonvolcanic tremor. Annual Review of Earth and PlanetarySciences , 271-296 (2011). 9. J. C. . Gomberg, B. W. Group, Slow-slip phenomena in cascadia from 2007 and beyond: A review. GSABulletin , 963–978 (2010).7. S. D. L. Rubinstein, J.L., W. L. Ellsworth,
Non-volcanic Tremor: A Window into the Roots of Fault Zones (Springer Science+Business Media B.V., 2010).8. B. Rousset, R. B¨urgmann, M. Campillo, Slow slip events in the roots of the san andreas fault.
Science Advances (2019).9. C. H. Scholz, The Mechanics of Earthquakes and Faulting (Cambridge University Press, 2019), pp. i–i, thirdedn.10. X. Gao, K. Wang, Rheological separation of the megathrust seismogenic zone and episodic tremor and slip.
Nature , 416 (2017).11. W. B. Frank, Slow slip hidden in the noise: The intermittence of tectonic release.
Geophysical ResearchLetters , 10,125–10,133 (2016).12. Y. Ito, R. Hino, M. Kido, H. Fujimoto, Y. Osada, D. Inazu, Y. Ohta, T. Iinuma, M. Ohzono, S. Miura,M. Mishina, K. Suzuki, T. Tsuji, J. Ashi, Episodic slow slip events in the japan subduction zone before the2011 tohoku-oki earthquake. Tectonophysics , 14 - 26 (2013).13. S. Ruiz, M. Metois, A. Fuenzalida, J. Ruiz, F. Leyton, R. Grandin, C. Vigny, R. Madariaga, J. Campos, Intenseforeshocks and a slow slip event preceded the 2014 iquique mw 8.1 earthquake.
Science , 1165–1169(2014).14. V. M. Kuna, J. L. N´abˇelek, J. Braunmiller, Mode of slip and crust–mantle interaction at oceanic transformfaults.
Nature Geoscience , 138 (2019).15. C. Hulbert, B. Rouet-Leduc, P. A. Johnson, C. X. Ren, J. Rivi`ere, D. C. Bolton, C. Marone, Similarity of fastand slow earthquakes illuminated by machine learning. Nature Geoscience , 69 (2019).106. C. Marone, Laboratory-derived friction laws and their application to seismic faulting. Annual Review of Earthand Planetary Sciences , 643–696 (1998).17. B. M. Kaproth, C. Marone, Slow earthquakes, preseismic velocity changes, and the origin of slow frictionalstick-slip. Science , 1229–1232 (2013).18. M. Scuderi, C. Marone, E. Tinti, G. Di Stefano, C. Collettini, Precursory changes in seismic velocity for thespectrum of earthquake failure modes.
Nature geoscience , 695 (2016).19. B. Rouet-Leduc, C. Hulbert, P. A. Johnson, Continuous chatter of the cascadia subduction zone revealed bymachine learning. Nature Geoscience (2019).20. G. Rogers, H. Dragert, Episodic tremor and slip on the cascadia subduction zone: The chatter of silent slip.
Science , 1942–1943 (2003).21. A. G. Wech, K. C. Creager, T. I. Melbourne, Seismic and geodetic constraints on cascadia slow slip.
Journalof Geophysical Research: Solid Earth (2009).22. A. G. Wech, N. M. Bartlow, Slip rate and tremor genesis in Cascadia.
Geophysical Research Letters ,392-398 (2014).23. J. C. Hawthorne, M. G. Bostock, A. A. Royer, A. M. Thomas, Variations in slow slip moment rate associatedwith rapid tremor reversals in Cascadia. Geochemistry, Geophysics, Geosystems , 4899–4919 (2016).24. H. Kao, S.-J. Shan, H. Dragert, G. Rogers, J. Cassidy, K. Ramachandran, A wide depth distribution of seismictremors along the northern Cascadia margin. Nature , 841-844 (2005).25. T. Perol, M. Gharbi, M. Denolle, Convolutional neural network for earthquake detection and location.
ScienceAdvances (2018).26. Z. E. Ross, M.-A. Meier, E. Hauksson, P wave arrival picking and first-motion polarity determination withdeep learning. Journal of Geophysical Research: Solid Earth , 5120-5129 (2018).117. I. W. McBrearty, A. A. Delorey, P. A. Johnson, Pairwise association of seismic arrivals with convolutionalneural networks.
Seismological Research Letters , 503–509 (2019).28. Z. E. Ross, Y. Yue, M.-A. Meier, E. Hauksson, T. H. Heaton, Phaselink: A deep learning approach to seismicphase association. Journal of Geophysical Research: Solid Earth , 856-869 (2019).29. B. K. Holtzman, A. Pat´e, J. Paisley, F. Waldhauser, D. Repetto, Machine learning reveals cyclic changes inseismic source spectra in geysers geothermal field.
Science Advances (2018).30. P. G. Silver, Y. Bock, D. C. Agnew, T. Henyey, A. T. Linde, T. V. McEvilly, J.-B. Minster, B. A. Romanowicz,I. Sachs, R. B. Smith, et al. , A plate boundary observatory. Iris Newsletter , 3 (1999).31. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer,R. Weiss, V. Dubourg, et al. , Scikit-learn: Machine learning in python. Journal of Machine Learning Research , 2825–2830 (2011).32. B. Rouet-Leduc, C. Hulbert, D. C. Bolton, C. X. Ren, J. Riviere, C. Marone, R. A. Guyer, P. A. Johnson,Estimating fault friction from seismic signals in the laboratory. Geophysical Research Letters , 1321–1329(2018).33. K. Gao, R. A. Guyer, E. Rougier, C. X. Ren, P. A. Johnson, From force chains to acoustic emission. PhysicalReview Letters (2019).34. C. X. Ren, O. Dorostkar, B. Rouet-Leduc, C. Hulbert, D. Strebel, R. A. Guyer, P. A. Johnson, J. Carmeliet,Machine learning reveals the state of intermittent frictional dynamics in a sheared granular fault.
GeophysicalResearch Letters , 7395-7403 (2019).35. B. Rouet-Leduc, C. Hulbert, I. McBrearty, P. A. Johnson, Deep learning slow earthquakes in cascadia anduniversal frictional characteristics. arXivarXiv