Studies of seismic deconvolution andlow-frequency earthquakesbyAlexandra Ame´lie RoyerB. Sc. Universite´ Joseph Fourier, 2007M. Sc. Universite´ Joseph Fourier, 2009A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Geophysics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2014c© Alexandra Ame´lie Royer, 2014iiAbstractDeconvolution of seismic data is an important component of signal processing thataims to remove the seismic source from seismograms, thereby isolating the Green’sfunction. By considering seismograms of multiple earthquakes from similar locationsrecorded at a given station and that therefore share the same Green’s function, we inves-tigate a system of equations where the unknowns are the sources and source durations.Our solution is derived using direct linear inversion to recover the sources and Newton’smethod to recover source durations. For the short seismogram durations considered,we are able to recover source time functions for noise levels at 1% of the direct P -waveamplitude. However, the nonlinearity of the problem renders the system expensive tosolve and sensitive to noise; therefore consideration is limited to short seismograms withhigh signal-to-noise ratio (SNR). When SNR levels are low, but a large multiplicity ofseismograms representing a common source-receiver path are available, we can apply adifferent deconvolution approach to recover the Green’s function. In an application totectonic tremor in northern Cascadia, we implement an iterative blind deconvolutionmethod that involves correlation, threshold detection and stacking of 1000’s of low fre-quency earthquakes (LFEs) that form part of tremor to generate templates that canbe considered as empirical Green’s functions. We exploit this identification to computehypocentres and moment tensors. LFE hypocentres follow the general epicentral dis-tribution of tremor and occur along tightly defined surfaces in depth. The majorityof mechanisms are consistent with shallow thrusting in the direction of plate motion.We analyze the influence of ocean tides on the triggering of LFEs and find a spatiallyvariable sensitivity to tidally induced up-dip shear stress (UDSS), suggesting that tidalsensitivity must partially depend on laterally heterogeneous physical properties. Themajority of LFEs fail during positive and increasing UDSS, consistent with combinedcontributions from background slow slip and from tides acting directly on LFEs. Weidentified rapid tremor reversals in southern Vancouver Island with higher sensitivityto UDSS than the main front and which at least partially explains an observed increasein LFE sensitivity to UDSS with time.iiiPrefaceAll of the work presented herein was carried out at The University of BritishColumbia by me, with advice, supervision and editorial authorship by my supervi-sor, Michael Bostock. My thesis committee (Catherine Johnson, Eldad Haber andAndrew Calvert) provided comments and suggestions during committee meetings andin personal discussions.Chapter 2 has been published as Royer, A. A., M. G. Bostock and E. Haber (2012),Blind deconvolution of seismograms regularized via minimum support, Inverse problems,28, 127,010, doi:10.1088/0266-5611/28/12/125010. It is reproduced unmodified, ex-cept for formatting changes and minor corrections for clarity. I was the lead author andwas responsible for the analysis of the data and composition of the paper. Michael Bo-stock and Eldad Haber provided feedback, advice and editorial supervision throughoutthe development of this paper.Chapter 3 has been published as Royer, A. A. and M. G. Bostock (2013), A compar-ative study of low frequency earthquake templates in northern Cascadia, Earth Planet.Sci. Lett., doi:10.1016/j.epsl.2013.08.040. It is reproduced unmodified, except forformatting changes and minor corrections for clarity. I was the lead author and wasresponsible for the analysis of the data and composition of the paper. Michael Bostockprovided feedback, advice and editorial supervision throughout the development of thispaper.Chapter 4 has been submitted in July 2014, entitled “Tidal modulation of lowfrequency earthquakes and triggering of secondary events in northern Cascadia”. I amthe lead author and I am responsible for the analysis of the data and composition ofthe paper. Co-author Amanda Thomas from Stanford University, California, provideddata, technical advice and editorial supervision. Michael Bostock provided feedback,advice and editorial supervision throughout the development of this paper.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Overview of deconvolution approaches . . . . . . . . . . . . . . . . . . 21.3 Tectonic context of the northern Cascadia subduction zone . . . . . . 41.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.2 Tremor, slow slip and low frequency earthquakes . . . . . . . . . 71.3.3 Earth tides and low frequency earthquakes . . . . . . . . . . . . 91.3.4 Rapid tremor reversals . . . . . . . . . . . . . . . . . . . . . . . . 102 Blind deconvolution of seismograms regularized via minimumsupport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3.1 Source duration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3.2 Introduction of source durations as unknowns . . . . . . . . . . . 162.3.3 Combined system . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.4 Separation of variables . . . . . . . . . . . . . . . . . . . . . . . . 182.3.5 Stopping criterion . . . . . . . . . . . . . . . . . . . . . . . . . . 212.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283 A comparative study of low frequency earthquake templatesin northern Cascadia . . . . . . . . . . . . . . . . . . . . . . . . . 353.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.2.1 Washington state . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.2.2 Vancouver Island . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3 Data Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3.1 Network autocorrelation . . . . . . . . . . . . . . . . . . . . . . . 393.3.2 Waveform-correlation cluster analysis . . . . . . . . . . . . . . . 393.3.3 Network cross-correlation . . . . . . . . . . . . . . . . . . . . . . 403.4 LFE templates as empirical Green’s functions . . . . . . . . . . . . . . 403.5 LFE template locations . . . . . . . . . . . . . . . . . . . . . . . . . . 443.6 LFE template moment tensors . . . . . . . . . . . . . . . . . . . . . . 47v3.7 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . . 494 Tidal modulation of low frequency earthquakes and triggeringof secondary events in northern Cascadia . . . . . . . . . . . . . 584.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.3.1 Tidal model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.3.2 LFE correlation with tidal stress . . . . . . . . . . . . . . . . . . 644.3.3 Time dependent tidal sensitivity . . . . . . . . . . . . . . . . . . 714.3.4 Tidal sensitivity of rapid tremor reversals vs the main front . . . 734.3.5 Diurnal variation of LFE occurrence . . . . . . . . . . . . . . . . 774.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.4.1 LFE correlation with tidal stress . . . . . . . . . . . . . . . . . . 784.4.2 Time dependent tidal sensitivity . . . . . . . . . . . . . . . . . . 794.4.3 Tidal sensitivity of RTRs vs the main front . . . . . . . . . . . . 794.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.5.1 Spatial distribution of tidal sensitivity . . . . . . . . . . . . . . . 804.5.2 The role of secondary events . . . . . . . . . . . . . . . . . . . . 834.5.3 Tidal sensitivity and LFE recurrence intervals . . . . . . . . . . . 854.5.4 Phase of LFE asperity failure relative to tidal stress . . . . . . . 874.5.5 Evolution of tidal sensitivity through slow slip events. . . . . . . 894.5.6 24 hour periodicity of LFE detections . . . . . . . . . . . . . . . 894.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.1 Key contributions and limitations . . . . . . . . . . . . . . . . . . . . 925.1.1 Blind deconvolution . . . . . . . . . . . . . . . . . . . . . . . . . 925.1.2 Comparative study of LFE templates in northern Cascadia . . . 935.1.3 Tidal modulation of LFEs . . . . . . . . . . . . . . . . . . . . . . 945.2 Future research directions . . . . . . . . . . . . . . . . . . . . . . . . . 955.2.1 Multigrid strategies in deconvolution . . . . . . . . . . . . . . . . 955.2.2 Green’s function recovery applications . . . . . . . . . . . . . . . 955.2.3 High resolution detection and location of LFEs using cross-stationapproaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.2.4 Structural studies using LFE templates . . . . . . . . . . . . . . 965.2.5 Tidal modulation of RTRs and streaks and implication for plateboundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111A Supplementary material for Chapter 3 . . . . . . . . . . . . . . . 112B Supplementary material for Chapter 4 . . . . . . . . . . . . . . . 116viList of Tables2.1 Inversion parameters employed in numerical experiments. . . . . . . . . . . 262.2 Schematic representation of the 3 classes of experiment. . . . . . . . . . . . 274.1 RTRs documentation for 2003, 2004, 2005, 2008. . . . . . . . . . . . . . . . 754.2 RTRs documentation for 2010, 2011, 2012 . . . . . . . . . . . . . . . . . . . 76A.1 Parameters employed in hypoDD inversion . . . . . . . . . . . . . . . . . . 113B.1 UDSS and FNS Nex values for RTRs and non-RTR events . . . . . . . . . . 120viiList of Figures1.1 Map of Cascadia subduction zone . . . . . . . . . . . . . . . . . . . . . . . . 61.2 Picture of the subducting Juan de Fuca plate . . . . . . . . . . . . . . . . . 71.3 2012 Tremor episode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1 Problem configuration in the context of global seismology . . . . . . . . . . 152.2 Weighting function w(t, τ) for different values of a and τ = 4 . . . . . . . . 172.3 P- and S-component of Green’s function employed in numerical experiments 252.4 P- and S-seismograms generated with random Gaussian source . . . . . . . 252.5 P- and S-seismograms generated with band-limited source . . . . . . . . . . 262.6 Inversion using random Gaussian sources and P-components . . . . . . . . . 292.7 Inversion using random Gaussian sources and S-components . . . . . . . . . 302.8 Inversion using random Gaussian sources and P- and S-components . . . . 312.9 Inversion using band-limited sources and P-components . . . . . . . . . . . 322.10 Inversion using band-limited sources and S-components . . . . . . . . . . . 332.11 Inversion using band-limited sources and P- and S-components . . . . . . . 343.1 Distribution of stations used in the study of LFE templates . . . . . . . . . 383.2 North, east and vertical components of 3 LFE templates . . . . . . . . . . . 413.3 Three-component LFE template waveforms . . . . . . . . . . . . . . . . . . 433.4 Maps of LFE locations computed using the Hyp2000 software . . . . . . . . 453.5 Depth profiles of seismicity in Vancouver Island and Washington state . . . 473.6 Example of waveform matching in moment tensor inversion . . . . . . . . . 533.7 Map of double couple mechanisms determined from moment tensor inversion 543.8 Map of double couple mechanisms for NVI subarray . . . . . . . . . . . . . 553.9 Map of double couple mechanisms for SW subarray . . . . . . . . . . . . . . 563.10 Map of double couple mechanisms for NW subarray . . . . . . . . . . . . . 574.1 Map of LFE family locations . . . . . . . . . . . . . . . . . . . . . . . . . . 614.2 2007 tremor episode time series of tidally-induced normal and shear stresses 634.3 Map of tidal stress amplitudes . . . . . . . . . . . . . . . . . . . . . . . . . 654.4 Map of Nex values relative to tidal stresses and stress rates . . . . . . . . . 674.5 Nex values relative to tidal stresses and stress rates plotted against each other 684.6 Correlation between Nex values and mean envelope amplitude . . . . . . . . 694.7 LFE rate plots as a function of the magnitude of the 4 tidal components . . 704.8 Phase histograms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.9 Nex UDSS histograms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.10 Along-strike profile and time-distance plots for 2003 and 2012 tremor episodes 744.11 LFE rate plots as a function of hour of day (Pacific time) . . . . . . . . . . 814.12 Cumulative Distribution functions . . . . . . . . . . . . . . . . . . . . . . . 86A.1 Map of LFE locations from hypoDD inversion . . . . . . . . . . . . . . . . . 113A.2 Depth profiles A-A’ to J-J’ of seismicity in northern Cascadia . . . . . . . . 114A.3 Condition number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114A.4 Histograms of strike, dip and rake distribution . . . . . . . . . . . . . . . . 115B.1 hypothesis test for spurious correlation . . . . . . . . . . . . . . . . . . . . . 121viiiList of AbbrevationsCNSN Canadian National Seismograph NetworkdFNS Fault Normal Stress ratedUDSS Up Dip Shear Stress rateETS Episodic Tremor and SlipFNS Fault Normal StressGCV Generalized Cross ValidationGPS Global Positionning SystemLFE Low Frequency EarthquakeLVZ Low Velocity ZoneNVI Northern Vancouver IslandNW Northern WashingtonPBO Plate Boundary ObservatoryPC Principal ComponentPCA Principal Component AnalysisPNSN Pacific Northwest Seismic NetworkRTR Rapid Tremor ReversalSNR Signal to Noise RatioSSE Slow Slip EventSVI Southern Vancouver IslandSW Southern WashingtonTA Transportable USArrayUDSS Up Dip Shear StressixList of SymbolsChapter 1c(t) cross-correlation functionF {δ(t)} filtered Dirac delta functiong(t) Green’s functiongˆ(t) Green’s function estimatesi(t) ith source signature time seriesui(t) ith seismogram time seriesVP elastic P -wave velocityVS elastic S-wave velocityδ(t) Dirac delta functionΦˆg(t) estimate of the true autocorrelation of the Green’s function g(t)⋆ convolution product⊗ cross-correlation productChapter 2a steepness parameterA matrix comprising a combination of Ui, Uj and zero blocksB matrix containing the derivatives ∂wi(t,τ )∂τie scaling constraintg(t) Green’s functionH(t) Heaviside functionH Hessian matrixL Lagrange functionp slownesssi(t) ith source signature time serieui(t) ith seismogram time serieUi Toeplitz matrice containing the seismogram ui(t)xVP elastic P -wave velocityVS elastic S-wave velocitywi(t, τ ) windowing functionW diagonal weighting matrix containing wi(t, τ )x source function vectory “pre-weighted” source function vectorY diagonal matrix containing yα constantβ Tikhonov parameterδ(t) Dirac delta functionλ Lagrange multiplier∇L gradient of Langrange functionρ densityΣ covariance matrixτ durationτ 0 initial duration⋆ convolution productChapter 3g(t) Green’s functiongˆ(t) Green’s function estimatesi(t) ith source signature time seriesui(t) ith seismogram time seriesΦˆg(t) estimate of the true autocorrelation of the Green’s function g(t)⊗ cross-correlation productChapter 4a friction coefficientD total amount of slip accumulated over a tremor episodexiDc characteristic length scaleerr,eθθ,eλλ vertical and horizontal strainsLm length of main slip frontLpatch LFE patch dimensionLs along-dip distanceNex degree of correlation of LFE occurrence with tidal stressR sample correlation coefficientT tidal periodTθ period that defines time scale for the evolution of state variableV creep rate of the fault surrounding an LFE patchVL long-term displacement averaged creep rateVss steady state velocity during SSE∆V velocity perturbationVP elastic P -wave velocityVS elastic S-wave velocity∆τsd stress drop∆τss tidal shear stress amplitudeµ shear modulusν Poisson’s ratioσe effective normal stressτ shear resistanceτ0 shear resistance at creep rate VLθ phase advance wrt peak UDSSChapter 5VP elastic P -wave velocityVS elastic S-wave velocityy “pre-weighted” source function∆τ positive change of source duration with iterationsxiiAcknowledgementsFirst and foremost, I would like to express my profound gratitude and utmost re-spect to my research advisor Michael Bostock. I benefited from his exemplary guidance,monitoring and constant encouragements throughout the course of this thesis. I willnot forget the opportunities Michael gave me to travel and share my work with the ”bigshots” of seismology.I would like to thank the co-authors of my publications Eldad Haber and AmandaThomas. I am very grateful for their contributions to the development of my research.Their participation, discussions and our meetings have been of greatest benefit.My gratitude also goes to the members of my PhD committee, Eldad Haber, Cather-ine Johnson and Andrew Calvert from Simon Fraser University, for discussion and guid-ance during committee meetings.I gratefully acknowledge the helpful comments of Nick Beeler, Roland Burgmann,Jessica Hawthorne, Allan Rubin and David Shelly on chapter 4.I would also like to show my appreciation to the feedback of Professor Simon Pea-cock on the development of my findings and extend my thanks to the computer techni-cians of the EOAS department for their help in running my code.I thank my colleagues Ben Postlethwaite, Gian Matharu, Genevie`ve Savard andBrendan Smithyman for their participation in many ways.The completion of this Ph.D is the culmination of 10 years of university studiesin France and Canada. I will always remember my professors Alexandre Fournier,Andrea Walsperdorf, Etienne Jaillard, Jean-Robert Grasso and Helle Pedersen fromthe Universite´ Joseph Fourier for their encouragements and support to reach this levelof understanding.I would like to give my loving thanks to my family and particularly to my aunt anduncle Annick and Pierre Spierer and to my grandparents Jacques and Jacqueline Royer,for their support and love. I am very grateful for their encouragements to find my ownpath far from them and honored to join the exclusive Ph.D club in our family! Thankyou to my family of heart, Babeth, Roland, Coraline, Cle´ment, Fabien, Lisa, Marilou,Alice, Nathalie and Mustapha. They have been there for me and with me every dayand I am very grateful for having them in my life.Lastly, I thank Ralf Hansen. His support in the good moments, but also in the moredifficult ones, has been a precious help. His presence became essential, particularly thispast year, for the final climb to the summit. Thank you.Je de´die cette the`se a` mon de´funt pe`re Pierre Royer. En gravant les plus hautssommets enneige´s, tu m’as inculque´ le gouˆt de l’effort et de l’ambition. Tu m’as renducurieuse des forces de la nature et aujourd’hui, je ne serais pas la` sans toi.Je de´die cette the`se a` ma de´funte me`re Martine Royer. Comme une e´toile dans le ciel,tu me guides et me soutiens. Tu m’as appris a` combattre l’adversite´, a` croire en moiet a` re´aliser mes reˆves. Aujourd’hui, je ne serais pas la` sans toi.Merci a` tous les deux. Merci pour tout.“You can’t stop the waves,but you can learn to surf”Joseph GoldsteinChapter 1. Introduction 1Chapter 1Introduction1.1 Research objectivesI present in this thesis three different research topics that can be read indepen-dently of each other. However, it is possible to relate them through the general conceptof deconvolution to estimate a Green’s function. The first topic relates to the blinddeconvolution of seismograms regularized via minimum support. Through an alge-braic analysis, I have built a mathematical model to recover the earthquake sourcesignature in terms of amplitude and source duration from recorded seismograms. Thesecond topic exploits a semi-empirical deconvolution approach involving iterative corre-lation–detection–stacking to generate low frequency earthquake (LFE) templates fromtremors recorded in northern Vancouver Island and Washington state, thereby extend-ing a previous LFE catalogue for southern Vancouver Island [Bostock et al., 2012]through much of northern Cascadia. I provide empirical and semi-analytical argu-ments justifying the identification of LFE templates with “Green’s function” sectionscorresponding to moment tensor point sources. I proceed to compare the distributionand excitation of LFE sources across northern Cascadia as revealed by LFE templatesand discuss implications for plate boundary structure and LFE genesis. The third topicof this thesis relates to the influence of the solid Earth and ocean tides on a catalog ofLFEs as derived from LFE templates for southern Vancouver Island and Washingtonstate. The spatial and temporal sensitivity of LFEs provides insight into the sourceChapter 1. Introduction 2region (i.e plate boundary) properties such as friction and pore fluid pressure, but alsoon the mechanism of displacement and conditions of triggering.1.2 Overview of deconvolution approachesIn seismology, the Green’s function describes the impulse response of the Earththrough which seismic waves travel. It describes the effects of reflections and refrac-tions of seismic waves at interfaces between rocks of different properties, the phasedelays/advances and geometrical spreading due to lateral changes in velocity structure,and the influence of attenuation due to anelasticity. Knowledge of the Green’s functioncan aid significantly in imaging the internal structure of the Earth. Deconvolution ofseismic data is an important component of signal processing that aims to remove theseismic source from seismograms, thereby isolating the Green’s function.The two deconvolution approaches I employ in this thesis start from the same math-ematical description in which a seismic record u(t) is defined as the convolution of aseismic source s(t) with a Green’s’ function g(t), expressed in the time domain as:u(t) = s(t) ⋆ g(t), (1.1)where ⋆ denotes the convolution product. Blind deconvolution, the main topic of chapter2, refers to the isolation of either s(t) or g(t) without prior knowledge of the other. Inthis approach, I consider seismograms of multiple earthquakes from similar locationsrecorded at a given station and that therefore share the same Green’s function. We canwrite a linear relation in the time domain:uj(t) ⋆ si(t)− ui(t) ⋆ sj(t)= g(t) ⋆ si(t) ⋆ sj(t)− g(t) ⋆ sj(t) ⋆ si(t)= 0,(1.2)where ui(t) is the seismogram for the ith source and sj(t) is the jth unknown source.From two or more seismograms, we obtain a system where the unknowns are the sources.Our approach is to identify the Green’s function as the common convolutional elementwithin a set of seismograms, but the approach might equally well apply to a commonChapter 1. Introduction 3source, as implied by the symmetry of s(t) and g(t) in 1.1. Although the topic ofblind deconvolution has been addressed by many authors, the originality of this firstapproach is that it requires only multiple seismograms recording the same sources (orreceivers), and regularization by minimum support. However, the nonlinearity of thisanalysis renders the system expensive to solve and sensitive to noise, and therefore itcan be applied to only short duration seismograms with high SNR.Non-volcanic or tectonic tremors are weak signals with low SNR originating frommajor convergent or transform plate boundaries. It has been shown that this tremorcomprises swarms of many small, repeating earthquakes known as LFEs [Shelly et al.,2007]. Since the SNR levels are low we must rely on a different deconvolution approachthat exploits the multiplicity of repeating sources. The second blind deconvolutionapproach described in this thesis employs an iterative strategy that involves correlation,threshold detection and stacking to isolate the Green’s function. It correlates a smallwindowed section of a Green’s function estimate (e.g. a single LFE) at a number ofchannels with the corresponding tremor time series. Times for which the summedcorrelation coefficient exceeds a specified threshold are considered further detections ofLFEs representing the same Green’s function. An individual tremor episode typicallylasting 2-3 weeks may contain more than 1000 LFE detections emanating from the samelocation. These detections can be aligned and stacked to produce a high SNR template.The template can be used as an updated Green’s function estimate and the same processis then iterated until there is no further improvement in SNR. The interpretation ofthe template as an empirical Green’s function can be justified as follows. Consider thecross correlation c(t) of a Green’s function estimate gˆ(t) with a seismogram u(t). Wecan write, using 1.1:c(t) = gˆ(t)⊗ u(t) = gˆ(t)⊗ g(t) ⋆ s(t) = Φˆg(t) ⋆ s(t), (1.3)where ⊗ denotes cross-correlation. Φˆg(t) is an estimate of the true autocorrelation ofg(t) and so will peak near lag t=0. The timing of the maximum of c(t) will then occurat or near the time t = tmax at which s(t) possesses its maximum. Times for which c(t)exceeds a certain threshold are considered as detections and the stack of correspondingChapter 1. Introduction 4seismograms ui(t), each shifted by the corresponding detection time tmaxi , can be writtenas:∑iui(t− tmaxi ) = g(t) ⋆∑isi(t− tmaxi ). (1.4)As the number N of detections increases, the sum of shifted band-limited source func-tions si(t − tmax) will tend toward a scaled, filtered delta function (since LFEs andtremor are inherently band limited to 1-10Hz),limN→∞N∑isi(t− tmaxi ) → F {δ(t)} , (1.5)and a scaled, band-limited version of g(t) can be approximated by summing the shiftedui(t). The detection of LFEs that repeat 1000’s of times from the same location overperiods of a few years allow us to assemble fully 3-D empirical Green’s functions, incontrast with previous studies (stacking of phase-normalized long-period body waves[Shearer , 1991] and broadband teleseismic P−waves [Kumar et al., 2010]) for whichthe relative infrequency of regular seismicity limits Green’s function retrieval to 1-Destimates.1.3 Tectonic context of the northern Cascadia subductionzone1.3.1 OverviewThe Cascadia subduction zone is located along the western coast of North America,where the oceanic Juan de Fuca, Explorer and Gorda plates are being subducted be-neath the continental North American plate. It stretches over 1000 km from northernVancouver Island to northern California. For the purpose of this thesis, we will focuson the northern part of Cascadia from northern Vancouver Island to southern Wash-ington state (Figure 1.1). The Juan de Fuca plate subducts at a rate of ∼40 mm/yearbeneath the North American plate and, due to its relatively young age (≤ 11Myr), itis associated with warmer temperatures than most other subduction zones.Chapter 1. Introduction 5The analysis of seismic waves affords several different means for studying the in-ternal structure of the subduction zone and the geological processes that are activelyreshaping it. Deep slab structure (> 100 km depth) has been primarily imaged us-ing teleseismic P -wave traveltime tomography. Below Washington and western BritishColumbia, the slab is characterized by a quasi-planar, high P -velocity slab evident toat least 400 km depth [Audet et al., 2008]. The slab dip varies from about 50◦ to 60◦to 45◦ below southern Vancouver island, northern Washington and southern Washing-ton, respectively [Bostock and Vandecar , 1995, Burdick et al., 2008, Chu et al., 2012,Mercier et al., 2009, Michaelson and Weaver , 1986, Obrebski et al., 2011, Rasmussenand Humphreys, 1988]. It dips to the northeast below Vancouver Island and to the eastbelow Washington state. Slab structure at shallower levels (<100 km depth) beneathsouthern British Columbia and Washington state has been investigated using seismicreflection and receiver functions analysis, that reveal a landward-dipping, low-velocityzone (LVZ), first observed by Langston [1977, 1981] and interpreted as the subductingoceanic crust. It coincides with a reflective and conductive layer originally identifiedin seismic reflection studies [Clowes et al., 1987, 1990, Green et al., 1986] and laterin magnetotelluric surveys [Kurtz et al., 1990, Soyer and Unsworth, 2006], known asthe “E”-layer on the southern Vancouver Island Lithoprobe transect. The LVZ hasbeen most recently interpreted to be upper oceanic crust [Bostock , 2013, Hansen et al.,2012]. A schematic cartoon illustrating important tectonic elements of the Cascadiasubduction zone is provided in Figure 1.2.In this study I will use 2 different models for the Juan de Fuca subducting plate,one constructed by Audet et al. [2010] and the other by McCrory et al. [2012], to anal-yse LFE signals. Audet et al. [2010] identified the plate boundary with the top ofthe LVZ using scattered teleseismic P arrivals and by assuming fixed P -velocity of theoverriding plate with locally variable VP/VS ratio. In contrast, McCrory et al. [2012]used intraplate earthquake locations and regional seismic velocity profiles to synthesizeinformation of the depth structure. The authors identified the top of the Juan de Fucaslab near the upper surface of intraplate seismicity and weighted these seismicity con-straints more highly than structural information derived from velocity profiles in areaswhere both sources of information were available. The two slab models are plotted inChapter 1. Introduction 6Figure 1.1.−130˚ −129˚ −128˚ −127˚ −126˚ −125˚ −124˚ −123˚ −122˚39˚40˚41˚42˚43˚44˚45˚46˚47˚48˚49˚50˚51˚GORDAPLATEPACIFICPLATEJUAN DE FUCAPLATEEXPLORER PLATEExplorer ridgeJuan de Fuca ridgeBlanco faultNootkafaultCascadia trenchNORTHERN CASCADIA SUBDUCTION ZONE40 mm/yearMendocino faultFigure 1.1: Map of Cascadia subduction zone. Northern Cascadia subduction zoneis delimited by a dashed square. Cyan and orange lines indicate the 20, 30 and 40 kmdepth contours to the top of the subducting Juan de Fuca plate modeled by Audetet al. [2010] and McCrory et al. [2012], respectively.Chapter 1. Introduction 7Upper oceanic crust = LVZOceanic mantle Continental Moho40kmLOCKED ZONESLOW TRANSIENT SLIP ZONEarc volcanoMainlandGeorgia StraitVancouver Isl.OceanOverpressuredoceanic crustSerpentinized mantle wedgePermeableplate boundaryFigure 1.2: Schematic cartoon of subducting plate in northern Cascadia, modifiedfrom Audet et al. [2009]1.3.2 Tremor, slow slip and low frequency earthquakesThe Juan de Fuca plate does not subduct continuously. At shallow depths (< 25km depth), thermal and deformation studies [Dragert et al., 1994, Wang et al., 1994]indicate a locked zone approximately 60 km wide in which strain energy is accumu-lated without rupture over a period of few hundred years (Figure 1.2). Using GlobalPositioning System (GPS) monitoring, slow slip events equivalent to a M ∼6-7 mag-nitude earthquake have been detected farther downdip (depths of 25-40 km, referredas “slow transient slip zone” in Figure 1.2) beneath Vancouver Island and Washing-ton state. These slow slip events occur over a period of a few weeks [Dragert et al.,2001] and repeat at intervals of 13 to 16 months [Miller et al., 2002]. First thought tobe seismically silent, the slip is accompanied by pulsating, tremorlike seismic signals[Rogers and Dragert , 2003], that is, the tremor signals described in section 1.2. Thesetectonic tremors are weak, semi-continuous and characterized by a frequency spectrumthat is peaked between 1 and 10 Hz. Figure 1.3 show the 2012 tremor episode locationsin northern Cascadia. Tremor locations in Cascadia can been viewed on an interactivewebpage (http : //tunk.ess.washington.edu/map display), built by Aaron Wech [Wechand Creager , 2008].Improvements in seismic monitoring networks and in signal processing techniquesled to the discovery that LFEs form part of tectonic tremor [Shelly et al., 2007]. De-spite their low magnitudes and limited bandwidth, these signals can be detected at lowChapter 1. Introduction 8Figure 1.3: 2012 episodic tremor and slip event (picture from Wech’s interactivewebpage displaying automatic tremor locations, referred in text).SNR thresholds and with high temporal precision using powerful network correlationtechniques [Brown et al., 2008, Gibbons and Ringdal , 2006, Shelly et al., 2006] as sum-marized in section 1.2 and further developed in chapter 3.The fact that LFE templates can be considered as empirical Green’s functions isuseful in the identification of arrivals for location and for waveform matching as requiredin moment tensor inversion. Previous spatio-temporal studies of LFEs in Shikoku andCascadia suggest that tremor and LFEs are the result of slow slip at the plate interface[Bostock et al., 2012, Shelly et al., 2007, Wech and Creager , 2007] and the combina-tion of tremor and slow slip is often referred to as Episodic Tremor and Slip (ETS). Acomparative study of LFEs across northern Cascadia is developed in chapter 3 of thisthesis.Chapter 1. Introduction 91.3.3 Earth tides and low frequency earthquakesThe combined action of the Sun and the Moon on the Earth produces periodic tidalstress variations of the order of 4 kPa in the Earth’s crust to depths of 10’s of km. Theamplitudes of the tidal stress variations are up to three orders of magnitude smallerthan the stress drop produced by earthquakes. However, tidal stress rate is generallyhigher than the tectonic stress accumulation on a fault. Therefore, tidal stresses couldtrigger earthquakes on a weak fault close to the critical rupture level. The correlationbetween tides and earthquake occurrence has been investigated for the past century andit has transpired that most studies have not been able to highlight a correlation betweenearthquakes and tides [Knopoff , 1964, Rydelek et al., 1992, Schuster , 1897, Shudde andBarr , 1977, Simpson, 1967, Vidale et al., 1998]. However, more recent studies on largenumbers of earthquakes with input from laboratory experiments [Beeler and Lockner ,2003] show the existence of a weak correlation between the two phenomena [Cadicheanuet al., 2007, Enescu and Enescu, 1999, Stavinschi and Souchay , 2003, Tanaka et al.,2002, 2006].In regions where tectonic tremor and LFEs occur, the presence of high pore fluidpressure as revealed by teleseismic receiver function and regional tomographic studies[Audet et al., 2009, Hansen et al., 2012, Shelly et al., 2006], suggests that fluids may playan important role in allowing slow slip to occur by reducing the effective normal stressto near zero. The peculiar physical conditions represented within the transition fromplate locking to continuous creep have prompted new investigations into the sensitivityof ETS to the tides. It has become increasingly evident in recent years that slow slip,tectonic tremor and LFEs are sensitive to stresses from tidal forces during periods ofETS [Hawthorne and Rubin, 2010, 2013, Klaus, 2012, Lambert et al., 2009, Nakataet al., 2008, Rubinstein et al., 2008, Shelly et al., 2007, Thomas et al., 2012, 2013].Tremors analysis in northern Cascadia shows a clear pulsing of activity with a periodof 12.4h, corresponding to the M2 tide [Rubinstein et al., 2008], with a peak at timesof maximum tidal shear stress [Lambert et al., 2009]. Simulations from laboratoryexperiments [Hawthorne and Rubin, 2010, 2013] are able to reproduce a quasi-sinusoidaltidal modulation of the slip rate, with a maximum moment rate close to the maximumapplied tidal stress. The phase between peak LFE excitation and peak tidal stress is anChapter 1. Introduction 10important factor as it may reflect how the background loading rate influences the timingof LFEs in response to the tidal load. Beeler et al. [2013] proposed a model for whichLFE sources are small seismic patches that fail at a threshold stress on a otherwisecreeping fault. From this model, LFE rate can be explained by adding stressing ratecontributions from the plate motion, the background creep and the tides directly actingon the LFE patches. In chapter 4 of this thesis I analyze the influence of tidal stresseson LFE catalogues in southern Vancouver Island and Washington state. The 1000’s ofdetections that constitute an LFE template and their locations allow a high-resolutionanalysis of the spatial and temporal sensitivity of LFEs to tidally induced up-dip shearstress (UDSS), fault normal stress and their derivatives. We use our observations andthe Beeler et al. [2013] model to constrain fault rheology, such as the effective normalstress on the LFE patches and their dimension, and the characteristic slip for frictionto evolves between two steady state.1.3.4 Rapid tremor reversalsObserved during ETS, Rapid tremor reversals (RTRs) are linear streaks of tremormigration (or LFE detections) moving at high apparent velocities in the “opposite”direction to the main front [Houston et al., 2011]. Thomas et al. [2013] note that RTRsin northern Washington occur almost exclusively during times of positive UDSS. Con-sequently, RTRs may play an important role in slip processes and their characterizationmay lead to improvements in our overall understanding. In chapter 4 of this thesis, Iidentify RTRs in southern Vancouver Island and examine their sensitivity to UDSS.Chapter 2. Blind deconvolution of seismograms regularized via minimum support 11Chapter 2Blind deconvolution ofseismograms regularized viaminimum support2.1 IntroductionThe blind deconvolution problem has a long history in seismology owing to theimportance of separating effects of the source (or “wavelet”) from details of wave prop-agation (the Earth’s “Green’s function”, “impulse response” or “reflectivity”) in studiesof earthquake rupture, earth structure and hydrocarbon exploration. A range of ap-proaches has been considered that includes i) homomorphic deconvolution, ii) higher-order statistics (cumulants), iii) representation theorems and iv) multichannel blinddeconvolution.Homomorphic deconvolution is a non-linear deconvolution technique that satisfiesa generalized principle of superposition.Ulrych [1971] and Tribolet [1978] applied thistechnique to the analysis of seismic wavelets by considering the complex cepstrum ofa seismogram, defined as the complex logarithm of its z -transform. The structure ofthe complex cepstrum is such that it contains the additive contribution of the seismicwavelet and of the impulse response. These contributions may in some instances beseparated by windowing the cepstrum. A narrow window around the quefrency ori-gin, for example, tends to retain the contribution from the source pulse. Using thisChapter 2. Blind deconvolution of seismograms regularized via minimum support 12approach, these authors were able to produce reasonable source estimates but otherssuch as Buttkus (Buttkus [1975]) have noted that the success of wavelet estimation bylow-quefrency windowing is largely dependent on the signal-to-noise ratio of the inputdata, the degree of overlapping of the wavelet cepstrum and impulse response cepstrum,and the choice of the window function.In cumulant methods, as proposed by Sacchi and Ulrych [2000], higher-order cu-mulants of seismic signals are employed in conjunction with spectral factorization torecover seismic wavelets. Higher-order cumulants are obtained by multiplying the nthorder moment of the seismic wavelet by the central lag of the nth order cumulant ofthe reflectivity. Spectral factorization of higher-order cumulants enables estimation ofthe wavelet through the recovery of its minimum-phase and non-minimum-phase com-ponents. This procedure requires the assumption of a non-Gaussian, white reflectivity.Moreover, the authors note that the estimation of higher-order cumulants from data isnot simple. The procedure is computationally intensive and requires long records thatare not always available when working with real data.Weglein and Secrest [1990] proposed a wave-theoretical source estimation methodbased on the Lippmann-Schwinger equation and Green’s theorem for seismic reflectiongeometries. These two equations lead to a relationship between the Fourier transformof the source and the Green’s function that characterizes the propagation within themedium. Results show that this multi-dimensional, deterministic technique allows oneto estimate either an acoustic or an elastic earth wavelet, without any assumptionsconcerning Earth structure. However, it does require that both the field and its normalderivative be measured at the surface.In multichannel blind deconvolution, one exploits the redundancy of multiple seis-mograms sharing a common convolutional element in either sources or Green’s func-tions, and the problem can be reduced to finding the greatest common divisor of twopolynomials (e.g., Gogonenkov [1990]). Rietsch [1997a,b] exploited this relation inimplementing an algorithm that allows one to estimate the source wavelet shared byseveral seismic traces. A matrix is constructed from the autocorrelations and crosscorrelations of the traces. The number of zero eigenvalues of this matrix is equal to thenumber of samples of the wavelet and the eigenvectors associated with these eigenvaluesChapter 2. Blind deconvolution of seismograms regularized via minimum support 13are related to the reflection coefficients. Numerical experiments show that the methodworks well for data that have reasonably met all assumptions (the data share the samewavelet, incorporate different reflectivities, and are not contaminated by noise). In thepresence of noise, however, the performance of this algorithm deteriorates and is diffi-cult to assess.Another method proposed by Mazzucchelli and Spagnolini [2001] extracts the Gr-een’s function by solving an overdetermined linear system in the least squares sensewithin the frequency domain. The results demonstrate that multichannel deconvolutioncan correctly recover the Green’s function when starting from noisy traces convolvedwith a mixed-phase wavelet, without any assumptions on the source’s signature and onthe whiteness properties of the Green’s function. However, the authors remain unclearregarding the implementation of a scaling constraint and indicate that the estimationof the source signature remains difficult when processing large datasets.Here we reexamine the multichannel approach in exploiting the commonality of aconvolutional element within a set of seismograms. Our examples identify the Green’sfunction as this element, but our approach might equally apply to a common source.The extraction of the complementary element (the source in our case) will be accom-plished by solving a system of equations comprising differences of convolutions betweendata and unknowns sources with regularization through minimum support.2.2 Problem formulationWe will cast the blind deconvolution problem in the context of teleseismic body-wave, receiver-side scattering as depicted in Figure 2.1, commonly referred to as receiverfunction analysis [Bostock , 2007, Langston, 1977, Vinnik , 1977]. We consider separaterecordings (or “seismograms”) ui(t) of multiple earthquake sources in close proximity,made at a single, distant receiver. Each earthquake is characterized by a distinct sourcesignature si(t) but the propagation path is the same, so that the seismograms share acommon convolutional element in the Green’s function g(t), that is,ui(t) = si(t) ⋆ g(t), (2.1)Chapter 2. Blind deconvolution of seismograms regularized via minimum support 14where ⋆ denotes convolution, and index i enumerates the source. In actual fact, earth-quake recordings are typically made in 3 orthogonal components of particle velocity.However, for simplicity we shall, for the time beginning, consider ui(t) to represent onlyone component.This component could be one of the vectorial components, e.g. the horizontal compo-nent in the vertical plane that includes both source and receiver, also known as the“radial” component, or it might represent a linear combination of vectorial componentsdesigned to isolate energy within a single scattering mode (i.e. P or S, [Bostock , 2007]).By considering seismograms from any pair of earthquakes si(t), sj(t), we may writethe following expression:ui(t) ⋆ sj(t)− uj(t) ⋆ si(t) = 0. (2.2)We further assume that si(t) should be independent of sj(t); that is, their z-tranformsshould not share zeros, or equivalently, the convolutional elements common to ui(t) anduj(t) should reside only within g(t). This expression can be recast as a homogeneousmatrix system:A x = 0. (2.3)In the case of 2 sources, the system can be written in block form as:[U2 −U1]s1s2= 0, (2.4)where U1 and U2 are Toeplitz matrices containing the seismograms u1(t) and u2(t).If more than 2 sources are included within the system, that matrix A will comprise acombination of blocks Ui, Uj and zero blocks. By imposing a scaling constraint, forexample, by setting the first sample of the first source equal to 1, as s1(1) = 1, thesystem Equation 2.3 becomes:Ax = 0eTx = 1(2.5)Chapter 2. Blind deconvolution of seismograms regularized via minimum support 15PSPpM p Pp M s PsM s P PM sMantleCrust✸✸MantleOuter CoreInner Core✸ +PsM pFigure 2.1: Problem configuration in the context of global seismology. Top panelshows three earthquakes (stars) in close proximity (location separation exaggerated)characterized by source functions s1, s2, s3. P-waves generated by these sources traveleffectively the same propagation path through the Earth’s mantle to a receiver on theEarth’s surface. At the receiver side (inset, bottom panel), the incident P-wave (quasiplanar wavefronts within the mantle indicated in light gray) interact with materialproperty contrasts, e.g. the boundary between crust and mantle or “Moho”, throughthe generation of a suite of scattered phases some of which are initially reflected atthe Earth’s free surface and some of which involve mode conversion from P (solid raypaths) to S (dashed ray paths).Chapter 2. Blind deconvolution of seismograms regularized via minimum support 16with eT = [1 0 · · · 0 0]. The linear system will, in general, be underdetermined andthere are therefore infinitely many solutions. To solve this system, we need to imposemeaningful information about the sources and, in particular, their durations (i.e numberof points in the time series) which are not known a priori.2.3 Methodology2.3.1 Source durationIf the earthquake source durations sj(t) are specified as too short, i.e too few points,we have an overdetermined, inconsistent system. For example, at one extreme if we setthe source durations equal to 1, i.e s1(t) = δ(t) and s2(t) = αδ(t), then u2(t) = αu1(t)which cannot generally be true since u1(t) and u2(t) are required to be independentand so not equal to scaled versions of one another. If the source durations are specifiedas too long, we have an underdetermined, inconsistent system. For example, if we setthe durations equal to the durations of the corresponding seismograms, one solutionis s1(t) = u1(t) and s2(t) = u2(t), since u2(t) ⋆ u1(t) − u1(t) ⋆ u2(t) = 0 . Moreoverany function f(t) convolved with both u1(t) and u2(t) will produce time series that alsosatisfy the system. Consequently we will select as our target the shortest solution whichsatisfies the system Equation 2.5. The solution is thereby regularized as “minimumsupport” and we assume that all the sources are “minimum support”.2.3.2 Introduction of source durations as unknownsWe choose to introduce the source durations as unknowns to be solved for as partof our solution since we do not know the length of the sources a priori. Consequently,it will prove useful to construct a diagonal weighting matrix W(τ) which has for eachsource i the windowing function wi(t, τi) defined as follows:wi(t, τi) =1 for t ≤ τi1− t−τia − 12pia cos2pi(t−τi)a for τi < t ≤ τi + a0 for t > τi + a(2.6)Chapter 2. Blind deconvolution of seismograms regularized via minimum support 17The integer a controls the steepness of the cutoff such that wi(t, τi) tends toward theHeaviside function H(t−τ) as a goes to zero. Note also that this function is continuouswith continuous 1st and 2nd derivatives. Figure 2.2 shows an example of the discretizedweighting function w(t, τ) for different values of a and τ = 4.0 5 10 15 20 25 30 35 40−0.200.20.40.60.811.2time (s)window function w(t,τ) a=1a=5a=10a=15a=20a=25Figure 2.2: Weighting function w(t, τ) for different values of a and τ = 4.The matrix W(τ ), for e.g. 2 sources, can thus be written as:W(τ ) = diag [w1(t, τ1), w2(t, τ2)] . (2.7)The weighting matrix allows us to define a new vector y where x = W(τ )y, such thatthe matrix system in Equation 2.5 is recast as:AW(τ )y ≈ 0eTW(τ )y = 1. (2.8)Here we purposely choose x and y to have dimensions that are considerably longer thanthe anticipated source durations (for example, we may set them to have the length ofu1(t) and u2(t)).Chapter 2. Blind deconvolution of seismograms regularized via minimum support 182.3.3 Combined systemWe now have a combined system where our unknowns are the source durations τand what we will refer to as the “pre-weighted” source functions y(= W−1x). Our aimis to obtain y such that AWy = 0, subject to the constraint eTWy− 1 = 0. We writethe Lagrange function asL(y, λ, τ ) = 12‖ AW(τ )y ‖2 +β2yTΣ−1y − λ(eTW(τ )y − 1) (2.9)where λ is a Lagrange multiplier. Since the matrixWT (τ )ATAW(τ ) is singular or closeto singular, we add a regularization term with Tikhonov parameter β. We introduce themodel prior covariance matrix Σ = E[yyT ] into the regularization term in order to takeinto account the band limited nature of the seismograms. The term 12 ‖ AW(τ )y ‖2refers to the error misfit, β2 yTΣ−1y to the model variance and λ(eTW(τ )y− 1) to thescaling constraint. We shall minimize the Lagrange function L with respect to y, λ andτ .2.3.4 Separation of variablesIt is difficult to minimize the objective function in Equation 2.9 since there is a non-linear dependence on τ . We will use the separation of variables approach suggested byGolub and Pereyra [1973] (see also review in Golub and Pereyra [2003]) to solve thenon-linear problem. Our approach will be to alternately minimize a modified objectivefunction with respect to the linear parameters y and λ, and then with respect to thenon-linear parameters τ . The approach is particularly attractive since the problem isquadratic with respect to y. Thus, given τ we solve a quadratic programing problemfor y and λ,W(τ )TATAW(τ ) +Σ−1 W(τ )eeTW(τ ) 0y(τ )λ(τ )= −∇L (2.10)Chapter 2. Blind deconvolution of seismograms regularized via minimum support 19and rewrite the Lagrange function in terms of τ aloneL(y(τ ), λ(τ ), τ ) = 12 ‖ AW(τ )y(τ ) ‖2 +β2y(τ )TΣ−1y(τ ) (2.11)−λ(τ )(eTW(τ )y(τ )− 1). (2.12)The solution is then obtained by finding a stationary point of the Lagrange functionwith respect to τ .To be more specific, to solve the linear system Equation 2.10 we note that we cansolve for y as a function of λy(τ ) = λ(WT (τ )ATAW(τ ) + βΣ−1)−1WTe. (2.13)Setting yˆ = yλ yieldsyˆ(τ ) = (WT (τ )ATAW(τ ) + βΣ−1)−1WT e, (2.14)and, by multiplying the both sides of Equation 2.13 by eTW from the left, we obtainusing the second equation in Equation 2.10λ = 1eTWyˆ. (2.15)The non-linear problem in τ is addressed by differentiating L with respect to theelements of τ , i.e.dLdτ =∂y∂τ∂L∂y +∂λ∂τ∂L∂λ +∂L∂τ , (2.16)exploiting ∂L∂y = 0 and∂L∂λ = 0 from the linear problem and setting Equation 2.16 tozero, to obtainBT (τ )YT (ATAW(τ )y − λe) = 0, (2.17)Chapter 2. Blind deconvolution of seismograms regularized via minimum support 20where Y = diag(y) and B is a matrix containing the derivatives ∂wi(t,τ )∂τi along itscolumns. For example, for 2 sources, B(τ ) is written as follows:B(τ ) =∂w1(t,τ1)∂τ100 ∂w2(t,τ2)∂τ2. (2.18)However, since the unknowns τi enter this system through W(τ ) in a non-linear fashionwe must solve it using non-linear optimization. In particular, we use Newton’s methodto solve the linear system:H(τ )∆τ = −g(τ ), (2.19)withg(τ ) = (AYB)T (AW(τ )y) − λBTYe). (2.20)The Hessian H(τ ) equalsH(τ ) = C(τ ) +D(τ ), (2.21)whereC(τ ) = (AYB)TAYB+ diag(∑ ∂2∂τ 2 (w(τ , t))YTATAWy), (2.22)andD(τ ) = −λdiag(∑ ∂2∂τ 2 (w(τ , t))Ye). (2.23)Using separation of variables, we proceed iteratively to converge to a solution. Ateach iteration we first recover y and λ by solving the linear system Equation 2.10 andthen solve the linearized system Equation 2.19 to update τ . This procedure is repeateduntil a stopping criterion is satisfied.The separation of variables approach is useful and efficient in our case since theChapter 2. Blind deconvolution of seismograms regularized via minimum support 21number of unknowns in the non-linear problem Equation 2.17 (which is difficult tosolve) is small relative to the linear problem Equation 2.13 (which is relatively easyto solve). Indeed, the number of unknowns in the non-linear problem corresponds tothe number of sources we are trying to obtain whereas the number of unknowns in thelinear problem corresponds to the number of points of all the sources. The alternativeof solving one very large non-linear optimization problem would be extremely difficult.2.3.5 Stopping criterionFrom one perspective, the elements of τ can be regarded as regularization parame-ters that are to be chosen such that the source time functions embodied in x are as shortas possible while still solving the system in Equation 2.8. While the τi are too short,the data misfit (or, more accurately the first term in the LHS of Equation 2.9) is large.When the τi are larger than some set of “optimum” values, the data misfit remainseffectively unchanged. In this sense our solution is one of regularization by minimumsupport. Adopting this rationale, we consider a stopping criterion that invokes Gener-alized Cross Validation (GCV) with τ as the regularization parameter. Note that β inEquation 2.13 is held fixed in our work at a small value that is just sufficient to renderthe linear system non-singular and does not enter the primary regularization enforcedby τ . To construct the GCV parameter GCVmod(τ ) we first recast the constrainedminimization problem as an equivalent least squares problem. Consider, for example,the case where our scaling in Equation 2.5 is applied to the first point of the first source,that is s1(1) = 1. We rewrite the system in Equation 2.5 as:A1s1(1) +A′s′ = 0, (2.24)in whichA =[A1 A′]and x =s1(1)s′. (2.25)Chapter 2. Blind deconvolution of seismograms regularized via minimum support 22Setting s1(1) = 1 and introducing s′ = W′y′, the system Equation 2.24 becomesA′W′y′ = −A1, (2.26)where W′ is equivalent to W with the first column and first row removed. It can beeasily shown, using the leave-one-out lemma [Aster et al., 2005] that the modified GCVparameter employing W′(τ ) as the regularization agent, can be written as:GCVmod(τ ) =m||A′W′y′ +A1||2(trace(I −A′W′A#))2(2.27)withA# = (W′TA′TA′W′ + βΣ′−1)−1W′TA′T , (2.28)where Σ′ is equivalent to Σ with the first column and first row removed, and m is thecolumn dimension of A′. The numerator in Equation 2.27 represents the misfit in theregularized system whereas the denominator represents the effective degrees of freedomsquared [Hansen, 1998].As τ increases, GCVmod(τ ) decreases, though sometimes irregularly, but eventuallybottoms out as opposed to reaching a minimum and subsequently increasing, as inclassical problems. From visual inspection, acceptable solutions are recovered for τcorresponding to the corner of this “L”-shaped GCVmod(τ ) function.Chapter 2. Blind deconvolution of seismograms regularized via minimum support 232.4 Algorithm1. Select a set of initial (short) lengths τ 0 for the sources.2. Solve linear problem:y(τ 0) = λ(WT (τ 0)ATAW(τ 0) + βΣ−1)−1WT (τ 0)e3. Compute threshold GCVmod(τ ) = m||A′W′y′+A1||2(trace(I−A′W′A#))24. While stopping criterion not satisfied4.1 Solve non linear part⇒ ∆τ = −H−1(τ k)g(τ k)4.2 Update τ k+1 = τ k +∆τ4.3 Repeat parts 2. and 3.5. Solution x = Wy if the stopping criterion is satisfied2.5 ResultsThe algorithm explained in previous sections has been tested on several syntheticdatasets. These examples are meant to be simplified representations of teleseismic bodywave seismograms [Bostock , 2007]. Primary body wave phases, such as teleseismic P,are recorded from distant (> 30◦) earthquakes on three-component seismometers withgeometries like those indicated in Figure 2.1. The dominant scattering contributionsarise from structures in the near-source and/or near-receiver regions (Earth’s hetero-geneity tends to be more pronounced in the crust and upper mantle than at greaterdepths). Receiver-side structure which is of interest here is conventionally accessedthrough so-called receiver function analysis [Langston, 1977, Vinnik , 1977] that reliesChapter 2. Blind deconvolution of seismograms regularized via minimum support 24heavily upon the minimum phase nature of the seismogram component in the incidentmode [Bostock , 2007]. In P-receiver function analysis, the S-components of motion aredeconvolved by the P-components to produce a source free impulse response (i.e. thereceiver function) that approximates the S-component of the Green’s function and canbe interpreted in terms of discontinuities in various shear related properties (e.g. S-velocity, S-impedance). However, this approach sacrifices information residing withinthe P-component that pertains to different (compressional) modulli and so a more gen-eral approach that retrieves the full three-component Green’s function is desired. Weexamine multichannel deconvolution to this end. We may consider three-componentseismograms as representing three different components of the Green’s function eachconvolved with the same source, thereby creating a redundancy of up to a factor of 3over the single component system represented in Equation 2.3.The P-component and S-component Green’s functions used in the numerical exper-iments are shown in Figure 2.3 and represent the leading order response of a simplifiedcrustal model comprising a 36 km thick layer (elastic wave velocities and density equalto VP = 6000 m/s, VS = 3464 m/s, ρ = 2700 kg/m3) bounded by free surface aboveand half-space mantle (elastic wave velocities and density equal to VP = 8000 m/s,VS = 4619 m/s, ρ = 3300 kg/m3) below, to an impulsive plane P-wave with slownessp = 0.06 s/km. The large arrival on the P-component at time t = 0 s represents the inci-dent P-wavefield whereas the remaining arrivals on both P and S components comprisefirst-order scattering interactions formed through either direct P-to-S mode conversion(the PMs arrival at time t = 4 s on the S-component) or through reverberation at thefree surface (PpMp at t = 12 s, PsMp at t = 16 s on the P-component and PpMs att = 16 s, PsMs at t = 19 s on the S-component, see Figure 2.1). Because the model is1-D and isotropic we consider only the P-component and the sagittal S-component thatcouples to the incident P-wavefield, since the transverse S-component is zero. Figure 2.4and Figure 2.5 display the convolution of these Green’s functions with 2 sets of 4 sources.The first set of sources comprises random Gaussian time series of length 35.5 seconds(Figure 2.4), whereas the second set is slightly more realistic in that the sources areband limited and extracted from real (decimated) P-wave seismograms (Figure 2.5).Both sets of seismograms are used as input to the blind deconvolution algorithm uponChapter 2. Blind deconvolution of seismograms regularized via minimum support 250 2 4 6 8 10 12 14 16 18 20−0.500.51time (s)amplitudeP−Green’s function0 2 4 6 8 10 12 14 16 18 20−0.100.10.2time (s)amplitudeS−Green’s functionPpMp PsMpPMsPpMsPsMsPFigure 2.3: P-component of Green’s function (upper panel) and S-component (lowerpanel) of Green’s function employed in numerical experiments. Phase labels refer toscattered wave geometries illustrated in lower panel of Figure 2.1.addition of random white noise at 1% of the peak amplitude of P-component seismo-grams (or, equivalently, approximately 10% of the scattered mode amplitudes).The parameters used to perform the inversion are summarized in Table 2.1. The0 5 10 15 20 25 30 35−505P−seismogram from random Gaussian sourceamplitudetime (s)0 5 10 15 20 25 30 35−0.500.51S−seismogram from random Gaussian sourceamplitudetime (s)Figure 2.4: Example P- and S-seismograms generated from the convolution of P-component Green’s function (upper panel) and S-component Green’s function (lowerpanel) with random Gaussian source number 1.Chapter 2. Blind deconvolution of seismograms regularized via minimum support 260 5 10 15 20 25 30 35−1012P−seismogram from band−limited sourceamplitudetime (s)0 5 10 15 20 25 30 35−0.4−0.200.2S−seismogram from band−limited sourceamplitudetime (s)Figure 2.5: Example P- and S-seismograms generated from convolution of P-component Green’s function (upper panel) and S-component Green’s function (lowerpanel) with band-limited source number 1.fig. initial system noise (% ofpeak P)β constraint τ 0(1 : 4) a iter.8 P Gaussian 1 1e-3 e(15)=1 15 20 9109 S Gaussian 1 1e-3 e(15)=1 15 16 206010 P+S Gaussian 1 1e-3 e(15)=1 15 16 120011 P band-limited 1 5e-5 e(23)=1 23 20 123012 S band-limited 1 1e-5 e(23)=1 23 16 90013 P+S band-limited 1 1e-4 e(23)=1 23 20 355Table 2.1: Inversion parameters employed in numerical experiments.β parameter is held fixed in our work at a small value that is just sufficient to renderthe linear system non-singular and does not enter the primary regularization enforcedby τ . Because the P-component in the Green’s function is minimum phase and, inparticular, dominated by the direct arrival at t = 0 s, we may choose our scaling con-straint si(t) = 1 to correspond to a point t on the observed P-component seismogramthat exhibits large relative amplitude. In addition, this point must be chosen so as tofall within a time interval that does not exceed the length of the corresponding sourcefunction. Accordingly, the parameters τ 0 should be chosen to be consistent with thechoice of scaling point, that is they should correspond to times equal to or greater thanthe timing of the scaling constraint. Needless to say, for computational efficiency sake,the τ 0 should be chosen to be as long as possible while remaining shorter than theChapter 2. Blind deconvolution of seismograms regularized via minimum support 27true source durations. The parameter a controls the steepness of the window functionwi(t, τi), (see Equation 2.6). Choice of too small a value for a forces ∆τ to be smallat each iteration and substantially increases computation time. If the a is chosen toolarge, the effective support of x will exceed the true length of the sources.Figure 2.6 to Figure 2.11 display results from 3 classes of experiments summarizedschematically in Table 2.2 using the 2 sets of sources described above. These experi-ments include inversion using P-components alone, S-components alone, and P- and S-components simultaneously. Figure 2.6(a), Figure 2.7(a) and Figure 2.8(a) show resultsexperiment 1 experiment 2 experiment 3A =P-seismogramsA =S-seismogramsA =P-seismogramsS-seismogramsTable 2.2: Schematic representation of the 3 classes of experiment.for the 4 random Gaussian sources. The solutions (dashed grey curves) are retrievedat the corners of the GCVmod parameter versus iteration number curves (black trian-gles in Figure 2.6(b), Figure 2.7(b) and Figure 2.8(b)) corresponding to 910, 2060 and1200 iterations for the 3 experiments, and are compared to the true solution (solid blackcurves). Here we note that the GCVmod(τ ) parameter defined in Equation 2.27 behavesdifferently from that used in more standard (i.e. regularized least-squares) problems asis evident in Figure 2.6(b), Figure 2.7(b) and Figure 2.8(b). As we proceed from lowlevels of regularization (small τ and short source time functions) to high levels (large τand long source time functions), GCVmod(τ ) decreases, though sometimes irregularly,but eventually bottoms out as opposed to reaching a minimum and subsequently in-creasing, as in classical problems. This behavior is presumably due to the fact thatthe misfit (m||A′W′y′ + A1||2) in Equation 2.27 remains constant after the optimum(shortest) solution is reached regardless of the values of τ . We have found from visualinspection that acceptable solutions are recovered for τ corresponding to the corner ofChapter 2. Blind deconvolution of seismograms regularized via minimum support 28this “L”-shaped GCVmod(τ ) function. The GCVmod curve in our experiments did notincrease as τ increased. We believe that this behaviour is due to the special structureof the noise associated with our problem, notably that it enters the RHS of the systemin Equation 2.3. Indeed, in experiments where random noise was added to the righthand side of Equation 2.26, the GCV behaved “normally” decreasing at first and thenincreasing as β increased. Stopping at the corner of the GCVmod curve implies that weselect the solution with the least degrees of freedom that generates a small generalizedcross validation measure. In all experiments source waveforms, durations and relativeamplitudes are well recovered. Figure 2.9(a), Figure 2.10(a) and Figure 2.11(a) showresults of 3 experiments using the set of 4 band-limited sources presented earlier in thesection. The solutions (dashed-grey curve) are retrieved using GCVmod at 1230, 900and 355 iterations (black triangles in Figure 2.9(b), Figure 2.10(b) and Figure 2.11(b)),and again are compared to the true solutions (solid black curves). In this case wehave employed a model prior covariance matrix Σ computed from a large number ofdecimated real P-wave seismograms. Again, source waveforms, durations and relativeamplitudes are well recovered. The 3 classes of experiments yield good results, with asmaller number of iterations required to arrive to a solution for the P+S combination.However, the latter case involves an increase in the number of rows of the matrix A bya factor of 2 and thus increasing the computational burden.2.6 ConclusionWe have described a method that is capable of solving the blind deconvolutionproblem when 2 or more seismograms are available that share a common convolutionalelement (source or Green’s function). For the time series lengths considered here weare able to recover source time functions for noise levels at 1% of the direct P waveamplitude (time 0) (or approximately 10% of scattered wave amplitudes). One mainchallenge in future work will be to improve computational efficiency (execution’ speed)so that realistic seismograms comprising time series with 103−104 points can be accom-modated. At present updating of ∆τ in the non-linear problem proceeds slowly (fewChapter 2. Blind deconvolution of seismograms regularized via minimum support 290 5 10 15 20 25−1−0.500.51Source #1time (s)amplitude 0 5 10 15 20 25−1−0.500.511.5Source #2time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #3time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #4time (s)amplitude true solutionmodeltrue solutionmodeltrue solutionmodeltrue solutionmodel(a)0 500 1000 1500 2000 250010−610−510−4GCVmod parameter# iterationGCVmod (iteration)(b)Figure 2.6: (a)Experiment 1: Inversion using random Gaussian sources and P-components alone. The true solution is shown in solid black and the recovered solutionin dashed grey. (b) GCVmod(τ ) on random Gaussian sources. This example is set for4 random Gaussian source-time functions, a noise-level equal to 1% of the peak P am-plitude, β = 1e − 3, τ 0 = 15, the scaling constraint e(15) = 1, a = 20. At the cornerof the GCVmod(τ ) function (black triangle), a good solution is provided (panel (a)).hours) through 1000’s of iterations to achieve a short duration source solution. Multi-level methods that work first at low frequencies through to higher frequencies may offersolutions in this regard.Chapter 2. Blind deconvolution of seismograms regularized via minimum support 300 5 10 15 20 25−1−0.500.51Source #1time (s)amplitude 0 5 10 15 20 25−1−0.500.511.5Source #2time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #3time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #4time (s)amplitude true solutionmodeltrue solutionmodeltrue solutionmodeltrue solutionmodel(a)0 500 1000 1500 2000 2500 3000 350010−810−710−610−5# iterationGCVmod(iteration)GCVmod parameter(b)Figure 2.7: (a)Experiment 2: Inversion using random Gaussian sources and S-components alone. The true solution is shown in solid black and the recovered solutionin dashed grey. (b) GCVmod(τ ) on random Gaussian sources. This example is set for4 random Gaussian source-time functions, a noise-level equal to 1% of the peak P am-plitude, β = 1e − 3, τ 0 = 15, the scaling constraint e(15) = 1, a = 16. At the cornerof the GCVmod(τ ) function (black triangle), a good solution is provided (panel (a)).Chapter 2. Blind deconvolution of seismograms regularized via minimum support 310 5 10 15 20 25−1−0.500.51Source #1time (s)amplitude 0 5 10 15 20 25−1−0.500.511.5Source #2time (s)amplitude 0 5 10 15 20 25−1−0.500.51Source #3time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #4time (s)amplitude true solutionmodeltrue solutionmodeltrue solutionmodeltrue solutionmodel(a)0 500 1000 1500 200010−710−610−510−4# iterationGCVmod(iteration)GCVmod parameter(b)Figure 2.8: (a)Experiment 3: Inversion using random Gaussian sources and P- and S-components simultaneously. The true solution is shown in solid black and the recoveredsolution in dashed grey. (b) GCVmod(τ ) on random Gaussian sources. This exampleis set for 4 random Gaussian source-time functions, a noise-level equal to 1% of thepeak P amplitude, β = 1e − 3, τ 0 = 15, the scaling constraint e(15) = 1, a = 16.At the corner of the GCVmod(τ ) function (black triangle), a good solution is provided(panel (a)).Chapter 2. Blind deconvolution of seismograms regularized via minimum support 320 5 10 15 20 25−1−0.500.51Source #1time (s)amplitude 0 5 10 15 20 25−1−0.500.51Source #2time (s)amplitude 0 5 10 15 20 25−0.500.51 Source #3time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #4time (s)amplitude true solutionmodeltrue solutionmodeltrue solutionmodeltrue solutionmodel(a)0 500 1000 1500 2000 2500 3000 3500 4000 4500 500010−910−810−710−6# iterationGCVmod(iteration)GCVmod parameter(b)Figure 2.9: (a)Experiment 1: Inversion using band-limited sources and P-components alone. The true solution is shown in solid black and the recovered solutionin dashed grey. (b) GCVmod(τ ) on band-limited sources. This example is set for 4band-limited sources, a noise-level equal to 1% of the peak P amplitude, β = 1e − 5,τ 0 = 23, the scaling constraint e(23) = 1, a = 16. At the corner of the GCVmod(τ )function (black triangle), a good solution is provided (panel (a)).Chapter 2. Blind deconvolution of seismograms regularized via minimum support 330 5 10 15 20 25−1−0.500.51Source #1time (s)amplitude 0 5 10 15 20 25−1−0.500.51Source #2time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #3time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #4time (s)amplitude true solutionmodeltrue solutionmodeltrue solutionmodeltrue solutionmodel(a)0 500 1000 1500 2000 2500 300010−1010−910−810−7# iterationGCVmod(iteration)GCVmod parameter(b)Figure 2.10: (a)Experiment 2: Inversion using band-limited sources and S-components alone. The true solution is shown in solid black and the recovered solutionin dashed grey. (b) GCVmod(τ ) on band-limited sources. This example is set for 4band-limited sources, a noise-level equal to 1% of the peak P amplitude, β = 5e − 5,τ 0 = 23, the scaling constraint e(23) = 1, a = 16. At the corner of the GCVmod(τ )function (black triangle), a good solution is provided (panel (a)).Chapter 2. Blind deconvolution of seismograms regularized via minimum support 340 5 10 15 20 25−1−0.500.51Source #1time (s)amplitude 0 5 10 15 20 25−1−0.500.51Source #2time (s)amplitude 0 5 10 15 20 25−0.500.51 Source #3time (s)amplitude 0 5 10 15 20 25−1−0.500.51 Source #4time (s)amplitude true solutionmodeltrue solutionmodeltrue solutionmodeltrue solutionmodel(a)0 100 200 300 400 500 600 700 80010−910−810−710−6# iterationGCVmod(iteration)GCVmod parameter(b)Figure 2.11: (a)Experiment 3: Inversion using band-limited sources and P- and S-components simultaneously. The true solution is shown in solid black and the recoveredsolution in dashed grey. (b) GCVmod(τ ) on band-limited sources. This example is setfor 4 band-limited sources, a noise-level equal to 1% of the peak P amplitude, β = 1e−4,τ 0 = 23, the scaling constraint e(23) = 1, a = 20. At the corner of the GCVmod(τ )function (black triangle), a good solution is provided (panel (a)).Chapter 3. A comparative study of LFE templates in northern Cascadia 35Chapter 3A comparative study of lowfrequency earthquake templatesin northern Cascadia3.1 IntroductionThe discoveries of non-volcanic tremor by Obara [2002] just over a decade ago andits associations with slow slip [Rogers and Dragert , 2003] and low-frequency earthquakes[Shelly et al., 2006] have opened many new avenues of study in seismology. Unlike reg-ular earthquakes, the LFEs that constitute tremor may repeat 1000’s of times overperiods of a few years. Despite low magnitudes (M < 3) and limited bandwidth (1-10Hz), these signals can be detected at low signal-to-noise ratio (SNR) thresholds and withhigh temporal precision using powerful network correlation techniques [Brown et al.,2008, Gibbons and Ringdal , 2006, Shelly et al., 2006]. Analysis of LFEs, and tremormore generally, can be facilitated through the generation of high SNR LFE templatesassembled by stacking multitudes of aligned, repeating waveforms [e.g. Bostock et al.,2012, Frank et al., 2013, Nowack and Bostock , 2013, Shelly and Hardebeck , 2010]. Thetemplates accentuate impulsive body wave arrivals allowing traveltimes to be measuredmore accurately, polarities to be determined with greater confidence and waveform dis-tortion and scattered signals related to structure to be identified.Chapter 3. A comparative study of LFE templates in northern Cascadia 36In this paper, we generate LFE templates from tremor recorded in northern Van-couver Island and Washington state, thereby extending a previous LFE catalogue forsouthern Vancouver Island [Bostock et al., 2012] through much of northern Cascadia.We then provide empirical and semi-analytical arguments justifying the identificationof LFE templates with “Green’s function” sections corresponding to moment tensorpoint sources exhibiting a step-function time dependence in displacement. We proceedto compare the distribution and excitation of LFE sources across northern Cascadia asrevealed by LFE templates and discuss implications for plate boundary structure andLFE genesis.3.2 DataData employed in this study were collected from a variety of EarthScope sourcesthat include Plate Boundary Observatory (PBO) short-period borehole installations,the transportable USArray (TA), and the Flexible Array experiments CAFE and FACES.In addition, we also employ data from permanent stations of the Canadian NationalSeismograph Network (CNSN) and the Pacific Northwest Seismic Network (PNSN), aswell as portable POLARIS deployments on Vancouver Island. Data from southern Van-couver Island have been described previously by Bostock et al. [2012] and are includedin this work for comparative purposes. The full suite of stations are shown in figure3.1. Data were divided into 4 subarrays, namely northern Washington (NW), southernWashington (SW), northern Vancouver Island (NVI) and southern Vancouver Island(SVI), that were processed independently.3.2.1 Washington stateThe majority of data employed in the analysis of LFEs in Washington state werecollected as part of the Flexible array CAFE experiment [Abers et al., 2009, Calkinset al., 2011]. We separate the dataset into NW and SW components. The formercomprises 25 stations skirting the eastern flanks of the Olympic Peninsula to the westof Puget Sound, whereas the latter includes 25 stations running approximately east-westto the southwest of Puget Sound. Figure 3.1d displays the distribution of these stations.Chapter 3. A comparative study of LFE templates in northern Cascadia 37For the NW dataset, major tremor episodes in 2007, 2008, 2009, 2010 and 2011 wererecorded, but only data from the 2007 and 2008 tremor episodes were available for theSW subarray. Episodic tremor and slip episodes such as these typically last up to ∼3weeks in Washington and southern Vancouver Island [Rogers and Dragert , 2003].3.2.2 Vancouver IslandTremor data for the NVI dataset were assembled from a portable POLARIS de-ployment on northern Vancouver Island that comprised 27 broadband seismometersdeployed along two mutually perpendicular arms [Audet et al., 2008], as shown in fig-ure 3.1b. Two tremor episodes in 2006 and 2007 were available for analysis, each episodelasting under 1 week in duration. Analysis of LFEs for SVI relied heavily on a line ofsome 10 stations extending from the west coast of Vancouver Island (PFB) to the GulfIslands (SNB) but also incorporated an additional >20 stations from surrounding areas.Tremor episodes in 2003, 2004 and 2005 were recorded at a majority of the stations,but data from subsequent episodes through 2012 have been incorporated for stationsas available.3.3 Data ProcessingAll data were divided into 24 hour-long segments, band-pass filtered between 1 and8 Hz, resampled to 40 sps, and subjected to three distinct data processing steps to obtainLFE templates. These steps comprised i) network autocorrelation to identify pairs ofrepeat LFEs [Brown et al., 2008], ii) waveform-correlation cluster analysis [Rowe et al.,2002] to sort LFE detections into initial templates based on waveform similarity, andiii) network correlation [Gibbons and Ringdal , 2006, Shelly et al., 2006] and stackingto increase detections and improve LFE template SNR. Bostock et al. [2012] providedetails on processing of the SVI dataset, so we consider only NW, SW and NVI datasetsbelow.Chapter3.AcomparativestudyofLFEtemplatesinnorthernCascadia38 125 oW 30’ 124 oW 30’ 123oW 30’ 1 2 2o W 30’ 47 oN 30’ 48 oN STW B007B 04 ALRIVSQM OPC B001BLN B S1 1 FACBW 02 0B943B 01 3W030C04ADOSEHDW W040GNW PL11W060SMW N050N080N060N015N020NLWAN030WISHD 03 AS015N040S 02 0B017B019E03AW090W080S030W070S040B018S050S055N070 20’ 127oW 40’ 20’ 126oW 50’ 50 oN 10’ 20’ 30’ VI10VI11VI12TASBMAYBVI52VI53VI54VI55VI56VI57WOSBVI01VI02VI03VI04VI05VI06VI07VI08VI09NCRBVW01VW02VW03VI30VI50 126 oW 125 oW 124oW 123oW 1 2 2o W 48 oN 20’ 40’ 49 oN 20’ 40’ B012OZB BMSBB 92 8THABB927MGB TFRBYOUBB926SHB BIB GHNBNLLBGOBBPHYBBPCBSHDBSSIBA04AENGBSNB GOWBSILBB011PGC KELBM GC BTWKBLZB TSJBTWBBTWGBPFB S OK B SH VBKHVBVGZ SQM B001OPC LRIVB007B 04 AB003B004CPLBLCBC JRBCGLBC 126oW 124oW 1 2 2o W 1 2 0o W 130 oW 128 oW 46 oN 47 oN 48 oN 49 oN 50 oN 51 oN 52 oN N oo t kaF a ul tJUAN DE FUCAPLATE4 0 mm /y ea rBDCA BC DFigure 3.1: Distribution of stations used in the study of LFE templates in Washington state and Vancouver Island. Top left panel a) showsnorthern Cascadia region including plate boundary and color-coded station locations (NVI-orange, SVI-yellow, NW-red, SW-green, shared NW/SW-purple). Remaining panels b), c), d) show map insets of station locations for the NVI, SVI and collective NW and SW subarrays, respectively.Chapter 3. A comparative study of LFE templates in northern Cascadia 393.3.1 Network autocorrelationFor the NW dataset, we applied network autocorrelation to 4 tremor episodes(January 17-30 2007, May 01-23 2008, August 07-24 2010 and August 04-23 2011) usingthe 7 three-component stations B001, SQM, BS11, W020, W040, GNW and PL11.For the SW dataset, we performed network autocorrelation using 2 tremor episodes(January 14-February 01 2007, May 01-25 2008) and employed the 7 three-componentstations N050, S030, W070, N060, PL11, S040 and N070. Network autocorrelation forthe NVI dataset was applied to two tremor episodes (September 05-11 2006, June 13-182007) for different combinations of 7 of the 10 three-component stations VI10, VW03,VI11, VW02, VI52, VI53, VI05, VI08, VI06, VI04.Each of these data sets was independently analysed on an hour-by-hour basis. Eachhour-long segment was divided into 15 s windows lagged by 0.5 s to produce a total of7170 individual windows per hour. Each window was correlated with all other windowsin the same 1-hour record for an individual (station-component) channel. The resultingtime-series from all stations and all components corresponding to the same hour werestacked to create “network autocorrelation” records. When the network autocorrelationexceeded 8 times the median absolute deviation [Shelly et al., 2006], we registered thecorresponding window pairs as LFE detections. This catalogue was then culled toexclude overlapping detections and retain only those with high SNR. After processingall data, we had selected 4915, 3306 and 4267 initial detections for further analysis fromthe NW, SW and NVI datasets, respectively.3.3.2 Waveform-correlation cluster analysisWe proceeded to combine all channels corresponding to a given detection into a sin-gle “super” trace and cross-correlated all such traces against one another to determinemaximum correlation coefficients and corresponding lags. The correlation coefficientswere used to populate a similarity matrix employed within a hierarchical cluster anal-ysis, allowing the detection waveforms to be grouped into clusters. Waveforms forall channels (stations/components) available for detections within a given cluster wereshifted and stacked to produce an initial LFE template. This procedure resulted in 224,70 and 54 initial templates for the NW, SW and NVI datasets, respectively.Chapter 3. A comparative study of LFE templates in northern Cascadia 403.3.3 Network cross-correlationIterative network cross-correlation and stacking was used to register further detec-tions and improve template SNR. We performed network cross-correlation by choosingsubsets of 7 to 10 high SNR, three-component stations available within a template andscanning through all available tremor dates. As before, when the summed networkcross-correlation coefficient exceeded 8 times the median absolute deviation, now for a24-hour period, we logged detections. New templates were formed by stacking wave-forms (normalized to unit maximum amplitude across three components) for all newlyregistered detections. After several iterations of network cross-correlation and stacking,we obtained final sets of templates suitable for location and waveform analysis. Ourfinal suites of templates number 122, 54 and 47 for the NW, SW and NVI datasets,respectively. Each template comprises 100’s to 1000’s of independent detections. Fig-ure 3.2 shows examples of LFE templates from each of the NW, SW, NVI datasets. Sarrivals can be clearly seen on both horizontal and vertical components and P arrivalscan be clearly seen on the vertical component.3.4 LFE templates as empirical Green’s functionsIn figure 3.3 we plot P and S waveforms at stations located near the center ofNVI (top panels), NW (middle panels) and SW (bottom panel) arrays. In each panel,waveforms are ordered in increasing epicentral distance (ranging from 1 to 71 km, seefollowing section) and aligned with respect to the dominant P or S arrival. Bostocket al. [2012, Figure 6] present a similar figure for the SVI data. In contrast to thisprevious study, we have applied a 90◦ phase shift (Hilbert transform) to the templatewaveforms that accomplishes partial transformation from particle velocity to particledisplacement without altering the amplitude spectrum. The effects of post-critical scat-tering interactions that induce complex waveform distortions are minimized for thesestation selections because the majority of the propagation paths are near-vertical [Boothand Crampin, 1985].As remarked by Bostock et al. [2012], the P and S waveforms for the full selectionChapter 3. A comparative study of LFE templates in northern Cascadia 41Time [s]NVI #040 NORTH0 20 40 60VI57VI56VI55VI30VI50MAYBVI52VI54VI53TASBVI12VI11VI10VW03VW02VW01VI08VI09WOSBVI07VI06VI05VI04VI03VI02NCRBVI01Time [s]NVI #040 EAST0 20 40 60VI57VI56VI55VI30VI50MAYBVI52VI54VI53TASBVI12VI11VI10VW03VW02VW01VI08VI09WOSBVI07VI06VI05VI04VI03VI02NCRBVI01Time [s]NVI #040 VERTICAL0 20 40 60VI57VI56VI55VI30VI50MAYBVI52VI54VI53TASBVI12VI11VI10VW03VW02VW01VI08VI09WOSBVI07VI06VI05VI04VI03VI02NCRBVI01Time [s]SW #053 NORTH0 20 40 60N015N020NLWAN030WISHD03AS015N040S020B017B019E03AW090W080N050SMW W060S030W070N060PL11S040B018S050S055Time [s]SW #053 EAST0 20 40 60N015N020NLWAN030WISHD03AS015N040S020B017B019E03AW090W080N050SMW W060S030W070N060PL11S040B018S050S055Time [s]SW #053 VERTICAL0 20 40 60N015N020NLWAN030WISHD03AS015N040S020B017B019E03AW090W080N050SMW W060S030W070N060PL11S040B018S050S055Time [s]NW #210 NORTH0 20 40 60STW B007B04ALRIVSQM OPC B001BLN BS11FACBW020B943B013W030C04ADOSEHDW W040GNW PL11W060SMW N050N080N060Time [s]NW #210 EAST0 20 40 60STW B007B04ALRIVSQM OPC B001BLN BS11FACBW020B943B013W030C04ADOSEHDW W040GNW PL11W060SMW N050N080N060Time [s]NW #210 VERTICAL0 20 40 60STW B007B04ALRIVSQM OPC B001BLN BS11FACBW020B943B013W030C04ADOSEHDW W040GNW PL11W060SMW N050N080N060ACBFigure 3.2: North (left), east (center) and vertical (right) components of 3 LFEtemplates from the a) NVI (template 040), b) NW (template 210) and c) SW (template053) subarrays.of templates at a given station display remarkable uniformity across the range of epi-central distances. Furthermore, the dipolar pulses observed in the earlier study and thebandlimited zero-phase pulses evident in figure 3.3 imply that the LFE templates canbe considered as empirical Green’s functions originating from a moment tensor pointsource with a step-function time dependence in displacement. Hilbert transformation ofChapter 3. A comparative study of LFE templates in northern Cascadia 42the particle velocity records is practically useful because it aids in the identification of(band-limited) arrivals for location (section 3.5) and for waveform matching as requiredin moment tensor inversion (section 3.6) and recovery of structure from scattered waves[Nowack and Bostock , 2013].The assertion that LFE templates produced by the iterative correlation-detection-stacking procedure can be considered as empirical Green’s functions can be further justi-fied as follows. Consider the single-channel cross correlation c(t) of a (e.g. band-limited)Green’s function estimate gˆ(t) with a seismogram u(t) that is itself the convolution ofa source s(t) with the corresponding true Green’s function g(t):c(t) = gˆ(t)⊗ u(t) = gˆ(t)⊗ g(t) ⋆ s(t) = Φˆg(t) ⋆ s(t),where ⊗ and ⋆ denote correlation and convolution, respectively. Φˆg(t) is an estimateof the true autocorrelation of g(t) and so will peak near lag t = 0. The timing ofthe maximum (or maxima) of c(t) will depend on s(t) but, owing to the timing ofmaximum in Φˆg(t), will occur at or near the time(s) at which s(t) possesses its maximum(maxima), say t = tmax. Note that Φˆg(t) may contain significant subsidiary maximaat, for example, times equal to the S-P time if both P- and S-waves project onto thegiven channel with comparable amplitude and polarity. The contaminating influence ofsecondary peaks, any bias of the autocorrelation estimate maximum away from t = 0and other forms of noise will be mitigated in the selection of tmax when the normalizedcross correlations c(t)/√∫gˆ(τ)2dτ∫u(τ)2dτ for many channels are combined withinnetwork correlation detection [Gibbons and Ringdal , 2006, Shelly et al., 2006].Once the detections are logged as described in section 3.3, the stack of correspondingseismograms ui(t) for a given channel, each shifted by the corresponding detection timetmaxi , can be written as∑iui(t− tmaxi ) = g(t) ⋆∑isi(t− tmaxi ).As the number of detections increases, the sum of band-limited source functions si(t−tmax), shifted such that their maxima align but assumed to be otherwise random, willtend toward a scaled, filtered delta function. The sum of shifted ui(t) thus becomes aChapter 3. A comparative study of LFE templates in northern Cascadia 43Time [seconds]SW−N060 − N: S10 20 30 40 501011121314151617181920SW−N060 − E: S10 20 30 40 501011121314151617181920SW−N060 − Z: S10 20 30 40 501011121314151617181920Increasing Epicentral DistanceTime [seconds]SW−W070 − N: P10 20 30 40 501011121314151617181920SW−W070 − E: P10 20 30 40 501011121314151617181920SW−W070 − Z: P10 20 30 40 501011121314151617181920Time [seconds]NW−W020 − N: S20 60 1001011121314151617181920NW−W020 − E: S20 60 1001011121314151617181920NW−W020 − Z: S20 60 1001011121314151617181920Time [seconds]NW−B013 − N: P20 60 1001011121314151617181920NW−B013 − E: P20 60 1001011121314151617181920NW−B013 − Z: P20 60 1001011121314151617181920Time [seconds]NVI−VI11 − N: P20 401011121314151617181920NVI−VI11 − E: P20 401011121314151617181920NVI−VI11 − Z: P20 401011121314151617181920Time [seconds]NVI−VW01 − N: S20 401011121314151617181920NVI−VW01 − E: S20 401011121314151617181920NVI−VW01 − Z: S20 401011121314151617181920abcFigure 3.3: Three-component LFE template waveforms for individual stationsaligned on direct P/S phases and plotted as functions of epicentral distance. a) NVIwaveforms aligned on P for station VI11 and S for station VW01. b) NW waveformsaligned on P for station B013 and S for station W020. c) SW waveforms aligned onP for station 2070 and S for station N060. Red/blue polarities are positive/negative,respectively. Note simple zero-phase signatures of direct P and S.scaled, band-limited approximation to g(t) that can be used as an improved estimategˆ(t) to log further detections.Chapter 3. A comparative study of LFE templates in northern Cascadia 44Approximate deconvolution through stacking of phase-normalized seismograms hasbeen applied previously to long-period body waves [Shearer , 1991] and broadband tele-seismic P-waves [Kumar et al., 2010] using global earthquakes; however the relativeinfrequency of regular seismicity limits Green’s function retrieval to 1-D estimates. Incontrast, LFEs that repeat 1000’s of times over periods of a few years allow fully 3-Dempirical Green’s functions to be assembled.3.5 LFE template locationsA large proportion of stations represented within LFE templates display unam-biguous zero-phase, impulsive P and S arrivals that can be timed and used to establishrepresentative locations. We perform two locations, one using a standard linearized in-version [Hyp2000, Klein, 2002] and one using the double difference location algorithm[hypoDD, Waldhauser , 2001], both distributed by the United States Geological Survey.Figure 3.4 shows a map of LFE Hyp2000 epicenters (see supplementary figure A1 inAppendix A for hypoDD epicenters). Superimposed on this map are the 20, 30 and 40km depth contours to the top of the subducting Juan de Fuca plate modelled by Audetet al. [2010] and McCrory et al. [2012].Using the initial Hyp2000 locations, we apply the hypoDD algorithm to a combi-nation of ordinary phase picks from our LFE catalog and differential travel times fromphase correlation of P- and S-waves. The parameters used to perform the inversion aresummarized in supplementary Table A1 in Appendix A. These parameters are set toproduce a dynamic weighting scheme to optimize the least-squares solution. Solutionsare found by iteratively adjusting the vector difference between nearby hypocentralpairs, with the locations and partial derivatives updated after each iteration. Eventslacking close neighbours are automatically removed in this procedure such that 91%,92%, 77% and 72% of the original events remain in hypoDD solutions for the NVI,SVI, NW and SW datasets, respectively. Mean, nominal uncertainties in horizontaland vertical locations as supplied by Hyp2000 for the individual subarrays are 1.8 km,2.16 km (NVI), 1.1 km, 1.9 km (SVI), 2.0 km, 1.3 km (NW), and 3.7 km, 3.5 km (SW).Chapter3.AcomparativestudyofLFEtemplatesinnorthernCascadia45 130 oW 128 oW 126oW 124oW 1 2 2o W 1 2 0o W 46 oN 47 oN 48 oN 49 oN 50 oN 51 oN 52 oN C′D′E′AA′B′HH′II′JJ′CD F′GG′BEFJUAN DE FUCAPLATE4 0 mm /y ea rN oo t kaF a ul t 20’ 127oW 40’ 20’ 126oW 50’ 50 oN 10’ 20’ 30’ 126 oW 125 oW 124oW 123oW 1 2 2o W 48 oN 20’ 40’ 49 oN 20’ 40’ 125 oW 30’ 124oW 30’ 123oW 30’ 1 2 2o W 30’ 47 oN 30’ 48 oN A BC DFigure 3.4: Maps of LFE locations computed using the Hyp2000 software organized as in figure 3.1. Orange, red and green diamonds are LFElocations from the NVI, NW and SW arrays, respectively. Yellow diamonds are locations of LFEs in southern Vancouver Island described byBostock et al. [2012]. Cyan and blue lines indicate the 20, 30 and 40 km depth contours to the top of the subducting Juan de Fuca pate modelledby Audet et al. [2010] and McCrory et al. [2012], respectively. Black dots represent regular earthquake locations for the period 1985-2012.Chapter 3. A comparative study of LFE templates in northern Cascadia 46LFE epicenters in northern Cascadia, shown in figures 3.4 and S1 fall within thegeneral tremor epicenter distributions previously mapped by Kao et al. [2009, Figures7,9] and Wech et al. [2009, Figures 1,3]. In the southern Vancouver Island (SVI) andWashington state (NW and SW) regions, LFE epicenters are bounded by the 25 and38 km slab depth contours for the Audet et al. [2010] model and the 32 and 45 kmcontours for the McCrory et al. [2012] model. In northern Vancouver Island the LFEtemplates map close to the 30 km contour for both models. As noted by Bostock et al.[2012] for southern Vancouver Island LFEs and by Kao et al. [2005] for tremor moregenerally, LFEs epicenters tend to avoid regions with higher levels of regular seismicity,shown as black dots in figure 3.4 for the period 1985-2012.Figure 3.5 plots the hypoDD locations of 10 depth profiles defined in figure 3.4(analogous plots for the hyp2000 locations can be found in supplementary figure A2in Appendix A). Note that A-A’ and B-B’ profiles below southern Vancouver Islandwere previously presented by Bostock et al. [2012]. Profiles are constructed using binsthat extend ± 25 km to either side, and include the Audet et al. [2010] slab modelquadratically interpolated through the 20,30 and 40 km contours and the McCrory et al.[2012] model linearly interpolated through 5 km depth intervals between 20 and 80 km.This figure further emphasizes the segregation of LFEs from regular earthquakes. LFEstend to lie several km on average above intraplate earthquakes where their epicentersoverlap. The two slab models bracket the LFE hypocenters below southern VancouverIsland and northern Washington (profiles A,B,C), the Audet et al. [2010] model fromabove and the McCrory et al. [2012] model from below. Hypocenters beneath SVI rangebetween 29-40 km depth whereas those below Washington state fall between 31-46 kmdepth (NW) and 32-49 km depth (SW). Hypocenters along the SW profiles (E,F,G)coincide more closely with the McCrory et al. [2012] model, whereas hypocenters forLFEs for NVI profiles I,J are better aligned with the Audet et al. [2010] model wherehypocenters fall between 33-39 km depth.Chapter 3. A comparative study of LFE templates in northern Cascadia 47−100 −50 0 50 100 150 200020406080Depth [km]H H′−100 −50 0 50 100 150 200020406080Depth [km]I I′−100 −50 0 50 100 150 200020406080Depth [km]J J′−100 −50 0 50 100 150 200020406080Depth [km]A A′−100 −50 0 50 100 150 200020406080Downdip Distance [km]Depth [km]B B′−100 −50 0 50 100 150 200020406080C C′−100 −50 0 50 100 150 200020406080D D′−100 −50 0 50 100 150 200020406080E E′−100 −50 0 50 100 150 200020406080F F′−100 −50 0 50 100 150 200020406080Downdip Distance [km]G G′Figure 3.5: Depth profiles of seismicity in Vancouver Island and Washington state.Profile locations are identified in figure 3.4a. Black dots represent regular earthquakelocations for the period 1985-2012. Orange, yellow, red and green diamonds are LFEtemplate locations determined using hypDD for NVI, SVI, NW and SW arrays, re-spectively. Dashed cyan and blue lines represent depth estimates to the top of thesubducting Juan de Fuca plate from Audet et al. [2010] and McCrory et al. [2012]models, respectively.3.6 LFE template moment tensorsIn our previous effort [Bostock et al., 2012] to determine focal mechanisms of LFEson southern Vancouver Island, we employed P-polarity estimates that were derivedthrough correlation of vertical component seismograms with a reference pulse. Thatwork suggested that LFE mechanisms comprised a mixture of thrust and strike slipfaulting. We improve our analysis in the present study in several ways. First, 90◦ phaserotation of particle velocity seismograms provides a more robust means of producingzero-phase pulses and ascertaining the polarity of P- and S-arrivals at smaller epicen-tral distances than correlation with a reference pulse. Second, we have incorporateddata from several previously unavailable stations along the southern coast of VancouverChapter 3. A comparative study of LFE templates in northern Cascadia 48Island (CPLB,LCBC,JRBC,GLBC) that significantly improve coverage in the south-western quadrant. Finally, we now exploit the identification of the LFE templates withempirical Green’s functions to perform a moment tensor inversion that, in addition toP-wave polarities, leverages constraints from S-wave polarities and relative amplitudesof both P- and S-waves.Our moment tensor inversion procedure incorporates elements from Kikuchi andKanamori [1991] and involves the following steps: i) each 3-component LFE templateseismogram is normalized to unit, maximum absolute amplitude and phase rotated by90◦ to aid in identication and picking of zero-phase, primary arrivals; ii) locations aredetermined from traveltime picks corresponding to those P- and S-arrivals judged todisplay unambiguous, zero-phase waveforms, iii) ray theoretic synthetic seismograms[Cerveny et al., 1987] comprising direct P- and S-arrivals are generated for a basis of 6independent moment tensors and normalized/phase rotated/filtered to match the spec-tral properties of the data; iv) synthetic and observed P and S waveforms are alignedbased on their amplitude extrema to account for unmodelled velocity structure; v) least-squares inversion is performed for the coefficients of the moment tensor basis that bestexplain the relative amplitudes and polarities of the primary LFE arrivals; and vi) thebest double couple mechanism is extracted from the general moment tensor solutionusing Newton optimization.Figure 3.6 shows an example of the waveform matches derived for SVI template005 with a double couple solution. The effect of a steeply inclined nodal plane roughlyparallel to strike is evident in the change in P-wave polarity that occurs near stationLZB on the vertical component and the large, negatively polarized S-arrivals on thenorth and east components at e.g., stations TWBB to PGC that straddle the plane.Best-fit double couple solutions are plotted in figure 3.7 for a geographically representa-tive subset of SVI templates. SVI solutions are generally well modelled as double couplesources as indicated by F-tests that reveal no significant improvement in fit affordedby deviatoric or full moment tensor solutions. The majority of double couple solutionsare characterized by thrust mechanisms oriented in a northeasterly direction, althoughthere are two groupings of 3-4 LFE locations each near 48◦12’ N, 123◦6’ W and 48◦42’N, 124◦12’ W that appear to show somewhat consistent strike-slip components.Chapter 3. A comparative study of LFE templates in northern Cascadia 49In addition to the correct identification of wavelets on LFE templates as direct Pand/or S, our ability to constrain the moment tensor depends on event-station geome-try as quantified by the condition number of the normal equations. The extensive arealaperture of the NVI array and a localized distribution of LFEs near its center lead tocondition numbers (mean=17, median=22 for the deviatoric solution, see supplemen-tary figure A3 in Appendix A) that are less than those for SVI (mean=68, median=33),although this measure of solution quality does not account for the template SNR that issuperior for SVI. Double couple solutions for the NVI templates (figure 3.8) like thosefor SVI tend to be shallow thrust in nature. The same characterization holds for theSW solutions (figure 3.9) with mean and median condition numbers of 17 and 27, re-spectively although there is greater variability in their orientation. In contrast, momenttensor solutions for NW templates (figure 3.10) are marked by large condition numbers(mean=366, median=89) that manifest the quasi-linear distribution of stations, anddisplay still greater variability in orientation. Each of figures 3.7-3.10 contains an insetdisplaying the best double couple solution (i.e. with linear vector dipole componentremoved) determined from the average of individual double couple solutions for alltemplates within the respective subarray. The corresponding averaged plate motionvector determined from data in McCaffrey et al. [2007] is superposed on these plots.Supplementary figure A4 in Appendix A provides an alternative representation of thefocal mechanism distributions as histograms of strike, rake and dip for the individualsubarrays assuming that the nodal plane with shallower dip represents the fault plane.3.7 Discussion and conclusionsThrough the generation of LFE templates using data from EarthScope sources(Transportable array, PBO, FlexArray), POLARIS deployments and permanent net-works (CNSN, PNSN), we have extended documentation of Cascadia LFEs from south-ern Vancouver Island south into Washington state and north to northern VancouverIsland. Not surprisingly, LFE epicenters in northern Cascadia fall within the tremorepicentral distributions previously mapped by, e.g. Kao et al. [2009], Wech et al. [2009],Chapter 3. A comparative study of LFE templates in northern Cascadia 50and tend to avoid regions of denser, regular seismicity as reported previously for tremor[Kao et al., 2005] and for LFEs below southern Vancouver Island [Bostock et al., 2012].In particular, LFEs typically define a surface several km above the upper envelope ofintraplate earthquakes and several km below the lower envelope of overlying crustalseismicity where their epicenters overlap.LFE hypocenters generally parallel but do not coincide precisely with either of tworecent models for the plate boundary in Cascadia. In fact, the hypocenters frequentlylocate between the plate boundary estimates, lying below the Audet et al. [2010] modeland above the McCrory et al. [2012] model. The two models are based on differentassumptions and data sets, and both may be subject to bias. The Audet et al. [2010]model relies on the identification of the plate boundary with the top of a pronouncedlow-velocity zone (LVZ) that occurs throughout Cascadia and which has been mostrecently interpreted to be upper oceanic crust [Bostock , 2013, Hansen et al., 2012].This model was generated for the purposes of cross-Cascadia comparisons between thegeometry of the LVZ and tremor epicentral distributions [e.g. Wech and Creager , 2011,figure 1]. It was derived using the timing of scattered teleseismic P-phases and thesimplifying assumption of a homogeneous overriding plate with fixed P-velocity of 6.5km/s but locally variable VP /VS ratio as opposed to the 1-D, fixed VP /VS model usedfor LFE location. The differences in underlying velocity models may result in depthbiases that are locally significant.The McCrory et al. [2012] model was constructed by synthesizing depth informationfrom intraplate earthquake locations and regional seismic velocity profiles. The authorsinferred the top of the Juan de Fuca slab to lie near the upper surface of intraplateseismicity where present, and weighted these seismicity constraints more highly thanstructural information derived from velocity profiles in areas where both sources of in-formation were available. Plate boundary depth estimates based on this approach maybe biassed deep by ∼ 7 km if intraplate seismicity resides near the base of the subduct-ing crust as has been inferred by Shelly et al. [2006] for southwest Japan. Comparison ofhypocenters from southwest Japan with those from northern Cascadia [Bostock et al.,2012, figure 8; see also figures 3.5,A2 herein] points to similar geometrical relationsbetween LFEs, LVZs and intraplate seismicity, and implies that the structural controlsChapter 3. A comparative study of LFE templates in northern Cascadia 51on seismogenesis in the two regions are the same.The only region where the Audet et al. [2010] model maps the plate boundary todeeper levels than the McCrory et al. [2012] model, is along profiles I,J in northernVancouver Island. In this region, the Explorer microplate is interpreted to be detach-ing from the Juan de Fuca plate along the Nootka fault and its NE landward extension[e.g Audet et al., 2008, Braunmiller and Nabelek , 2002]. The majority of LFEs herelie along a relatively flat trajectory between 35-38 km depth coinciding with the Audetet al. [2010] model along profiles I,J, although they lie significantly deeper than thatmodel on profile H immediately to the north. The LFE hypocenters are also slightlydeeper than the 25-35 km depths quoted for tremor in this region by Kao et al. [2009].We have outlined simple arguments justifying the identification of LFE templateswith empirical Green’s functions, thereby facilitating their treatment in waveform in-versions for, e.g., moment tensor solutions. Note that we do not mean to imply thatall LFEs are characterized by single point-source, step displacements in time. Exami-nation of individual LFEs detected via network correlation often does reveal impulsive,albeit noisy, signals resembling the template waveforms. Just as frequently, however,the signals display more complex temporal dependencies consistent, for example, withrapid tremor streaking that has been documented in northern Washington using beam-forming techniques [Ghosh et al., 2010].Our examination of source mechanisms from the SVI subarray using moment ten-sor inversion improves on our analysis of P-wave polarities in Bostock et al. [2012] byincluding more stations from the previously poorly sampled southwestern quadrant, byemploying S-wave polarities and relative amplitudes across individual station compo-nents, and by the use of a 90◦ phase rotation to facilitate phase identification within thebandlimited signals. SVI moment tensor solutions are generally well constrained andthe large majority (∼ 90%) of mechanisms are consistent with shallow thrust faultingin the direction of relative plate motion. The variability in focal mechanisms increasesprogressively through the NVI, SW and NW subarrays. We suspect that this variabil-ity may be due partly to poorer station coverage (for NW in particular) and partly tolower SNR resulting from smaller numbers of contributing detections. Diminished SNRrenders it difficult to accurately isolate phases, especially S, within bandlimited data.Chapter 3. A comparative study of LFE templates in northern Cascadia 52S is particularly problematic because, at larger epicentral distances, generation of post-critically reflected P at the free surface significantly distorts S-waveforms [Booth andCrampin, 1985] such that only the SH component retains its original source signature.In addition, strong anisotropy is known to occur in some areas [Bostock and Chris-tensen, 2012] and may also contaminate waveforms through splitting. Notwithstandingvariability in focal mechanism consistency across the different subarrays, the averagedouble-couple solution for each (see insets in figures 3.7-3.10) is generally consistentwith shallow thrusting in the direction of plate motion.These subarray-averaged moment tensors and the preponderance of individualthrust mechanisms for the SVI, NVI (and to a lesser extent SW) subarrays approx-imately aligned with the plate motion direction leads us to suspect that shallow thrustfaulting may prevail throughout the northern Cascadia region, as has been argued forSW Japan by Ide et al. [2012] and central Mexico by Frank et al. [2013]. In so doing,we interpret variability in focal mechanisms in figures 3.7-3.10 (for the NW array inparticular) as due to variations in SNR and conditioning of the inverse problem. Weacknowledge, however, the possibility that, locally, LFE mechanisms may depart fromthe shallow thrust orientation dependent, for example, upon structure in the downgoingplate. The likelihood that a majority of, if not all, LFE mechanisms are shallow thrustwould weaken the argument by Bostock et al. [2012] that LFEs are distributed througha plate boundary shear zone coinciding with the (3-4 km thick) LVZ. Although ournominal depth location uncertainties (∼ ±2 km) do not allow us to address this issuedirectly, recent work by Nowack and Bostock [2013] employing scattered waves from aselection of templates requires LFEs to occur < 1 km below the top of the LVZ. Thisconstraint together with the moment tensor results presented here are consistent withan origin for LFEs as shear slip along a relatively sharp plate boundary atop the LVZ.Chapter 3. A comparative study of LFE templates in northern Cascadia 53Time [s] Z−Component5 10 15PFB TWGBTWBBTSJBLZB TWKBMGCBKELBPGC SSIBSILBSNB VGZ KHVBSHVBGLBCSOKBJRBCLCBCYOUBTime [s] N−Component5 10 15PFB TWGBTWBBTSJBLZB TWKBMGCBKELBPGC SSIBSILBSNB VGZ KHVBSHVBGLBCSOKBJRBCLCBCYOUBTime [s] E−Component5 10 15PFB TWGBTWBBTSJBLZB TWKBMGCBKELBPGC SSIBSILBSNB VGZ KHVBSHVBGLBCSOKBJRBCLCBCYOUBTime [s] N−Component5 10 15PFB TWGBTWBBTSJBLZB TWKBMGCBKELBPGC SSIBSILBSNB VGZ KHVBSHVBGLBCSOKBJRBCLCBCYOUBTime [s] E−Component5 10 15PFB TWGBTWBBTSJBLZB TWKBMGCBKELBPGC SSIBSILBSNB VGZ KHVBSHVBGLBCSOKBJRBCLCBCYOUBTime [s] Z−Component5 10 15PFB TWGBTWBBTSJBLZB TWKBMGCBKELBPGC SSIBSILBSNB VGZ KHVBSHVBGLBCSOKBJRBCLCBCYOUBABFigure 3.6: Example of waveform matching in moment tensor inversion for SVItemplate 005. a) Modelled data. Synthetics (red) are superimposed on data (black)for those channels selected for fitting; zero amplitude synthetic traces correspond tounused channels. S-waves are fit to horizontal components and P-waves to the verticalcomponent with shifts up to ±0.2 s applied to maximize correlation. b) Syntheticseismograms for full station complement corresponding to double couple solution de-rived from selected channels in a). Timing misalignments result from errors in velocitymodel used to generate synthetics.Chapter 3. A comparative study of LFE templates in northern Cascadia 54 40’ 20’ 12’ 24’ 36’ 48’ 20’ 49°N124°W 123°WSVIFigure 3.7: Map of double couple mechanisms determined from moment tensorinversion of a selection of LFE templates from the SVI subarray. Inset shows bestdouble couple mechanism from average over all individual moment tensors with arrowindicating corresponding averaged plate motion direction [McCaffrey et al., 2007].Chapter 3. A comparative study of LFE templates in northern Cascadia 55 20’ 40’ 20’ 40’ 50’ 10’ 20’ 30’ 50°N127°W 126°WNVIFigure 3.8: Map of double couple mechanisms for NVI subarray; see caption offigure 3.7 for explanation.Chapter 3. A comparative study of LFE templates in northern Cascadia 56 20’ 40’ 20’ 40’ 36’ 48’ 12’ 24’ 36’ 47°N124°W 123°W 122°WSWFigure 3.9: Map of double couple mechanisms for SW subarray; see caption ofFigure 3.7 for explanation.Chapter 3. A comparative study of LFE templates in northern Cascadia 57 40’ 20’ 40’ 20’ 15’ 30’ 45’ 15’ 48°N124°W 123°W 122°WNWFigure 3.10: Map of double couple mechanisms for NW subarray; see caption offigure 3.7 for explanation.Chapter 4. Tidal modulation of LFEs in northern Cascadia 58Chapter 4Tidal modulation of lowfrequency earthquakes andtriggering of secondary events innorthern Cascadia4.1 IntroductionIt has become increasingly evident in recent years that slow slip, tectonic tremorand LFEs are sensitive to small stress changes. Perturbations at levels of several 10’sof kPa, as generated by surface wave trains from teleseismic events, are capable oftriggering tremor [Miyazawa and Mori , 2005, Peng and Chao, 2008, Peng et al., 2008,Rubinstein et al., 2007, 2009] whereas tidal stresses at levels <10 kPa have been shownto modulate slip [Hawthorne and Rubin, 2010, 2013], tectonic tremor [Gallego et al.,2013, Klaus, 2012, Lambert et al., 2009, Nakata et al., 2008, Rubinstein et al., 2008,Thomas et al., 2013] and LFEs [Shelly et al., 2007, Thomas et al., 2012] during periodsof tremor and associated slip. The analysis of tidal modulation, in particular, has ledto important insights into the mechanism underlying the phenomena.The first report of tidal modulation of tectonic tremor in northern Cascadia wasmade by Rubinstein et al. [2008] who identified clear pulsing of tremor activity withChapter 4. Tidal modulation of LFEs in northern Cascadia 59periods of 12.4 and 24 to 25 hours, corresponding to the principal lunar (M2) and luniso-lar tides. A study of tremor below southern Vancouver Island by Lambert et al. [2009]demonstrated that peak tremor activity occurs at times of maximum tidal shear stressin the thrust direction suggesting that tidal tremor and slip are collocated. Similarmodulation of tremor by tidal shear stress has been observed in southwest Japan be-neath Shikoku [Ide, 2010, Nakata et al., 2008] and northern Washington [Klaus, 2012].To better understand the influence of tidal stresses on slow slip, Hawthorne and Rubin[2010] analyzed strainmeter records from northern Washington and southern Vancou-ver Island for episodic tremor and slip (ETS) events between 2007 and 2009. Theyobserved tidal modulation of slow slip-induced strain at the 12.4h M2 period com-patible with peak strain rate occurring somewhere between time of maximum shearstressing rate and maximum shear stress, depending on slip location. Hawthorne andRubin [2013] used a velocity weakening-to-strengthening friction law to model slow slipevents in the presence of periodic shear stress to simulate the tides. Their simulationsproduced quasi-sinusoidal tidal modulation of the slip rate, with the maximum mo-ment rate occurring close to the maximum applied stress.Thomas et al. [2013] notedthat large, high-amplitude rapid tremor reversals (RTRs, Houston et al. [2011]) andtremor streaks in northern Washington occur almost exclusively during times of posi-tive, thrust-encouraging tidal shear stress on the plate interface, leading them to suggestthat tidal stresses could trigger or intensify large, high-amplitude RTRs that propagatein reverse into previously ruptured and weakened zones on the fault. Beeler et al. [2013]used observations of LFE sensitivity to tidal stressing on the San Andreas fault to con-strain fault rheology. They assumed LFEs are generated on small patches that fail at athreshold stress on an otherwise creeping fault plane. These seismic patches are loadedtectonically, directly by the tides and also by time-dependent creep of the surroundingfault. The authors’ analysis dismissed dislocation creep and dislocation glide as toostrong to be mechanisms for fault creep where LFEs occur, although rate dependentfriction is permitted if effective normal stresses are low, implying near-lithostatic porepressures. Observations of high VP /VS ratios in the LFE source region in Cascadia [Au-det et al., 2009] and southwest Japan [Shelly et al., 2006] supply supporting evidencefor near-lithostatic pore pressures in the subduction zone setting.Chapter 4. Tidal modulation of LFEs in northern Cascadia 60The aforementioned studies clearly establish that slow slip and tremor are modu-lated by tidal forcing. In this paper, we analyze the sensitivity of LFEs to tidal forcesusing three LFE detection catalogues for northern Cascadia, beneath southern Vancou-ver Island and Washington state. LFE catalogues determined using network correla-tion hold advantages over tremor in that the LFE source locations and their temporalrecurrence are more precisely defined, enabling detailed studies of spatio-temporal vari-ability on discrete fault patches [Shelly et al., 2007, Thomas et al., 2012]. Followingthe methodology proposed by Thomas et al. [2012], we quantify the correlation of ourLFE families with tidal stresses, analyze the phase of LFE failure time relative to tidalload, compute the sensitivity to tidal shear stress as a function of time into an ETSevent and the sensitivity of LFEs near the leading edge of slip versus that of RTRs. Wecompare results with previous work in Cascadia, Parkfield and Shikoku and discuss theimplications of these observations for plate boundary properties.4.2 DataWe employ three distinct LFE detection catalogues in northern Cascadia. Thefirst LFE catalogue, assembled by Bostock et al. [2012] for southern Vancouver Island(SVI), was culled to 93 geographically independent LFE families with a total of 148,847detections from ETS episodes between 2003-2013. The two other LFE catalogues wereassembled by Royer and Bostock [2013] for Washington state (see Chapter 3). Thenorthern Washington (NW) catalogue was regrouped into 100 geographically indepen-dent LFE families and includes a total of 60,527 detections from ETS episodes between2007-2011. LFE locations for this catalogue span the eastern flanks of the OlympicPeninsula to the western reaches of Puget Sound. The southern Washington (SW) cat-alogue, regrouped into 49 geographically independent LFE families, includes a total of22,868 detections from ETS episodes during the period 2007-2008, from a region to thesouthwest of Puget Sound. Figure 4.1 shows LFE family locations for the three differentcatalogues. LFE detection methodology, location procedure, depth profiles and focalmechanisms for the 3 catalogues are presented in Bostock et al. [2012] and in Chapter3 [Royer and Bostock , 2013] .Chapter 4. Tidal modulation of LFEs in northern Cascadia 61−125˚00' −124˚30' −124˚00' −123˚30' −123˚00' −122˚30' −122˚00'46˚30'47˚00'47˚30'48˚00'48˚30'49˚00'4912324210 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240Family IDSVI dataset NW dataset SW datasetFigure 4.1: Map of LFE family locations in southern Vancouver Island (SVI detectioncatalogue) and Washington State (NW and SW detection catalogues) assembled byBostock et al. [2012] and Royer and Bostock [2013] (see Chapter 3), respectively. Cyanlines indicate the 20, 30 and 40 km depth contours to the top of the subducting Juande Fuca plate modeled by Audet et al. [2010]. Color bar legend indicates LFE familyIDs ordered by latitude, from north to south and the subdivision of the 3 detectioncatalogues is shown by black vertical lines. The three numbered families are used inplots of figure 4.2.Chapter 4. Tidal modulation of LFEs in northern Cascadia 624.3 Methodology4.3.1 Tidal modelWe compute time series of tidal stress using SPOTL [Agnew , 2012]. The currentlydistributed version of SPOTL calculates horizontal and shear strains from the solidEarth and load tides at the Earth’s surface. To compute strain at depth on a dippingfault requires two modifications. First, we use the Boussinesq solution for a pointload on a half space to compute depth-dependent Green’s functions for the load tidecalculations. These Green’s functions are calculated using a rigidity of µ =30 GPa andPoisson’s ratio, ν=0.25 [Farrell , 1972, Melchior , 1978]. Second, we assume that thestrains from the body tides do not vary substantially over the depth range of interestand are approximately equal to those at the surface. Assuming that shear stresses arezero at the Earth’s surface, the vertical strain can be expressed in terms of the twohorizontal strains as:err =−ν1− ν [eθθ + eλλ]. (4.1)We compute tidally induced fault-normal and up-dip shear stresses, FNS and UDSS,and their time derivatives dFNS and dUDSS, using stress tensor time series and esti-mates of dip and dip azimuth of the subducting Juan de Fuca plate based on a smoothedversion of the model of Audet et al. [2010]. Positive FNS indicates extension and posi-tive UDSS refers to thrust-promoting shear stress in the up-dip direction. These 4 tidaltime series are computed at sample intervals of 15 minutes for each tremor episode inthe period 2003-2013. Figure 4.2 shows an example of tidal time series FNS and UDSSfor LFE families 49 (SVI), 123 (NW) and 242 (SW) (see figure 4.1 for locations). Figure4.3 maps the time-averaged envelope amplitude of stress of each LFE family for the 4tidal components, computed as the magnitude of the analytic tidal signal over the cor-responding detection period. Envelope amplitudes for UDSS and dUDSS increase fromsouth to north, ranging from 481 Pa to 2276 Pa for UDSS, and between 2.6 Pa/minand 17 Pa/min for dUDSS. In contrast, envelope amplitudes for FNS and dFNS arelarger beneath the Strait of Juan de Fuca, where the ocean load tides have a strongerChapter 4. Tidal modulation of LFEs in northern Cascadia 63−2−1012Tidally induced stress (kPa)0 2 4 6 8 10 12 14 16 18 20 222007 tremor episode (days from 16/01/2007)UDSSFNSfamily SW−ID242R = −0.99−7−6−5−4−3−2−101234567Tidally induced stress (kPa)0 2 4 6 8 10 12 14 16 18 20 22UDSSFNSfamily NW−ID123R = −0.79−10−8−6−4−20246810Tidally induced stress (kPa)0 2 4 6 8 10 12 14 16 18 20 22UDSSFNSfamily SVI−ID49R = 0.00Figure 4.2: 2007 tremor episode time series of tidally-induced UDSS (red) and FNS(green) for LFE family 49 from the SVI detection catalogue (upper panel), LFE family123 from the NW detection catalogue (middle panel) and LFE family 242 from theSW detection catalogue (bottom panel). The correlation coefficient R between FNSand UDSS for the given family is indicated in each panel.Chapter 4. Tidal modulation of LFEs in northern Cascadia 64influence. The amplitudes vary between 95 Pa and 5293 Pa for FNS, and between 0.42Pa/min and 26 Pa/min for dFNS.For 50% of the LFE locations, FNS and UDSS time series are anticorrelated with acoefficient correlation lower than -0.5. We may explain this negative correlation betweenFNS and UDSS by considering the relative contributions from different tidal compo-nents. The ocean load tides exert a stronger influence on the stresses in the tremor zoneof northern Cascadia than the body tides by a factor that ranges from 2 to 10. Thisload tide influence can be divided into 3 components: the long wavelength load on theoceanic plate seaward of the trench, loading on the continental slope and shelf, and theload contribution from the inland seas, namely Georgia Strait and Puget Sound. Thefirst contribution is the least important on account of its distance (the trench lies some100 km southwest of the coast). As explained by Lambert et al. [2009], the two lattercontributions are nearly opposite in phase and so contribute constructively to UDSS anddestructively to FNS on that portion of the plate boundary where ETS occurs (roughlymidway between the two water bodies). The offshore coastal load produces oppositelysigned UDSS and FNS on the plate boundary and outweighs the contribution from theinland seas, with the result that FNS and UDSS are anticorrelated. In contrast, LFElocations beneath the Strait of Juan de Fuca (roughly 25% of the total) display thelargest FNS envelope due to the additional load contribution from the Strait of Juande Fuca resulting in weaker correlation between FNS and UDSS.4.3.2 LFE correlation with tidal stressFor each LFE family location, we compute the tidally induced stresses and theirderivatives at times corresponding to detections and compute the expected number ofdetections based on the amount of time spent in a particular loading condition (firstrelative to the sign of tidally induced stresses and stress rates, and later to the stress-magnitude range), assuming that LFE detections occur randomly in a given tremorepisode. We quantify the degree of correlation of LFE occurrence with tidal stress bycomputing the excess value:Nex =# observed detections−# expected detections# expected detections. (4.2)Chapter 4. Tidal modulation of LFEs in northern Cascadia 65dFNSC0 5 10 15 20 25 30Mean envelope amplitude (Pa/min)dUDSSDFNSA400 800 1200 1600 2000 2400 2800 3200 3600 4000 4400 4800 5200Mean envelope amplitude (Pa)UDSSBFigure 4.3: Map of tidal stress amplitudes. Mean tidal envelope amplitudes overdetection time intervals are plotted color coded at locations of LFE families: FNS(panel A), UDSS (panel B), dFNS (panel C) and dUDSS(panel D) time series.Chapter 4. Tidal modulation of LFEs in northern Cascadia 66Positive Nex values indicate a surplus of detections in the loading condition consid-ered, and negative Nex values indicate a deficit of detections. We include 95% and99% confidence intervals by generating 25,000 random catalogues of N detections each(the number N varies with the LFE family considered, between 386 and 3575 for SVI,between 141 and 1942 for NW, and between 204 and 895 for SW) to compute a distribu-tion of Nex values. Nex values that fall outside the 99% confidence interval are less than1% likely to occur randomly assuming that tides and LFEs are uncorrelated. Figure 4.4shows a map of LFE family locations, color coded by Nex values corresponding to the 4tidal components. We plot in figure 4.5 all possible combinations of Nex values againsteach other. We analyze correlation between Nex values and the mean stress envelopeamplitudes by plotting them against each other for each tidal component, as shown infigure 4.6.We consider how specific stress levels influence LFE occurrence by computing LFErates as a function of the magnitude of the applied stress components. We determinethe observed number of LFE detections for a particular family that falls within a givenstress range (40 bins over the total stress range, set independently for each of FNS,UDSS, dFNS and dUDSS time series). The expected number of detections in each binis then computed based on the null hypothesis that LFE origin times and tides areuncorrelated. The ratio# observed detections# expected detections= Nex + 1, (4.3)is computed for each bin for a given stress component and for all families. Ratios greaterthan one translate to a surplus of detections and ratios lower than 1 indicate a deficit ofdetections compared to a random detection distribution. Results are plotted in figure4.7.We assign an effective phase angle to each LFE detection of a particular familyusing a method similar to that of Tanaka et al. [2002]. For each of the two pairs ofcomponents UDSS, dUDSS and FNS, dFNS, the maximum and the minimum tidalstresses are assigned to 0◦ and ±180◦, respectively. However, in contrast to the methodof Tanaka et al. [2002], the maximum and minimum tidal stress rates are assignedChapter 4. Tidal modulation of LFEs in northern Cascadia 67dFNS99% Cl SVI [−17.0;16.8]%99% Cl NW [−28.2;26.9]%99% Cl SW [−24.1;23.6]%C−50 −40 −30 −20 −10 0 10 20 30 40 50Nex (%)dUDSS99% Cl SVI [−13.0;13.0]%99% Cl NW [−20.7;20.2]%99% Cl SW [−18.8;18.7]%DFNS99% Cl SVI [−20.0;19.1]%99% Cl NW [−31.1;32.0]%99% Cl SW [−28.1;27.3]%A5567187UDSS99% Cl SVI [−17.2;16.5]%99% Cl NW [−28.5;30.8]%99% Cl SW [−22.4;21.4]%BB A AB46˚30'47˚00'47˚30'48˚00'48˚30'49˚00'46˚30'47˚00'47˚30'48˚00'48˚30'49˚00'−125˚ −124˚ −123˚ −122˚ −125˚ −124˚ −123˚ −122˚Figure 4.4: Map of Nex values relative to the positive sign of tidally induced stressesand stress rates. LFEs are color-coded by Nex values corresponding to the tidal FNS(panel A), UDSS (panel B), dFNS (panel C) and dUDSS(panel D) components. 99%confidence intervals (CI) calculated for the family with the fewest detections are givenin the bottom left of each panel for each detection catalogue and significant familiesare indicated with black rims. The three numbered families in panel B are used in therecurrence analysis (section 4.5.3).Chapter 4. Tidal modulation of LFEs in northern Cascadia 68−60−50−40−30−20−10010203040−20 −10 0 10 20 30 40 50 60SVI R2=0.09NW R2=0.01SW R2=0.61−60−50−40−30−20−10010203040−60−50−40−30−20−10 0 10 20 30 40 50SVI R2=0.24NW R2=0SW R2=0−60−50−40−30−20−10010203040−20 −10 0 10 20 30 40 50SVI R2=0.01NW R2=0.16SW R2=0.35−20−100102030405060−60−50−40−30−20−10 0 10 20 30 40 50SVI R2=0.25NW R2=0.1SW R2=0.06−20−100102030405060−20 −10 0 10 20 30 40 50SVI R2=0.01NW R2=0.01SW R2=0.09−60−50−40−30−20−1001020304050−20 −10 0 10 20 30 40 50SVI R2=0NW R2=0.16SW R2=0.31UDSS Nex (%) dFNS Nex (%) dUDSS Nex (%)dFNS Nex (%)UDSS Nex (%)FNS Nex (%)Figure 4.5: Nex values relative to the positive sign of tidally induced stresses andstress rates are plotted against each other, for SVI (blue), NW (green) and SW (red)detection catalogues. A coefficient of determination R2 (in this case the square of thesample correlation coefficient) is provided in each panel for the 3 detection catalogues,distinguished by different colors (blue - SVI, green - NW, red - SW). High values ofR2 indicate a strong correlation or anticorrelation between two components.to 90◦ and -90◦, respectively, increasing the accuracy of our assigned effective phases.The effective phase is assigned by interpolating the time of the event between timeswith assigned phases. We assemble the observed numbers of detections that fall withina phase angle range (here we consider 36 bins of 10◦), and the expected numbers ofdetections in each bin are computed based on the null hypothesis that LFEs and tidesare randomly correlated. We compute again the ratio Nex + 1 for each bin of a givencomponent pair (FNS,dFNS and UDSS,dUDSS) and for all families. Rose histograms ofphases for the 3 different LFE detection catalogues are plotted in figure 4.8. A schematicphase plot with labeling is provided to aid in understanding the rose histograms. BinsChapter 4. Tidal modulation of LFEs in northern Cascadia 69−30−20−10010203040FNS Nex (%)1000 2000 3000 4000 5000FNS mean envelope amplitude (Pa)SVI R2=0NW R2=0.17SW R2=0.04−30−20−10010203040UDSS Nex (%)600 900 1200 1500 1800 2100UDSS mean envelope amplitude (Pa)SVI R2=0.04NW R2=0.03SW R2=0.01−40−30−20−1001020304050dFNS Nex (%)0 5 10 15 20 25 30dFNS mean envelope amplitude (Pa/min)SVI R2=0.06NW R2=0.02SW R2=0.05−30−20−1001020304050dUDSS Nex (%)0 5 10 15 20dUDSS mean envelope amplitude (Pa/min)SVI R2=0.13NW R2=0SW R2=0.05Figure 4.6: Nex values relative to the positive sign of tidally induced stresses andstress rates plotted against the mean envelope amplitude, for each tidal componentFNS (top -left), UDSS (top-right), dFNS (bottom-left) and dUDSS (bottom-right)and for detection catalogues SVI (blue), NW (green) and SW (red). A coefficientof determination R2 (in this case the square of the sample correlation coefficient) isprovided in each panel for the 3 detection catalogues, distinguished by different colors(blue - SVI, green - NW, red - SW).that contain the expected number of detections have a radius of 1 (100% of the expectedvalues assuming a random distribution shown by red contour). Confidence intervalscorresponding to the 95% and 99% levels are also provided on each polar histogram(blue contours).Chapter4.TidalmodulationofLFEsinnorthernCascadia70−25−20−15−10−5051015202530UDSS rate (Pa/min)0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240Family
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Studies of seismic deconvolution and low-frequency...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Studies of seismic deconvolution and low-frequency earthquakes Royer, Alexandra Amélie 2014
pdf
Page Metadata
Item Metadata
Title | Studies of seismic deconvolution and low-frequency earthquakes |
Creator |
Royer, Alexandra Amélie |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | Deconvolution of seismic data is an important component of signal processing that aims to remove the seismic source from seismograms, thereby isolating the Green’s function. By considering seismograms of multiple earthquakes from similar locations recorded at a given station and that therefore share the same Green’s function, we investigate a system of equations where the unknowns are the sources and source durations. Our solution is derived using direct linear inversion to recover the sources and Newton’s method to recover source durations. For the short seismogram durations considered, we are able to recover source time functions for noise levels at 1% of the direct P -wave amplitude. However, the nonlinearity of the problem renders the system expensive to solve and sensitive to noise; therefore consideration is limited to short seismograms with high signal-to-noise ratio (SNR). When SNR levels are low, but a large multiplicity of seismograms representing a common source-receiver path are available, we can apply a different deconvolution approach to recover the Green’s function. In an application to tectonic tremor in northern Cascadia, we implement an iterative blind deconvolution method that involves correlation, threshold detection and stacking of 1000’s of low frequency earthquakes (LFEs) that form part of tremor to generate templates that can be considered as empirical Green’s functions. We exploit this identification to compute hypocentres and moment tensors. LFE hypocentres follow the general epicentral distribution of tremor and occur along tightly defined surfaces in depth. The majority of mechanisms are consistent with shallow thrusting in the direction of plate motion. We analyze the influence of ocean tides on the triggering of LFEs and find a spatially variable sensitivity to tidally induced up-dip shear stress (UDSS), suggesting that tidal sensitivity must partially depend on laterally heterogeneous physical properties. The majority of LFEs fail during positive and increasing UDSS, consistent with combined contributions from background slow slip and from tides acting directly on LFEs. We identified rapid tremor reversals in southern Vancouver Island with higher sensitivity to UDSS than the main front and which at least partially explains an observed increase in LFE sensitivity to UDSS with time. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-08-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0165914 |
URI | http://hdl.handle.net/2429/49978 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2014-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2014_september_royer_alexandra.pdf [ 20.35MB ]
- Metadata
- JSON: 24-1.0165914.json
- JSON-LD: 24-1.0165914-ld.json
- RDF/XML (Pretty): 24-1.0165914-rdf.xml
- RDF/JSON: 24-1.0165914-rdf.json
- Turtle: 24-1.0165914-turtle.txt
- N-Triples: 24-1.0165914-rdf-ntriples.txt
- Original Record: 24-1.0165914-source.json
- Full Text
- 24-1.0165914-fulltext.txt
- Citation
- 24-1.0165914.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0165914/manifest