Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Drolet, Doriane
Doctor of Philosophy - PhD
University of British Columbia
Large earthquakes often dominate attention due to their significant impact, but the majority of seismic activity consists of smaller and more frequent earthquakes. These smaller earthquakes, when properly analyzed, can produce enormous datasets and greatly enhance ourunderstanding of Earth’s dynamic processes. Through three research projects, this thesis focuses on developing and applying methods to analyze small earthquakes (M < 4) using the relative moment tensor inversion, addressing the challenges of small seismic sources. The first project examines two significant M7.1 intraslab earthquakes in southern Alaska: the 2016 Iniskin (24/01/2016) and the 2018 Anchorage (30/11/2018) earthquakes. The study provides insights into rupture processes and stress evolution by relocating aftershocks, recovering moment tensors, and inverting for stress variation. The aftershock distribution of the Iniskin event suggests a rupture along a steeply dipping fault in the Pacific Plate, with evidence of pre-existing fault reactivation (likely an outer-rise fault) and fluid migration. The Anchorage sequence reveals complex aftershock patterns, likely due to simultaneous rupture on both nodal planes. The second project tackles the challenges of determining moment tensors for small earthquakes recorded by sparse regional networks. Building upon previous work, a simultaneous multiple-event relative moment tensor inversion scheme is developed, leveraging constraints from P-wave polarities, relative P- and S-wave amplitudes, and local magnitude estimates without requiring prior knowledge of a reference moment tensor. This method is validated with synthetic examples and applied to a seismic cluster near San Juan Island, Washington, revealing insights into the local stress regime and fluid migration from the subducting Juan de Fuca Plate. The final project focuses on the Haida Gwaii incipient subduction zone in Canada. The study uses the novel relative moment tensor inversion technique to analyze earthquake groups in this sparsely instrumented area. Moment tensor results for events near the expected Queen Charlotte Fault and thrust zone display reverse mechanisms due to underthrusting of the Pacific Plate. Overall, this thesis advances the understanding of seismic source processes and stress distributions across different tectonic settings, offering valuable insights into fault mechanics, role of geologic structure, fluid presence and regional seismotectonics through innovative inversion methodologies.
Advancements in Seismotectonic Analysis through Relative MomentTensor InversionbyDoriane DroletBaccalaure´at en physique, Universite´ de Sherbrooke, 2018A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Geophysics)The University of British Columbia(Vancouver)December 2024© Doriane Drolet, 2024The following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the thesis entitled:Advancements in Seismotectonic Analysis through Relative Moment Tensor Inversionsubmitted by Doriane Drolet in partial fulfillment of the requirements for the degree of Doctor ofPhilosophy in Geophysics.Examining Committee:Michael Bostock, Professor, Earth, Ocean and Atmospheric Sciences, UBCSupervisorSimon Peacock, Professor, Earth, Ocean and Atmospheric Sciences, UBCSupervisory Committee MemberCharles Sammis, Professor Emeriti, Earth, Ocean and Atmospheric Sciences, UBCSupervisory Committee MemberChristian Schoof, Professor, Earth, Ocean and Atmospheric Sciences, UBCUniversity ExaminerBratt Gladman, Professor, Physics and Astronomy, UBCUniversity ExaminerAbhijit Ghosh, Professor, Earth and Planetary Sciences, University of California, RiversideExternal ExamineriiAbstractLarge earthquakes often dominate attention due to their significant impact, but the majority of seis-mic activity consists of smaller and more frequent earthquakes. These smaller earthquakes, whenproperly analyzed, can produce enormous datasets and greatly enhance our understanding of Earth’sdynamic processes. Through three research projects, this thesis focuses on developing and applyingmethods to analyze small earthquakes (M < 4) using the relative moment tensor inversion, address-ing the challenges of small seismic sources.The first project examines two significant M7.1 intraslab earthquakes in southern Alaska: the2016 Iniskin (24/01/2016) and the 2018 Anchorage (30/11/2018) earthquakes. The study providesinsights into rupture processes and stress evolution by relocating aftershocks, recovering momenttensors, and inverting for stress variation. The aftershock distribution of the Iniskin event suggestsa rupture along a steeply dipping fault in the Pacific Plate, with evidence of pre-existing fault reac-tivation (likely an outer-rise fault) and fluid migration. The Anchorage sequence reveals complexaftershock patterns, likely due to simultaneous rupture on both nodal planes.The second project tackles the challenges of determining moment tensors for small earthquakesrecorded by sparse regional networks. Building upon previous work, a simultaneous multiple-eventrelative moment tensor inversion scheme is developed, leveraging constraints from P-wave polar-ities, relative P- and S-wave amplitudes, and local magnitude estimates without requiring priorknowledge of a reference moment tensor. This method is validated with synthetic examples and ap-plied to a seismic cluster near San Juan Island, Washington, revealing insights into the local stressregime and fluid migration from the subducting Juan de Fuca Plate.The final project focuses on the Haida Gwaii incipient subduction zone in Canada. The studyuses the novel relative moment tensor inversion technique to analyze earthquake groups in thissparsely instrumented area. Moment tensor results for events near the expected Queen CharlotteFault and thrust zone display reverse mechanisms due to underthrusting of the Pacific Plate.Overall, this thesis advances the understanding of seismic source processes and stress distribu-tions across different tectonic settings, offering valuable insights into fault mechanics, role of ge-ologic structure, fluid presence and regional seismotectonics through innovative inversion method-ologies.iiiLay SummaryWhile large earthquakes grab attention, small earthquakes occur more frequently and provide valu-able data when analyzed effectively. Through this thesis, I offer a detailed analysis of seismicsequences composed of earthquakes with magnitude smaller than 4 across western North America(Alaska (USA), British-Columbia (Canada) and Washington (USA)). Each sequence is examinedusing an advanced technique to recover the radiation pattern of individual earthquakes part of adense earthquake group based on nearby seismic wave recordings (seismograms). This methodhas been improved and developed to accurately determine the mechanism of small earthquakes,even when data is limited. Additionally, this study explores the role of fluids and regional forceswhich influence earthquake behaviour. Overall, the findings contribute to our knowledge of howearthquakes occur and evolve, offering new perspectives on the dynamics of seismic events acrossdifferent tectonic settings.ivPrefaceThis thesis represents the research conducted by Doriane Drolet under the guidance and supervisionof Professor Michael Bostock. It was written by myself and revised by Michael Bostock. The workpresented herein is composed of three research chapters, two has been accepted to peer-reviewedjournals and one that I consider submitting to a peer-review journal.Chapter 2 is a reproduction of the published paper :Drolet, D., M. G. Bostock, A. P. Plourde, and C. G. Sammis, Aftershock distributions, moment ten-sors and stress evolution of the 2016 Iniskin and 2018 Anchorage Mw 7.1 Alaskan intraslab earth-quakes, Geophysical Journal International, 231(1), 199–214, DOI: 10.1093/gji/ggac165, 2022.The research scope was developed by Michael Bostock and myself. I was responsible for thedata acquisition, processing, analysis and first draft. Alex Plourde provided essential code and sup-port for the data processing. The interpretation of the results was a collaborative effort, particularlybetween Michael Bostock and myself. Charles Sammis provided a valuable discussion on fluid dif-fusion after the Iniskin Earthquake.Chapter 3 is a reproduction of the published paper :Drolet, D., Bostock, M., and Peacock, S., Relative Moment Tensor Inversion for Microseismicity:Application to Clustered Earthquakes in the Cascadia Forearc. Seismica, 2(4). DOI: 10.26443/seis-mica.v2i4.1311, 2024.This paper extended work done by Alex Plourde on relative moment tensor inversion. Theproject was designed collaboratively by Michael Bostock and myself. I was responsible for writingthe first manuscript draft. Simon Peacock helped in the interpretation of real data. Python code writ-ten by myself of the developed algorithm is available on Zenodo DOI: 10.5281/zenodo.10850729as well as other supplementary documents. Wasja Bloch will further develop this tool as part of theMarie-Sklodowska-Curie Action funded by the European Union.Chapter 4’s project was designed by Michael Bostock and myself. I was responsible for allaspects of data acquisition and processing. The interpretation of the results was a joint effort amongvthe Michael Bostock, Simon Peacock and myself. This chapter was written by myself and revisedby Michael Bostock.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Objectives and Research Questions . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Background Knowledge and Methodology Overview . . . . . . . . . . . . . . . . 31.2.1 Subduction Zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.2 Seismic Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.3 Moment Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.4 Focal Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.5 Beachball Diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2.6 Moment Tensor Inversion Methods . . . . . . . . . . . . . . . . . . . . . 71.2.7 Stress Inversion Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Aftershock Distributions, Moment Tensors and Stress Evolution of the 2016 Iniskin& 2018 Anchorage Mw 7.1 Alaskan Intraslab Earthquakes . . . . . . . . . . . . . . 102.1 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10vii2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3.1 January 24, 2016, Iniskin Event . . . . . . . . . . . . . . . . . . . . . . . 122.3.2 November 30, 2018, Anchorage Event . . . . . . . . . . . . . . . . . . . . 152.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.4.1 Earthquake Relocations . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.4.2 Moment Tensor Inversions . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4.3 Temporal Stress Inversions . . . . . . . . . . . . . . . . . . . . . . . . . . 182.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5.1 Hypocenteral Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5.2 Moment tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.5.3 Stress Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.6.1 Iniskin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.6.2 Anchorage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343 Relative Moment Tensor Inversion for Microseismicity: Application to ClusteredEarthquakes in the Cascadia Forearc . . . . . . . . . . . . . . . . . . . . . . . . . . 353.1 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3.2 Relative Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.3.3 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.4 Synthetic Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.5 Real Case Example: San Juan Cluster . . . . . . . . . . . . . . . . . . . . . . . . 503.5.1 Data Set and Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.5.2 Stress Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604 Moment Tensor Analysis in the Southern Haida Gwaii: Insights into UnderthrustingProcesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.1 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.3 Data and Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69viii4.4.1 Moment Tensor Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 694.4.2 Stress Inversion Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.5.1 Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.5.2 Slip Partitioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.5.3 Methodology limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 755 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 765.1 Question 1: What information can we extract from small earthquakes? . . . . . . . 765.2 Question 2: How can we improve or adapt pre-existing methods to recover mecha-nisms for microseismic events given their challenges? . . . . . . . . . . . . . . . . 775.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80A Supporting Materials for Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 95A.1 Data Availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95A.2 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96A.3 Slab E115 vs Slab 2 models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102A.4 Decay of the pressure in a pressurized fault zone . . . . . . . . . . . . . . . . . . . 104A.4.1 Analytic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104A.4.2 Relation between R and P . . . . . . . . . . . . . . . . . . . . . . . . . . 105A.4.3 Relation between the Omori-Utsu law and pressure . . . . . . . . . . . . . 106B Supporting Materials for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 110B.1 Data Availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110B.2 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111B.3 Moment Tensor Results for 350 Earthquakes . . . . . . . . . . . . . . . . . . . . . 113C Supporting Materials for Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 124C.1 Data Availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124C.2 Supplementary Figures and Tables . . . . . . . . . . . . . . . . . . . . . . . . . . 125C.2.1 Misfit for the different velocity models . . . . . . . . . . . . . . . . . . . 125C.2.2 MT solutions for the different velocity models . . . . . . . . . . . . . . . 126C.3 Analysis of the earthquake group L1 . . . . . . . . . . . . . . . . . . . . . . . . . 131C.3.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131C.3.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135C.3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137ixC.3.4 Moment Tensors solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 138xList of TablesTable 2.1 Nodal plane angles associated with different groups or events within the Anchor-age sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Table 4.1 General information about each earthquake group . . . . . . . . . . . . . . . . 66Table 4.2 Properties of the final solutions for each group . . . . . . . . . . . . . . . . . . 70Table A.1 Value of parameters in the diffusion equation. . . . . . . . . . . . . . . . . . . 105Table C.1 Residual and polarity agreement from each group and velocity model . . . . . . 125Table C.2 General information about earthquake group L1 . . . . . . . . . . . . . . . . . 132Table C.3 Residual and polarity agreement with the different velocity models . . . . . . . 133xiList of FiguresFigure 1.1 Subduction zone diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.2 Example of beachball diagrams . . . . . . . . . . . . . . . . . . . . . . . . . 6Figure 1.3 Example of a beachball diagram for full MT . . . . . . . . . . . . . . . . . . . 7Figure 2.1 Overview of the study region . . . . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 2.2 Earthquakes distribution for the 2016 Iniskin Sequence and the 2018 Anchoragesequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Figure 2.3 Relocated events near the Iniskin earthquake . . . . . . . . . . . . . . . . . . 20Figure 2.4 3D representation of relocated events near the Anchorage earthquake . . . . . 21Figure 2.5 Moment tensor clusters for the Iniskin sequence . . . . . . . . . . . . . . . . . 22Figure 2.6 Waveform comparisons with corresponding stations projected onto focal sphere 23Figure 2.7 Moment tensor clusters for the Anchorage sequence . . . . . . . . . . . . . . 24Figure 2.8 Waveform comparisons with corresponding stations projected onto focal sphere 25Figure 2.9 Nodal plane angles of the 108 MTs that form the event plane defined by thenorthern group . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.10 Temporal stress inversion results for the Iniskin sequence . . . . . . . . . . . . 27Figure 2.11 Temporal stress inversion results for the Anchorage sequence . . . . . . . . . . 28Figure 2.12 Stress ratio recovered for different inversion window sizes . . . . . . . . . . . 31Figure 3.1 Schematic diagram summarizing the assumptions . . . . . . . . . . . . . . . . 38Figure 3.2 Summary of the relative moment tensor algorithm . . . . . . . . . . . . . . . . 40Figure 3.3 Summary of the four solving stages . . . . . . . . . . . . . . . . . . . . . . . 44Figure 3.4 Synthetic examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 3.5 Synthetic results for 1,000 realizations . . . . . . . . . . . . . . . . . . . . . . 50Figure 3.6 Seismic stations used for MT inversion of local seismicity near San Juan Island 51Figure 3.7 Temporal distribution of San Juan Island seismicity . . . . . . . . . . . . . . . 52Figure 3.8 Earthquake distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54Figure 3.9 Optimal solution selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55Figure 3.10 Stable and likely MTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56xiiFigure 3.11 Kaverian diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.12 Stress inversion results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.13 Seismicity Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59Figure 4.1 Schamatic of the tectonic setting for oblique convergence . . . . . . . . . . . . 63Figure 4.2 Overview of the study region . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 4.3 Overview of the earthquake groups analyzed in this study . . . . . . . . . . . . 67Figure 4.4 Cross-section showing the earthquake group locations . . . . . . . . . . . . . 71Figure 4.5 Moment tensor solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 4.6 Kaverian diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.7 Stress inversion results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure A.1 Anchorage sequence relocation comparison . . . . . . . . . . . . . . . . . . . 96Figure A.2 Choice of HASH solution weights for Iniskin and Anchorage sequences. . . . 97Figure A.3 Quality test on moment tensors of the Iniskin sequence . . . . . . . . . . . . . 98Figure A.4 Quality test on moment tensors of the Anchorage sequence . . . . . . . . . . . 99Figure A.5 Reference events used for the Iniskin moment tensor inversion. . . . . . . . . . 100Figure A.6 Reference events used for the Anchorage moment tensor inversion. . . . . . . 100Figure A.7 Intermediate deep events of Mw≥6 with moment tensor solutions . . . . . . . 101Figure A.8 Difference between E115 Slab and Slab 2 models . . . . . . . . . . . . . . . . 103Figure A.9 Fault zone geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104Figure A.10 Comparison of analytic and numerical solutions for the pressure distribution . . 107Figure A.11 Cumulative number of aftershocks as a function of time after the mainshock . . 108Figure A.12 Difference between the cumulative number of aftershocks and the Omori-Utsulaw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109Figure B.1 Synthetics tests with extra stations . . . . . . . . . . . . . . . . . . . . . . . . 111Figure B.2 Synthetics tests with extra events . . . . . . . . . . . . . . . . . . . . . . . . . 111Figure B.3 Earthquake locations before relocation . . . . . . . . . . . . . . . . . . . . . . 112Figure B.4 Colour code for MTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113Figure B.5 MT results from 2001/11/29 to 2004/8/13 . . . . . . . . . . . . . . . . . . . . 114Figure B.6 MT results from 2004/8/18 to 2007/4/16 . . . . . . . . . . . . . . . . . . . . . 115Figure B.7 MT results from 2007/4/17 to 2011/1/3 . . . . . . . . . . . . . . . . . . . . . 116Figure B.8 MT results from 2011/2/16 to 2012/7/4 . . . . . . . . . . . . . . . . . . . . . 117Figure B.9 MT results from 2012/7/4 to 2013/9/3 . . . . . . . . . . . . . . . . . . . . . . 118Figure B.10 MT results from 2013/9/11 to 2017/1/30 . . . . . . . . . . . . . . . . . . . . . 119Figure B.11 MT results from 2017/1/30 to 2018/8/20 . . . . . . . . . . . . . . . . . . . . . 120Figure B.12 MT results from 2018/9/23 to 2019/7/23 . . . . . . . . . . . . . . . . . . . . . 121xiiiFigure B.13 MT results from 2019/8/21 to 2020/7/11 . . . . . . . . . . . . . . . . . . . . . 122Figure B.14 Moment tensor solutions for stable events clustered by source type . . . . . . . 123Figure C.1 Legend for Figs C.2 to C.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Figure C.2 MTs results for group O1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127Figure C.3 MTs results for group O2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128Figure C.4 MTs results for group O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129Figure C.5 MTs results for group O3 (continuation) . . . . . . . . . . . . . . . . . . . . . 130Figure C.6 Overview of the earthquake group L1 . . . . . . . . . . . . . . . . . . . . . . 131Figure C.7 Relocation comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132Figure C.8 Moment tensor solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134Figure C.9 Time evolution of the earthquakes . . . . . . . . . . . . . . . . . . . . . . . . 135Figure C.10 Kaverian diagrams and stress inversion results . . . . . . . . . . . . . . . . . . 136Figure C.11 Colour-code legend for the best solution . . . . . . . . . . . . . . . . . . . . . 138Figure C.12 MTs results from 2013/4/26 to 2013/8/23 . . . . . . . . . . . . . . . . . . . . 139Figure C.13 MTs results from 2013/8/26 to 2013/9/13 . . . . . . . . . . . . . . . . . . . . 140Figure C.14 MTs results from 2013/9/14 to 2013/10/30 . . . . . . . . . . . . . . . . . . . 141Figure C.15 MTs results from 2013/11/6 to 2013/11/19 . . . . . . . . . . . . . . . . . . . 142xivList of AbbreviationsCMT Centroid moment tensorsDOI Document Object Identifier (see http://doi.org) (also referred as doi)FM Focal mechanismM MagnitudeML Local magnitudeMW Moment magnitudeM0 Seismic momentMT Moment tensorPCA Principal component analysisR Stress Ratio also known as shape ratio, relative stress magnitude or eigenvalue ratioUBC University of British Columbia1D One-dimensional2D Two-dimensional3D Three-dimensionalxvAcknowledgmentsCompleting this PhD thesis marks a significant milestone in my academic journey, and it would nothave been as rewarding and fulfilling without the support and guidance of numerous individuals. Iwould like to take this opportunity to express my heartfelt gratitude to them.First and foremost, I would like to express my most profound appreciation to my supervisor,Michael Bostock, for his kindness, unwavering support, invaluable guidance and sharing his passionthroughout my research. Your expertise, patience, and constructive feedback have been instrumentalin shaping this work and my development.I extend my sincere thanks to Simon Peacock and Charles Sammis, the members of my thesiscommittee, for their insightful comments, critical evaluations, and suggestions, which significantlyimproved the quality of this thesis.A special mention goes to my fellow researchers and colleagues in the Seismology Group, in-cluding Geena Littel, Raul Mendoza, Reid Merrill, Alexandre Pourde, Alex Hughes and WasjaBloch. Your collaboration and camaraderie have enriched my time at UBC.I also want to thank all the great friends I made during this journey. Your companionship mademy free time truly enjoyable, and I have cherished memories of our outdoor adventures, beach time,dinners, and countless shared moments. Your friendship has been an invaluable source of joy andsupport throughout this experience.I am indebted to my family for their unconditional love, patience, and understanding. To myparents, Marc Drolet and Nathalie Richard, thank you for your unwavering belief in me, for alwaysencouraging me to dream, and for your support. To my partner, Ben, thank you for your support,encouragement and being a source of stress relief. I am deeply grateful for your presence in mylife. I am also thankful to Ben’s family (Jane, Denis, Berlyn and Melissa) for warmly welcomingme into their home and treating me as one of their own.Lastly, I would like to acknowledge the funding agencies and organizations that supported myresearch financially, including Fonds de recherche du Que´bec - Nature et technologies, Four YearFellowships For PhD Students, Eric Walter Mountjoy Memorial Scholarship and Thomas and Mar-guerite MacKay Memorial Scholarship, without which this work would not have been possible.xviChapter 1IntroductionSeismology, the scientific study of earthquakes and the propagation of elastic waves through planets,is a crucial tool for probing the Earth’s interior. It utilizes both natural and artificial seismic sourcesto investigate the Earth’s structure, dynamics, and composition. Among these sources, earthquakes,known as passive sources, provide an abundance of data and allow us to probe processes and struc-tures to greater depths than anthropogenic (active) sources.Earthquakes occur when stored energy is suddenly released, usually as a result of stress ac-cumulation along faults or fluid overpressure. This release generates seismic waves that propagateoutward from the earthquake’s focus and generate ground shaking. By studying these waves, we candetermine the location (space and time), magnitude, rupture mechanisms, and much more, whichin turn provides insights into the local stress and strain distribution. Earthquakes can cause pro-longed and intense ground shaking, leading to fatalities and infrastructure damage, highlighting theimportance of studying these fascinating natural phenomena.Large earthquakes often attract the most attention due to their significant impact and the wealthof data they produce. However, the majority of seismic activity is composed of numerous smallerearthquakes, which occur much more frequently but are often overlooked. If properly analyzed,these smaller events have the potential to produce enormous datasets and greatly enhance our un-derstanding of the Earth’s dynamic processes.Despite their lower magnitudes (M), small earthquakes carry crucial information about stressdistribution, fault networks, geologic features, and small-scale deformation. They can help fill thespatial and temporal gaps between larger seismic events, offering a more continuous and detailedpicture of seismotectonic activity. They also provide insights into the processes that lead to largerearthquakes and aftershock sequences that follow and improve imaging of the Earth’s interior.However, extracting meaningful information from small earthquakes poses several challenges.Firstly, the low signal-to-noise ratio in the seismic data for these events makes it difficult to deter-mine their key attributes accurately. Small earthquakes produce weaker signals that can be easily1masked by background noise or other seismic activity. Hence analysis of small earthquakes re-quires the development of advanced techniques to enhance the signal quality and extract reliableinformation from the available data.Secondly, their low magnitudes generate low levels of ground shaking that is only detectableat nearby seismic stations. Thus, the sparse distribution of seismic stations can limit the amountof data available for small events. In many regions, especially remote or less instrumented areasseismic network sampling is not dense enough to provide comprehensive coverage, complicatingthe analysis.Thirdly, as a result of reduced fault dimensions, small earthquakes release their dominant energythrough high-frequency waves that are more prone to scattering as a result of Earth’s heterogeneityspectrum. Consequently, waveforms are more complex in addition to suffering higher levels of at-tenuation, rendering them more difficult to interpret. Advanced modelling and inversion techniquesare therefore required to account for these complexities and to retrieve accurate information aboutseismic sources.1.1 Objectives and Research QuestionsGiven these challenges, this thesis focuses on developing methods, analyzing and interpreting re-sults for small earthquakes (M < 4) and microseismic events (unfelt seismic events). By employingnewly developed methods of relative moment tensor inversion, I aim to overcome the difficulties as-sociated with low signal-to-noise ratios, sparse station distribution, and complex wave propagation.Through this research, I seek to unlock invaluable information contained in small seismic events andenhance our understanding of the underlying processes that govern seismic activity. Ultimately, thiswork aims to contribute to the creation of more comprehensive and detailed seismotectonic models,providing insights that are crucial for both scientific research and earthquake hazard assessment.Here are the two main research questions that I attempt to answer in this thesis, including severalsubsidiary queries:1. What information can we extract from small earthquakes ?(a) What do their spatial or temporal variations tell us?As a result of earthquake power-law size-frequency scaling, there are many more smallearthquakes than large ones. Thus, small earthquakes afford larger datasets, yieldinggreater spatial and temporal characterization of faulting and faulting processes.(b) Can we get insight on earthquake triggering mechanisms?This includes exploring the role of different geologic structures in the triggering process.Additionally, fluids are often implicated as a key factor in triggering different types ofearthquakes. By analyzing small earthquakes, we may uncover valuable insights intothese processes.22. How can we improve or adapt pre-existing methods to recover mechanisms for microseismicevents given the previously mentioned challenges?(a) How can we maximize constraints on source mechanisms?Due to the signal-to-noise challenges of characterizing small earthquakes, fewer datapoints per event are available, so including new constraints could improve the accuracyof numerical moment tensor inversions.(b) Can this be done computationally efficiently?Adding new, non-linear constraints will increase computational complexity. Thus, look-ing at ways to speed up numerical processes are important.Question 1 is addressed throughout the thesis, whereas Question 2 will be mainly answered inChapter 3 and briefly touched upon in Chapter 4.1.2 Background Knowledge and Methodology OverviewBefore addressing my research questions, I will introduce some important terms, concepts and meth-ods.1.2.1 Subduction ZoneAll study areas in this thesis concern subduction zones: the Aleutian subduction zone in southernAlaska (Chapter 2), the Cascadia subduction zone on the western Canada-United-States border(Chapter 3) and the nascent Haida Gwaii subduction zone in western Canada (Chapter 4).A subduction zone is a region where two tectonic plates converge, and one plate is forced be-neath the other into the mantle (see Fig. 1.1 for a diagram). This process typically involves anoceanic plate being driven under a continental plate or another oceanic plate. Subduction zonesare characterized by intense geological activity, including the formation of deep oceanic trenches,volcanic arcs, and frequent, often powerful, earthquakes. The latter are a primary motivation for myresearch.The deepest earthquakes on the planet occur within subduction zones and can reach depths asgreat as ∼ 700 km. As a plate subducts, pressures and temperatures increase, leading to the re-lease of fluids through metamorphic dehydration reactions, which can participate in seismic events.Understanding subduction zones is essential for assessing major seismic hazards and elucidatingEarth’s dynamic processes.1.2.2 Seismic WavesAn earthquake generates several types of seismic waves, each with distinct properties and be-haviours as they travel through the Earth. These waves can be categorized into two main groups:3Figure 1.1: Subduction zone diagram from U.S. Geological Survey [2016]body waves and surface waves. Body waves travel through the Earth’s interior, whereas surfacewaves travel parallel to the Earth’s surface.There are two types of body waves: P-waves (primary or pressure waves) and S-waves (sec-ondary or shear waves). P-waves are the fastest seismic waves and, thus, the first to be detectedby seismometers. They are compressional waves, meaning they move through the Earth by com-pressing and expanding the material in the same direction as the wave is travelling. P-waves cantravel through both solid and liquid layers of the Earth. S-waves are slower than P-waves. Theseare shear waves with particle motions perpendicular to the direction of wave propagation. They canbe divided into SH (horizontal motion) and SV (vertical motion) waves, manifesting 2 degrees ofpolarization compared to 1 for P-waves. Unlike P-waves, S-waves can only travel through solidmaterials, which can aid in identification of fluid presence.There are multiple types of surface waves, such as Love Waves, Rayleigh Waves, and Stoneleywaves. They travel more slowly than body waves. For seismograms of local earthquakes, surfacewaves are typically low amplitude due to the location of sources multiple wavelengths below theEarth’s surface and often indistinguishable from the preceding S-waves. However, they usuallyappear as the largest waves on long-period seismograms of distant earthquakes because such sourcesare, relatively speaking, closer to the surface and surface waves attenuate much more slowly withdistance than body waves [Natural Resources Canada, 2021b].In this thesis, we focus on small earthquakes recorded by local seismic networks. Thus, we donot consider surface waves although, in principle, the techniques described herein could be equallyapplied to them.41.2.3 Moment TensorsA moment tensor (MT) is a numerical representation of the radiation pattern of a seismic pointsource. It is a 3×3 matrix written as:M =M11 M12 M13M21 M22 M23M31 M32 M33 (1.1)where the subscripts 1, 2 and 3 correspond to the East, North and Up directions, respectively. EachMi j element refers to two force vectors of magnitude f , pointing in opposite directions along i,separated by a distance d in the j direction. The magnitude of element Mi j is defined by the productf d and is assumed constant as d goes to zero (limit of a point source) [Shearer, 2009]. Conservationof angular momentum requires that M is symmetric (i.e. Mi j = M ji). Therefore, M has only sixindependent elements and can somewhat arbitrarily be decomposed in 3 different components:M = MDC +MISO+MCLV D (1.2)The double-couple component, MDC, describes excitation from pure shear fractures. Mathemati-cally, it is characterized by a zero trace and determinant. The isotropic component, MISO, is zerowhen no volume change is observed otherwise, it is a scalar times the identity matrix. MCLV Dis known as the compensated linear vector dipole component. It explains extra complexity at thesource. When describing faulting from earthquake sources, the CLVD component is normally min-imized, giving rise to a unique decomposition of M.Moment tensors are key quantities in seismological analysis, providing valuable insights intoearthquake mechanics and the forces at play in a region. By analyzing moment tensors, we can gaindetailed information about rupture processes, including the direction and extent of slip-along faultsand other geological structures. When multiple moment tensors distributed throughout a regionare available, they provide fundamental information for advancing our comprehension of regionalseismotectonics, strongly motivating my research.1.2.4 Focal MechanismsFocal mechanisms (FM) are fault-plane solutions that characterize the deformation-generating seis-mic waves at an earthquake source. They offer a simplified binary representation of the momenttensor, focusing exclusively on the shear fracturing aspect encapsulated by the double-couple com-ponent (MDC). Due to the inherent symmetry of focal mechanisms, they yield information abouttwo possible fault orientations (and slip directions), known as nodal planes. These nodal planes areorthogonal to one other, with one plane representing the actual fault plane along which rupture oc-curred and the other its conjugate that is normal to the slip vector. Determining which of the nodal5planes is the true fault plane requires additional data or contextual information.1.2.5 Beachball DiagramsA beachball diagram is a graphical tool that is widely used in seismology to quickly convey the char-acteristics of seismic sources without reference to source magnitude. It is a stereographic projectionthat illustrates the orientation of the fault plane and the direction of slip during the seismic event,and displays a binary representation of the first-motion P-wave radiation pattern. From a far-fieldview (e.g. seismic station), the ray directions experiencing compression have positive P-wave firstmotions, whereas tension directions show negative first motions. For focal mechanisms, the beach-ball appears as a sphere divided into quadrants of two contrasting colours, which correspond to thecompressional (coloured) and dilatational (white) areas of the seismic wavefield. The boundariesbetween these quadrants represent the nodal planes. To represent the 3D beachball graphically in2D, we generally project the lower hemisphere of the focal sphere onto a horizontal plane. Fig. 1.2is a cartoon to help visualize strike-slip, normal and reverse faulting in the beachball representation.Figure 1.2: Example of beachball diagrams with their block model and half focal sphere(sphere surrounding the earthquake focus). The moment tensor tension and pressureaxes are identified by T and P, respectively. Modified from Rowan [2009].6For a general moment tensor, the focal sphere pattern is no longer a simple beachball; it is morecomplex due to the inclusion of nonzero isotropic and compensated linear vector dipole components(see Fig. 1.3 for an example).Figure 1.3: Example of a beachball diagram for full MT with its decomposition.It is important to note that the tension/compression axes (T/P axes) for a moment tensor aredefined to reside in the middle of the compressional/tensional quadrants, respectively, since thoseaxes refer to the forces experienced at the source, not in the far-field.1.2.6 Moment Tensor Inversion MethodsFor large events (magnitude ≥ 4), moment tensors are routinely and automatically calculated [e.g.,Ekstro¨m et al., 2012; Kao et al., 2012]. The main approach is based on centroid moment tensors(CMT), meaning the moment tensors are estimated based on the average spatial and temporal prop-erties of an often complex seismic rupture, approximated to unidirectional point source. CMTs arerecovered modelling long-period (≥∼ 20s) waveforms using Green’s functions, which describe theimpulse response of the medium between the seismic event and the seismometer. CMT solutionsfor small events are hard to estimate because low-frequency waves are typically buried in noise.There are several techniques for recovering the focal mechanisms for small seismic events.The classic and simplest one involves identifying the first motion polarities of P-waves across anetwork of stations encompassing diverse azimuth and take-off (angle from the down vertical formwhich a seismic ray leave the focal sphere) angles. Thereafter, two orthogonal planes demarcatingpositive and negative polarities are determined, often employing grid search methods. Althoughstraightforward, choosing the P-wave polarity can be ambiguous and time-consuming. Further,the station distribution strongly defines the uncertainties of the retrieved focal mechanisms. P-polarities are sometimes augmented with S-to-P amplitude ratios to provide further constraints onfocal mechanisms. However, Hardebeck and Shearer [2003] note minimal improvements in focalmechanism recovery, likely due to large measurement uncertainties from the poorly constrained Pand S-velocity and attenuation models.Recent contributions by Cheng et al. [2023] and Zhang [2023] introduced inversion methodsbased on S/P-wave amplitude ratios. These techniques determine focal mechanisms, a subset of all7possible seismic patterns and do not characterize magnitude. Choi et al. [2023] propose a supervisedmachine learning approach to retrieval of MTs. Their method seems promising; however, a dense,local and regular seismic network is required (e.g. 21 x 41 stations), and the technique has onlybeen tested on synthetic microseismic events of the same magnitude.Another recent approach for improving moment tensor estimation of small earthquakes is basedon relative amplitudes of body waves [Dahm, 1996; Plourde and Bostock, 2019]. It demonstratessignificant proficiency in recovering several MTs with minimal human intervention. This techniquerequires events to be close enough to one another that Green’s functions are similar and that vari-ations therein can be neglected. Under this assumption, the relative amplitudes of P and S-waveseismograms between events recorded at a specific station depend upon the relative scalar momentsof the events and their radiation patterns. Plourde and Bostock [2019] advocate the use of principalcomponent analysis (PCA) to calculate relative amplitudes accurately. However, the relative MTalgorithm requires at least one known MT to fix the polarity of all MTs since relative informationdoes not provide absolute values. This method has demonstrated good recovery accuracy for smallevents [Plourde and Bostock, 2019]. It is used in Chapter 2 and is improved in Chapter 3 whichprovide a detailed overview of the method.1.2.7 Stress Inversion MethodStress inversion is a powerful technique in seismology that aims to determine the state of stressin a localized region by analyzing the focal mechanisms of earthquakes. It can be viewed as aform of averaging tool for focal mechanisms. The general purpose of stress inversion is to estimatethe orientation and relative magnitudes of the principal stress axes that drive faulting and seismicactivity in a given region. We define the three principal axes as σ1, σ2 and σ3 where σ1 ≥ σ2 ≥ σ3.σ1 and σ3 refer to the maximum and least compression stress orientation, respectively.In this thesis, we have employed the stress inversion technique developed by Vavrycˇuk [2014],which does not require prior knowledge of fault plane orientation. This method accurately recoversthe principal stress axes and the stress ratio(R = σ1−σ2σ1−σ3), which characterizes the relative magni-tudes of the principal stress axes. It builds upon the Michael [1984] approach that involves solutionof a least-squares system to minimize the difference between the observed slip direction and thepredicted tangential traction while ensuring that the magnitude of this traction remains significantand consistent across the various fault planes. The tangential traction refers to the force acting par-allel to the fault plane that generates sliding motion. However, this technique requires knowledgeof the fault plane orientation. Vavrycˇuk [2014] uses an iterative process to identify the fault planeorientation as it inverts for stress using the Michael [1984] method.81.3 Thesis OutlineThis thesis is structured around three primary research projects (Chapters 2, 3, and 4) that I con-ducted during my PhD. Each research Chapter consists of manuscripts that have been or will besubmitted to peer-reviewed scientific journals (see Preface for further details).In Chapter 2, I investigate two significant intraslab aftershock sequences for large intraplateearthquakes in southern Alaska, USA: the M7.1 Iniskin 2016 earthquake and the M7.1 Anchorage2018 earthquake. By employing the relative moment tensor algorithm from Plourde and Bostock[2019] on hundreds of small earthquakes, I analyze the spatiotemporal distribution of aftershocksto gain insights into the rupture processes of the mainshocks. This Chapter is especially relevant toresearch question 1.Recognizing potential improvements in the relative moment tensor inversion method, I ded-icated Chapter 3 to developing an enhanced algorithm that incorporates additional constraints,specifically aimed at improving the analysis of microseismic events. The algorithm’s performanceis tested on synthetic examples. Subsequently, the inversion method is applied on a dense seismiccluster near San Juan Island in Washington, USA, which yields insight into the nature of this con-centration of seismic activity. This Chapter addresses research question 2 and provides insight intoquestion 1.In chapter 4, I apply the newly developed relative amplitude algorithm to seismic activity inthe Haida Gwaii region of British Columbia, Canada. This region presents significant challenges tostudy of microseismic events due to its sparse seismic network coverage. In this chapter, I explorethe algorithm’s limitations and investigate several earthquake groups within the incipient subductionzone of Haida Gwaii. This Chapter focusses on research question 1 but also touches upon question2.Finally, Chapter 5 provides a brief conclusion that revisits my research questions and summa-rizes the key findings of this thesis. Additionally, the appendix includes supplementary materialssuch as figures, discussions and a data availability statement for each research chapter.9Chapter 2Aftershock Distributions, MomentTensors and Stress Evolution of the 2016Iniskin & 2018 Anchorage Mw 7.1Alaskan Intraslab Earthquakes2.1 Chapter SummaryThis Chapter presents a detailed study of two Mw 7.1 intraslab earthquakes that occurred in southernAlaska: the Iniskin earthquake of January 24th, 2016, and the Anchorage earthquake of November30th, 2018. Moment tensors were relocated and recovered for hundreds of aftershocks followingboth events, and stress histories were inverted. The aftershock distribution of the Iniskin earthquakesuggests that the rupture propagated updip along a fault dipping steeply into the Pacific Plate andterminated at a stratigraphic horizon, inferred to be either the interface or Moho of the subductingslab. In addition, four earthquakes ruptured the main fault in the preceding two years and had simi-lar moment tensors to the mainshock. This evidence suggests that the mainshock likely reactivateda pre-existing, outer-rise fault. The Anchorage earthquake sequence is complex due to its locationnear the boundary of the subducting Yakutat and Pacific plates, as evidenced by the aftershock dis-tribution. Aftershock hypocenters form two main clusters that appear to correspond to orthogonal,conjugate faults, consistent with the two nodal planes of the dominant focal mechanisms. Bothgeographic groups display many focal mechanisms similar to the mainshock, which could indicatesimultaneous rupture on conjugate planes. The time dependence in stress ratio for the Iniskin se-quence can be interpreted in terms of pore-pressure evolution within the mainshock fault zone. Inparticular, the observations are consistent with a dehydration-assisted transfer mechanism wherefluids are produced during rupture through antigorite dehydration and raised to high pore pressures10through matrix collapse and/or thermal pressurization. The Anchorage sequence exhibits a morecomplex stress ratio evolution that may be associated with stress adjustments within a distributedfault network, or reflect a strongly heterogeneous stress field.2.2 IntroductionIntraplate earthquakes at ocean-continent convergent margins are a significant source of seismichazard to nearby coastal communities. Notable examples of damaging intraslab earthquakes includethe deadly Mw 7.1 2017 Mexico earthquake, the Mw 6.1 2001 Geiyo earthquake in Japan, the Mw6.8, 2001 Nisqually earthquake near Seattle, the Mw 7.7 2001 El Salvador earthquake in Guatemalaand the Mw 7.1, 1997 Punitaqui earthquake in central Chile. These events can be more destructivethan larger megathrust events due to their greater proximity to population centres, shorter recurrenceintervals and higher ratio of seismic energy to seismic moment [Choy et al., 2001; Izutani andKanamori, 2001; Kirby et al., 2002; Preston et al., 2003]The mechanisms of intraslab earthquakes deeper than 30 km are poorly understood and are thetopic of long-standing study. Proposed candidates for events between 30 km and 300 km depthinclude [Hasegawa and Nakajima, 2017]: (1) thermal shear instability of plastic deformation, alsocalled self-localizing thermal shear runaway, and (2) dehydration embrittlement. The first mecha-nism requires that plastic deformation proceeds at a faster rate than thermal dissipation such that theslab melts and slipping is encouraged. The second mechanism assumes that hydrous minerals arepresent in the slab as it subducts and that spontaneous dehydration elevates pore pressures, reducingeffective stress and thus enabling slip instability [Kirby, 1995]. The leading candidate dehydrationreactions are the transformations of metabasalt to eclogite and serpentinite to peridotite in crustand mantle lithosphere, respectively [Hasegawa and Nakajima, 2017; Peacock, 2001; Peacock andWang, 1999].Numerous studies [Abers et al., 2013; Hacker et al., 2003; Hasegawa and Nakajima, 2017;Shiina et al., 2013; Wang, 2002; Yamasaki and Seno, 2003] have appealed to the dehydration em-brittlement mechanism. However, recent laboratory experiments do not support the dehydrationembrittlement hypothesis. For example, despite the creation of conditions favourable to dehydra-tion instability within antigorite, Chernak and Hirth [2011] found no evidence for brittle deforma-tion, but instead observed slow-slip and stress relaxation events. Gasc et al. [2011] did not detectany acoustic emissions during their serpentinite dehydration experiments. More recently, Proc-tor and Hirth [2015] demonstrated that dehydrating serpentine deforms through semibrittle flow,and thus concluded that dehydration embrittlement is unlikely to be the cause of intermediate/deepevents. Moreover, antigorite dehydration at pressure-temperature conditions representative of thosefor intermediate-depth earthquakes is accompanied by volume decrease [Ferrand et al., 2017; Junget al., 2004] and so contradicts the basis of the dehydration embrittlement mechanism. To addressthis problem, Ferrand et al. [2017] propose an alternative model to classical dehydration embrit-11tlement, in which dehydration driven stress transfer, rather than fluid pressure increase, promotesfailure.Alternatively, numerical modelling by Kelemen and Hirth [2007] suggests that thermal shear in-stability can occur within a temperature range (600-800◦C) expected for intermediate-depth earth-quakes, and thus could afford a potential explanation for their genesis. Furthermore, John et al.[2009] reproduced field observations using fully coupled viscoelastic and thermal models by simu-lating self-localizing thermal runaway. However, Prieto et al. [2013] argued that thermal runawaymust occur after slip has initiated, and that a trigger mechanism remains to be identified. As of yet,there is thus no clear consensus as to which of the two prevailing mechanisms, dehydration embrit-tlement or thermal shear runaway, is to be preferred. Indeed, a combination of both mechanismscould be the key to intraslab earthquake generation. Dehydration reactions could produce fine-grained zones along e.g., preexisting faults [Peacock, 2001], upon which thermal shear instabilitymay occur [Kelemen and Hirth, 2007; Prieto et al., 2013]. Further research is clearly required.In this study, the attention is focused on two Mw 7.1 intraslab earthquakes in southern Alaskaand their aftershock sequences: the January 24, 2016 Iniskin and November 30, 2018 Anchorageearthquakes. Both events were followed by over a hundred aftershocks of Mw ≥ 3. A large num-ber of stations of the USArray Transportable Array, along with permanent stations from the USGSVolcano Observatory and the Alaska Earthquake Center, recorded both events thereby providing anexceptional opportunity to study large intraslab earthquakes in some detail. The analysis of the af-tershock distributions, moment tensor solutions using a relative moment tensor inversion approach,and the temporal evolution of local stress throughout the aftershock period is presented. The analy-sis yields insights into the rupture process and mechanism of the two earthquakes.2.3 DatasetWaveforms for all events were requested from the Incorporated Research Institutions for Seismol-ogy (IRIS) for over 200 mainly 3-component stations, as shown in Fig. 2.1. These stations liewithin a 2.5◦ epicentral radius from studied events and belong to several networks: Alaska Re-gional Network, National Tsunami Warning Center Alaska Seismic Network, Alaska Volcano Ob-servatory, Global Seismograph Network, United States National Strong-Motion Network, USArrayTransportable Array, US Geological Survey Networks, Alaska Amphibious Community SeismicExperiment and Southern Alaska Lithosphere and Mantle Observation Network.2.3.1 January 24, 2016, Iniskin EventThe Mw 7.1 Iniskin earthquake struck at 127 km depth below Iniskin Bay in southern Alaska. Thestrike-slip rupture propagated updip along a plane dipping at 66◦ with strike of 58◦, penetrating thePacific slab, as represented in the finite fault model from U.S. Geological Survey [2017e]. It was12Figure 2.1: Station distribution used in this study. Orange, red and black triangles representstations used for the Iniskin sequence, Anchorage sequence or both sequences, respec-tively. Orange and red stars represent the mainshocks of the Iniskin and Anchorage se-quence, respectively. Dots are earthquakes below 40km depth of magnitude 3 and greaterthat occurred between January 1, 2000 and December 31, 2021. The E115 Pacific platecontours [Jadamec and Billen, 2010] are displayed in green solid lines and the trench ismarked by the black, toothed line. The blue polygon in the inset identifies the regionshown.followed by more than 140 aftershocks of Mw ≥ 3, ranging in depth from 90 km to 155 km. Toanalyze this sequence, 142 aftershocks were combined with 44 previous earthquakes of magnitude≥ 3, for a total of 186 events, occurring between January 1, 2014, and November 30, 2018. Themainshock occurred within the Pacific plate which subducts at a rate of ∼60mm/yr relative to theNorth American plate. The Pacific slab dips at ∼44◦ and strikes at ∼204◦, as calculated from theE115 slab model of Jadamec and Billen [2010]. At the trench near the site of the Iniskin event, theage of the Pacific plate is estimated to be ∼45 Ma [Mu¨ller et al., 1997].13Figure 2.2: A) Iniskin sequence earthquake relocation. Rupture plane of the mainshock ishighlighted in red on its lower hemisphere focal mechanism. B) Anchorage sequenceearthquake relocation. The two ellipses indicate the two groups of aftershocks that alignapproximately with the two respectively coloured nodal planes on the lower hemispherefocal mechanism. For A and B, the E115 Pacific plate contours [Jadamec and Billen,2010] are displayed in green solid lines. C) Cross-section showing all events projectedon plane I1-I2 from A. D) Cross-section showing all events projected on plane A1-A2from B.142.3.2 November 30, 2018, Anchorage EventThe Anchorage earthquake occurred ∼15 km north of Anchorage at a depth of ∼52 km, and wasfollowed by more than 400 aftershocks of Mw≥3. In the analysis, 447 events are used, including418 aftershocks, of magnitude ≥3 between January 1, 2014, and December 31, 2020. The main-shock was characterized by normal faulting [U.S. Geological Survey, 2017d] and occurred eitherwithin the Yakutat plate or the Pacific plate. The teleseismic migration and receiver function anal-ysis of Kim et al. [2014] revealed that the location of the earthquake is in a transitional zone wherethe Yakutat and Pacific plates cannot be clearly distinguished. The mainshock centroid and its af-tershocks appear to lie near the western edge of the Yakutat plate. Furthermore, a major bend inthe subducting plate is evident in the vicinity of the sequence Kim et al. [2014]. The Yakutat plate,also called the Yakutat terrane or block, is subducting at ∼50mm/yr N23◦W Elliott et al. [2010],whereas the ∼35 Ma Pacific plate is thought to be subducting at a similar rate and direction at thislocation.Several studies have examined the mechanics of the mainshock [Guo et al., 2020; Liu et al.,2019; Ruppert et al., 2019; West et al., 2019]. Results from the first three studies are ambiguousas concerns which nodal plane of the (normal) mainshock mechanism hosted the rupture. Theambiguity stems from the fact that aftershocks concentrate in two separate groups (see Fig. 2.2B). The northern group favours the steep, west-dipping nodal plane of the mainshock, whereasthe southern group favours the east-dipping plane. Fig. 2.4 A displays the distribution in depthof the groups along the respective best-fit planes. Guo et al. [2020] interpreted the mainshock asa simultaneous rupture of the two conjugate planes. Using strong-motion waveform inversions,they proposed a kinematic rupture model where the west-dipping plane and deeper portions of theeast-dipping plane experience the most slip. They argue that the shallow part of the east-dippingplane, with a more diffuse distribution of aftershocks, slipped only slightly. The complexities of theearthquake sequence might be associated with its location near the Pacific-Yakutat plate boundaryand/or a plate bend in the slab [Kim et al., 2014].2.4 MethodologyThree main processing steps were performed to analyse the earthquake sequences: 1) hypocentralrelocation, 2) moment tensor inversion, 3) temporal stress inversion. In this section, each of thesesteps is described.2.4.1 Earthquake RelocationsBoth earthquake sequences were initially relocated using the Hypoinverse program [Klein, 2014]with manually picked P- and S- waves arrival times to secure consistent, absolute locations. Thenthe double-difference algorithm, HypoDD [Waldhauser and Ellsworth, 2000], was applied for more15accurate relative relocation. A combination of manual picks, cross-correlation delay times andthe Hypoinverse initial locations were used as input to HypoDD. The automatic delay times wereobtained using the principal-component-based, multichannel alignment method of Bostock et al.[2021]. Note that these relocations represent earthquake centroids rather than hypocentres due tothe use of correlation-based delay times, although the two location measures are close for smallevents. As part of the relocation process, 1D-velocity models were employed, extracted from theNayak et al. [2020] 3D-velocity model corresponding to the central location of each sequence,modified to remove low velocity zones.2.4.2 Moment Tensor InversionsMethodMoment tensors for hundreds of events with Mw ≥ 3 were recovered using the relative momenttensor inversion developed by Plourde and Bostock [2019] (see also Dahm [1996]). This methodemploys a principal component description of P- and S-waveforms recorded at 3-component stationsto compute relative amplitudes between aligned waveforms of pairs (P) and triplets (S) of eventsrecorded at a common station. Events must be close enough to one another that their Green’sfunction can be assumed to be the same. In this case, relative amplitudes will depend upon therelative scalar moments of the events and their radiation patterns. Relative waveform amplitudeswere calculated for the frequency band maximizing the signal-to-noise ratio for each pair or tripleof events from 25 frequency bands ranging between [0.1, 1]Hz to [4, 8]Hz. Further, the windowlength was adjusted with respect to the frequency band so each window encompassed a comparablenumber of cycles. The algorithm also requires at least one known focal mechanism, reffered to as thereference event(s), to fix the focal mechanism polarity. All reference event focal mechanisms weretaken from the USGS website to ensure high quality [Helffrich, 1997]. Moment tensor recoveryinvolves solution of a large system of equations that comprises primarily relative moment tensorconstraints coupled to a small number of reference event equations. These latter equations areassigned large weights, 1000 times higher than the relative constraints, to ensure that outliers do notsignificantly perturb the reference event MTs. The Anchorage system (5,482,412 equations × 2175MT elements) has roughly five times more equations than the Iniskin system (1,256,996 equations× 835 MT elements) due to the larger number of associated aftershocks. For both inversions,double-couple (DC) solutions were recovered by enforcing a zero determinant for each MT throughnon-linear optimization.In addition, 38 and 14 focal mechanisms were incorporated as additional constraints within theIniskin and Anchorage systems, respectively, using HASH [Hardebeck and Shearer, 2002]. P-wavefirst motion polarities for high-quality recordings at given stations were assigned after waveformalignment [Bostock et al., 2021]. The HASH solutions were determined through a grid search with16a resolution of 5◦, and assigned a weight of 100 with respect to the relative constraints, and werethus allowed to vary from their initial values. This weight was determined based on an L-curveanalysis; see Fig. A.2 in the Appendix A.Furthermore, the inversion included only those events that met the following selection criteria:a minimum of 10 P and 10 S comparisons, a minimum sum of 30 P plus S comparisons, anda minimum of 8 stations contributing to the comparisons. Moreover, waveform comparisons weremade only between events that are at most 12 km apart to ensure Green function similarity, and wereaccepted only if the absolute cross-correlation values for P-waves are greater or equal to 0.8. S-waveconstraints involve 3 waveforms due to the 2 degrees of freedom for S-polarizations. Accordingly,one waveform can be reconstructed as a linear combination of the other two in the absence ofdegeneracy. S-wave constraints are selected based on the requirement that the misfit and cross-correlation between a waveform and its reconstruction using its two principal components be lowerthan 0.5 and above 0.8/0.9 (Iniskin/Anchorage), respectively. The misfit is the norm of the differencebetween the waveform normalized to unit magnitude and its reconstruction.Quality TestA bootstrapping approach was employed to measure the stability of each of the recovered MTs.A total of 100/120 (Iniskin/Anchorage) bootstrapping iterations over the comparison equations wasperformed and solutions from the 70 iterations that fully converged to DC MTs were kept. Then, thebootstrapped MTs for each event were compared to the preferred (full data set) solution using theKagan angle [Kagan, 1991] as a metric. The Kagan angle represents the minimum angle required torotate one DC focal mechanism into another orientation. Its value ranges from 0◦ indicating absoluteagreement to 120◦ for complete disagreement [D’Amico et al., 2011]. The standard deviation of theKagan angles for each event is considered as a measure of uncertainty.The quality test results are displayed in the Appendix A (Figs A.3 & A.4). For the Iniskinsequence, only 17 events have a Kagan angle standard deviation above 5◦. Most of these events oc-curred before August 2015, when the station density was lower (the transportable array deploymentstarted in 2014 and terminated in 2017). For the majority of the events, the recovered MTs are stable.For the Anchorage sequence, 63 events of a total 435 have a Kagan angle standard deviation above5◦, with 40 above 10◦. For both sequences, the less stable events are those with fewer comparisonsand/or poorer coverage over the focal sphere. In addition, the Kagan angle differences between 2/7(Iniskin/Anchorage) events and corresponding Global Centroid Moment Tensor (GCMT) solutions[Dziewonski et al., 1981; Ekstro¨m et al., 2012] were computed. The differences are of the sameorder as those between the correspondng USGS and GCMT solutions.17Clustering MethodThe MTs were clustered using the unweighted pair group with arithmetic mean (UPGMA) method.This algorithm allows us to cluster by similarity without specifying the number of clusters. A max-imum of 25◦ Kagan angle difference between MTs and a minimum of 10 elements per cluster wererequired. A 25◦ Kagan angle provided a good compromise between similarity and a numbers of un-clustered elements. For representation in Figures, the mean DC MT of each cluster was computed.2.4.3 Temporal Stress InversionsMethodTo analyze possible stress changes due to and following the mainshocks, stress inversions wereperformed as functions of time with the recovered MTs for both earthquake sequences using theVavrycˇuk [2014] method. The mainshock focal mechanisms were not included in the inversions.Errors on focal mechanisms were assigned as the standard deviation in Kagan angle obtained fromthe bootstrapping quality tests, with a fixed minimum error of 5◦. Stress inversions employed amoving window in time of 30 events (i.e. adjacent windows share 29 events). Martı´nez-Garzo´net al. [2016] have shown that for focal mechanisms with noise levels up to 20◦/40◦, a minimumof 20/40 events is necessary to retrieve the stress field orientation with accuracy of 10◦. Thus, a30-event window ensures reasonable focal mechanism diversity as well as stable inverse solutionsand reasonable resolution.Quality TestsResults from the time-dependent stress inversions are sensitive to the technique used to select thefault plane from the two nodal plane solutions. The iterative instability method as developed byVavrycˇuk [2014] was employed. Specifically, for each iteration, the stress tensor is computed andthe nodal plane with lesser stability, according to the Mohr-Coulomb failure criterion, is adoptedas the fault plane for the following iteration. As previously demonstrated [Lund and Slunga, 1999;Martı´nez-Garzo´n et al., 2016], it was observed that the alternative approach of choosing the planewith lowest misfit angle between the focal mechanism slip vector and shear traction from the stresstensor yielded less stable stress results, as did the efforts to combine the two methods.To assess the quality of the stress inversions, the four metrics recommended by Hardebeck andHauksson [2001] were employed, namely resolution, uncertainty, diversity, and misfit. Temporalresolution is determined by the root mean square (RMS) time computed relative to the mediantime of those events contributing to the inversion. The uncertainty parameter is defined as U =(R∆σ1)+ (1−R)∆σ3 where R is the stress ratio that will be introduced in section 2.5.3 and ∆σiis the uncertainty in the ith principal stress direction corresponding to the 1-σ error over 500 noise18realisations. The noise level is taken as the± error assigned to each MT. This uncertainty calculationaccounts for the fact that during uniaxial compression (tension) the uncertainty in σ3 (σ1) is large,due to the corresponding axis being indistinguishable from that of σ2. Focal mechanism diversity iscomputed as the RMS Kagan-angular difference from the average mechanism. The average misfitis the average of the angles between the slip direction and the predicted tangential traction from thestress tensor.To ensure reliable results, Hardebeck and Hauksson [2001] demonstrate that diversity must be atleast∼35◦– 45◦ for focal mechanisms with 0◦– 20◦ errors. Michael [1991] suggests that an averagemisfit lower than ∼35◦ will accurately recover the uniform component of the stress tensor underthe assumption of a homogeneous stress field for focal mechanisms with errors of 10◦– 20◦. Bythese measures, the diversity of the recovered focal mechanisms falls in the low range of acceptablevalues, but the misfits and uncertainties indicate stable inversions. It is assumed that the results arereliable and accurately reflect stress changes.2.5 Results2.5.1 Hypocenteral DistributionsIniskinFig. 2.3 displays a 3D representation of the earthquake sequence locations. Note that most ofthe events before the mainshock (light blue) are distributed parallel to the plate interface. Moreover,approximately half of the aftershocks (69 events) and one foreshock define a∼50 km long lineation.The one foreshock is a Mw 4 event that ruptured one month before the mainshock at the north-eastend of the lineation, ∼30 km away from the mainshock hypocentre. The orientation of the linearstructure, as determined by principal component regression, lies within 4◦ of the intersection ofthe mainshock rupture plane, as reported by the U.S. Geological Survey [U.S. Geological Survey,2017e], with the E115 slab model. At these depths, the E115 model is defined on the basis ofintraplate seismicity and so it is conceivable that the true slab surface is dipping up to 10◦ lesssteeply, and that the top of the slab is therefore up to ∼10 km shallower. In the Appendix A.3, adiscussion about the choice of E115 model [Jadamec and Billen, 2010] vs Slab 2 model [Hayes,2018] is provided. Consequently, the lineation plausibly represents an intersection between therupture plane of the mainshock and stratification within the slab, for example, the slab surface or theoceanic Moho. An additional 46 events were observed to lie on or near the rupture plane, includingthree events that occurred before the mainshock (one in 2014 and two in 2015). The angle betweenthe rupture plane and the slab model is∼76◦. The plane defined by the aftershocks penetrates∼ 9.5km below the lineation. The mainshock centroid is located ∼7.6 km below the lineation.19−30−20−10010Along-dip distance (km) 0 20 40BAC −20020Distance (km)−20 0 20Horizontal distance (km)Horizontal distance (km)Figure 2.3: Relocated events near the Iniskin earthquake. A) 3D representation highlightingthe earthquake lineation (viewed from 37.5◦ azimuth and 30◦ elevation). The horizontalorigin is centered on the mainshock centroid (yellow star). The mainshock rupture planefor the USGS solution is shown in pink and the local orientation of the Jadamec andBillen [2010] slab model is in green. Event symbols in light and dark blue earthquakesdistinguish foreshocks and aftershocks, respectively. Events with purple rims were usedto fit the lineation shown by the purple line. A video of this 3D representation is avail-able in the Supplementary Material of Drolet et al. [2022]. B) View showing the eventsprojected on the rupture plane. Filled events lie within the rupture plane. C) View orthog-onal to the intersection of the rupture plane and slab model (viewed from -48.6◦ azimuthand -20.3◦ elevation).AnchorageThe Anchorage sequence has been previously relocated by Ruppert et al. [2019]. Both their studyand this one have very similar seismicity patterns and identify a northern and southern cluster.However, as shown in Fig. A.1 of the Appendix A, this study relocation average ∼3 km shallower,∼3.5 km eastward, and the northern cluster is more localized. These differences are likely due tothe velocity model and minor differences in measured delays.Fig. 2.2 B and D show the map and cross-section views of the sequence, respectively whereasFig. 2.4 displays the 3D distributions. Recall that two main event groups formed, the northern group(142 events) favouring the steep west-dipping nodal plane of the mainshock and the southern group20BAFigure 2.4: 3D representation of relocated events near the Anchorage earthquake. The hori-zontal origin is centred on the mainshock centroid (purple star). A) View showing thetwo nodal planes from a perspective directed along their axis of intersection, and associ-ated seismic events. Event symbols with pink and yellow rims are empirically associatedwith the southern and northern groups, respectively. Earthquake symbols in light anddark blue correspond to foreshocks and aftershocks, respectively. A video of this 3Drepresentation is available in the Supplementary Material of Drolet et al. [2022]. B) En-larged view of the northern group including the best-fit plane shown as a red line. Yellowsymbols with red rims represent events used for plane fitting. Yellow and pink eventsform part of the northern and southern groups, respectively, whereas blue events are notclearly associated with either group.(117 events) favouring the east-dipping plane. The alignment of the northern aftershock group withthe corresponding nodal plane is reasonable. As shown in Fig. 2.4 B, a total of 108 events (yellowcircles with red rims) in the northern group form a reasonably well-defined plane with dip 70◦±2◦ and strike 220◦± 2◦. The quoted uncertainties are 1σ as determined by bootstrapping over theevents that define the plane. Three events on the plane occurred between 2014 and 2016 before themainshock. Only one event from this group has a USGS moment tensor solution [U.S. GeologicalSurvey, 2017a]: a Mw 5 earthquake at the deeper end that occurred 14 hours after the mainshock.One of its nodal planes is similar to the fitted plane description with a strike of 204◦, dip of 73◦ andrake of 291◦. In contrast, the southern aftershock group is less concentrated without any particularlyremarkable geometrical attributes.2.5.2 Moment tensorsIniskinA total of 167 MT solutions were recovered, including 30 events preceding the mainshock, usingfour reference events and 38 HASH MT solutions. The reference event focal mechanisms are shownin Fig. A.5 in Appendix A. The solutions are mainly strike-slip mechanisms, similar to the main-21Strike0901802702040Dip0 3060902060Rake0901802702040090180270480 30609051009018027048090180270480 306090102009018027048153°30'W 153°15'W 153°00'W 152°45'W59°30'N59°45'N60°00'N#1#2RestNo MTBA Restx 93 x 18 x 56#1 #2Figure 2.5: Moment tensor clusters for the Iniskin sequence. A) Top row: Mean momenttensor solutions for each cluster (lower hemisphere representation). The green grouprepresents not clustered MTs, a random MT from this group is displayed. The numberof elements per group is written under each MT and their cluster number in the boxes.A) Lower rows: Distributions of strike, dip and rake angles of each nodal plane for MTsin each cluster. B) Locations of events in each MT cluster. Events with no recoveredmoment tensors are displayed in gray.shock, with varying degrees of thrust. Fig. 2.5 displays the mean moment tensor for each clusterand the corresponding nodal plane angles. Cluster #2, almost pure strike-slip, includes the Mw 4foreshock and events located primarily at the northern end of the lineation. The unclustered events(light green) are scattered throughout the region. Cluster #1 is characterized by mixed thrust/strike-slip mechanisms, represents the majority of the events, and includes the mainshock. A total of 70events in cluster #1 lie in a plane with a width less than 3.5 km parallel to the USGS rupture planewhich is interpreted to be the same. One of the cluster #1 nodal planes has a median strike of ∼54◦and a corresponding median rake of ∼26◦. The seismicity lineation (Fig. 2.3) is oriented at ∼25◦from the horizontal in the rupture plane frame (or ∼22◦ if the intersection between slab model andrupture plane is employed). The significance of this observation is discussed in section 2.6.1. Notethat three foreshocks are also included within cluster #1 and in proximity to the rupture plane. Fig.2.6 displays examples of waveform comparisons used for the MT inversion and the correspondingstation locations on the upper focal hemisphere.222015-05-25 22h40 1 2 3 4 5 6 72017-02-07 17h52 1 2 3 4 5 6 70 1 2 0 1 2 0 1 2Time [s]7 - AV.NCT6 - AK.SSN5 - AV.SPCP4 - AV.AUJA3 - AT.SVW22 - AK.CNP1 - AK.BRLKP-waves 0.94583P-waves 0.8634P-waves 0.82638P-waves 0.81769P-waves 0.94747S-waves 0.87206S-waves 0.907432017-02-07 17h52 1 2 3 4 5 6 72018-04-02 08h45 1 2 3 4 5 6 70 1 2 0 1 2 0 1 2Time [s]7 - TA.P19K6 - TA.P18K5 - TA.N19K4 - TA.M20K3 - TA.L19K2 - AV.SPBG1 - AV.ANCKP-waves 0.94346P-waves -0.95527P-waves 0.94946P-waves -0.82299P-waves -0.88627P-waves 0.87533P-waves 0.92454North VerticalEastNorth VerticalEastFigure 2.6: Waveform comparisons with corresponding stations projected onto focal sphere.Waveforms are filtered to maximize their signal-to-noise ratio and scaled so that relativeamplitudes of each pair is maintained. The cross-correlation coefficients between wave-forms are listed on the right side of each plot. Focal mechanisms are plotted in the upperhemisphere representation, so the P-wave polarities can be matched the radiation pattern.Top: Events separated by 2.6 km with different MTs. Bottom: Events separated by 3.4km with similar MTs. Note the February 7, 2017 (17h52) was used in both comparison.AnchorageA total of 435 moment tensors were recovered for the Anchorage sequence, including 24 eventspreceding the mainshock. U.S. Geological Survey [2017d] solutions for 49 MTs were incorporatedwithin the inversion. Four USGS MTs were selected as reference events to maximize the number ofcomparisons with other events and fault mechanism diversity. Those moment tensors are shown inFig. A.6 of Appendix A. Forty-five extra USGS MTs were added to the equation set and weightedas the 14 HASH focal mechanism constraints to avoid scaling issues while solving for DC solutions.Fig. 2.7 displays the mean MTs for each cluster along with their respective event locations. Nine23x 209 x 28 x 21x 19 x 18 x 14x 13 x 10 x 10x 93Rest#1 #2 #3#4 #5#7#6#8 #9150°15'W 150°00'W 149°45'W 149°30'W61°15'N61°30'N#1#2#3#4#5#6#7#8#9RestNo MTA BFigure 2.7: Moment tensor clusters for the Anchorage sequence. A) Mean moment tensorsolutions for each cluster (lower hemisphere representation). The green group representsnot clustered MTs, a random MT from this group is displayed. The number of elementsper group is written under each MT and their cluster number in the boxes. B) Locationsof events in each MT cluster. Events with no recovered moment tensors are displayed ingray.clusters were obtained that represent primarily normal faulting with minor strike-slip components,and that display greater diversity than the Iniskin solutions. MTs from 93 events were left unclus-tered (light green) and are geographically scattered. The largest cluster, #1, constitutes almost halfof all recovered MT events including the mainshock, the majority of the southern group (83 events),and part of the northern group (42 events). Cluster #3 (21 MTs) events tend to be located to thenorth-west of both groups and cluster #4 (19 MTs) is dispersed around both groups. Cluster #5 (18MTs) concentrates at the north end of the northern group, cluster #9 (10 MTs) falls along its south-ern end, and cluster #7 (13 MTs) occurs within its central portions. Cluster #6 (14 MTs) is locatedbetween the two groups. Finally, clusters #2 (28 MTs) and #8 (10 MTs) are scattered throughout theregion. Interestingly, focal mechanisms similar to the mainshock (cluster #1) are distributed acrossboth geographic groups. Fig. 2.8 displays examples of waveform comparisons used for the MTinversion and the corresponding station locations on the upper focal hemisphere. These examplesdemonstrate that some events with similar MTs from different geographic groups can display highlycorrelated waveforms, and that some events with different MTs belonging to the same group, can242018-11-30 18h30 1 2 3 4 5 6 72018-11-30 23h45 1 2 3 4 5 6 72018-12-24 20h46 1 2 3 4 5 6 72018-12-31 03h28 1 2 3 4 5 6 70 1 2 0 1 2 0 1 2Time [s]7 - GM.AD036 - AK.K2205 - TA.L20K4 - GM.AD093 - AV.SPBG2 - AV.RDJH1 - AK.SLKP-waves 0.97906P-waves 0.83156P-waves 0.94202P-waves 0.89818P-waves 0.84648S-waves 0.85975S-waves 0.851110 1 2 0 1 2 0 1 2Time [s]7 - AK.SSN6 - AK.K2155 - TA.M22K4 - TA.M20K3 - AV.SPBL2 - AK.K2231 - AK.K220P-waves -0.98472P-waves -0.89522P-waves 0.85354P-waves 0.92191P-waves 0.87177S-waves 0.83407S-waves 0.82551North VerticalEastNorth VerticalEastFigure 2.8: Waveform comparisons with corresponding stations projected onto focal sphere.Waveforms are filtered to maximize their signal to noise ratio and scaled so that relativeamplitudes of each pair is maintained. The cross-correlation coefficients between wave-forms are listed on the right side of each plot. Focal mechanisms are plotted in the upperhemisphere representation, so the P-wave polarities can be matched the radiation pattern.Top: Events from the northern group separated by 1.7 km with different MTs. Bottom:Events from north and south groups, separated by 10.0 km with similar MTs.display negative waveform correlations. Table 2.1 summarizes nodal plane orientations of the maingroups and cluster #1 along with those of the mainshock.It was previously demonstrated that most events within the northern group fall approximatelyalong a plane. Not surprisingly, the MTs recovered for these events share a common, similar nodalplane. Fig. 2.9 displays the nodal-plane orientation distributions for the events assigned to the plane.One foreshock, which occurred in February 2015, is part of cluster #1 and is located near the centre-bottom of the event plane. The Mw 5.0 (December 1, 2018) event which has a USGS MT solution25Table 2.1: Nodal plane angles associated with different groups or events within the Anchoragesequence. Note that each event group has a steep west-dipping (first plane) and a shalloweast dipping (second plane) plane, and that the west dipping planes are all within 10◦ ofone another.Event Group Strike Dip RakeMainshock 189◦ 62◦ 272◦6◦ 28◦ 267◦Northern Group(within event plane)197◦ 71◦ 284◦324◦ 29◦ 231◦Southern Group 199◦ 62◦ 272◦330◦ 29◦ 233◦Cluster #1 195◦ 69◦ 288◦332◦ 28◦ 231◦Strike0901802701020Dip0 3060902040Rake0901802701030Figure 2.9: Nodal plane angles of the 108 MTs that form the event plane defined by the north-ern group. The nodal plane angles in red correspond to the inferred rupture plane, withmedian angles of 197◦, 71◦ & 284◦ for strike, dip and rake, respectively.is part of cluster #1. The preferred strike of the event plane computed using principal componentregression on event locations differs by 23◦ from the median nodal plane strike of recovered MTs,whereas the difference in dip falls within the uncertainties.2.5.3 Stress AnalysisStress Axis DirectionsFig. 2.10 A displays the Iniskin earthquake temporal stress inversion results projected on stereonettogether with the fault plane (red), auxiliary plane (blue) and the local E115 model slab surface(green). Note that the principal stress directions remain relatively stable through the sequence. Inparticular, the minimum compression (σ3) axis aligns closely with the dip of the slab indicatingdowndip tension, whereas the intermediate (σ2) axis is slab-normal, and the maximum compression(σ1) axis is horizontal, thus favouring strike-slip motions. There is evidence for a slight variation inthe stress direction in time, corresponding to a minor and temporary rotation σ2 and σ3 around theσ1 axis ∼6 months after the mainshock. The affected events are highlighted with black rims in Fig.2.10 A.26TimeResolutionUncertaintyDiversityMis t0°10°20°30°40°50° 0 36 72108144180Days (Resolution)(Uncertainty, Diversity & Mis t)2015 2016 2017 2018Median Date00.20.40.60.81Stress Ratio (R)Contain MainshockUniaxialCompressionUniaxialTensionA BCFigure 2.10: Temporal stress inversion results for the Iniskin sequence. A) Stereonet showingprincipal stress directions colour-coded in time. The green plane shows the orientationof the slab, the red plane shows mainshock rupture plane orientation and the blue planeshows the auxiliary plane. B) Variation in time of the stress ratio. The time is assignedas the median origin time of all constituent events. The gray shaded interval representsall inversions that include the mainshock. Solutions with black rims in A and B are closeto uniaxial compression. The light blue shaded zone represents the 1-σ vertical error onR over 500 noise realizations. C) Colour-coded time and quality metrics as describedin section 2.4.3.Fig. 2.11 A illustrates the corresponding Anchorage principal-stress-direction variation in time.As for the Iniskin sequence, the σ3 axis indicates downdip tension along the slab surface. However,the σ1 axis is approximately perpendicular to the slab interface, and favours normal faulting. Thestress axis directions vary slightly in time, without obvious systematic variability.Relative Stress MagnitudeA useful means of representing changes in relative stress amplitudes is through the stress ratio R.The stress ratio (also called relative stress magnitude, eigenvalue ratio or shape ratio) characterizesthe relative magnitudes of the principal stress axes and is defined as:R =σ1−σ2σ1−σ3 , (2.1)27Median Date2018Nov30 5pm Dec 1 Dec20.40.50.60.70.80.9Stress Ratio (R)AllSouth ClusterNorth Cluster2016 2018 2019 2020Median Date00.20.40.60.81Stress Ratio (R)2018Nov30 Dec15 20190.30.40.50.60.70.80.91Stress Ratio (R)TimeResolutionUncertaintyDiversityMis t0°10°20°30°40°50° 0 36 72108144180Days (Resolution)(Uncertainty, Diversity & Mis t)ContainMainshockDA BCECompressed5 timesUniaxialCompressionUniaxialTensionFigure 2.11: Temporal stress inversion results for the Anchorage sequence. A) Stereonetshowing principal stress directions colour-coded in time. The green plane displays theslab orientation, and the yellow and pink planes are the mainshock fault and auxiliaryplanes, respectively. B,C & D) Variation in time of the stress ratio R. For each inversion,time is assigned as the median origin time of all constituent events. The gray zone repre-sents those inversions that include the mainshock time. Before November 30, 2018, thetime axis is compressed five times relative to the following. Blue shading envelopingthe stress ratio time series for all events represents the ±1-σ error in R over 500 noiserealisations. Pink/yellow lines on D are the ±1-σ errors in R of the southern/northerngroup respectively. E) Colour-coded time and quality metrics as described in section2.4.3.28where σ1, σ2 and σ3 are the principal stress magnitudes (σ1 ≤ σ2 ≤ σ3), and R can take valuesbetween 0 (uniaxial tension) and 1 (uniaxial compression). Figs 2.10 B and 2.11 B, C & D displaythe Iniskin and Anchorage sequence stress ratios as functions of time R(t).The time series R(t) for the two sequences display different behaviours. The Iniskin stress ratiois characterized by a significant and abrupt drop (∼ 0.72 to 0.23) near the time of the mainshockfollowed by a more gradual increase (∼ 0.23 to 0.83), that is superposed upon a generally decreas-ing trend (∼ 0.8 to 0.6). Conversely, the Anchorage R(t) time series exhibits strong oscillationsfollowing the mainshock. The frequency of oscillations is the highest following 5p.m. on Novem-ber 30, 2018 (median time, ∼30 min/cycle), decreases with time, and may be influenced by thediminishing seismicity rate. Fig. 2.11 D shows R(t) obtained for the northern and southern groups.The individual (north, south) geographic groups display less marked oscillations.2.6 Discussion2.6.1 IniskinThe moment tensor analysis has demonstrated that three of the foreshocks belonging to the Iniskinearthquake possess focal mechanisms similar to that of the mainshock (i.e., they belong to cluster#1 in Fig. 2.5). Moreover, the Mw 4 foreshock that forms part of cluster #2 also has a similarmechanism: if a maximum Kagan angle difference between MTs 26.5◦ had been used for clusteringinstead of the 25◦, clusters #1 and #2 would form a single larger cluster. As shown previously, theangle between the mainshock rupture plane and the local slab orientation is∼76◦ indicating that thefault plane penetrates steeply into the oceanic lithosphere. These observations strongly suggest thatan outer rise fault was reactivated during this earthquake sequence [Peacock, 2001]. Note that thestrike angle of an outer rise fault does not have to match the strike angle of the subducting PacificPlate, because these faults were formed prior to subduction and approximately parallel to the trenchwhose strike may contrast from the strike of the subducting plate at depth.Another interesting observation is that the rake angles of a majority of the events occurring onthe mainshock rupture plane coincide with the direction of the lineation formed by the intersectionbetween the E115 slab model and the rupture plane. Moreover, the USGS mainshock finite faultmodel [U.S. Geological Survey, 2017e] indicates that, during the ∼20 km updip rupture, slip anglechanged orientation from∼68◦ (reverse-slip) at greater depths (earlier times) to∼21◦ (strike-slip) atshallower depths (later times) to parallel the aftershock lineation. The abrupt transition from modeII to mode III failure implies a gradient in stress orientation into the slab. Owing to their location,the majority of aftershocks reflect stress conditions at shallower levels within the slab, and theirconcentration along a lineation that parallels the intersection of rupture plane with the local E115slab geometry implies that the rupture termination was controlled by a stratigraphic horizon withinthe slab.29Since there is more significant uncertainty in the absolute depth of the E115 slab model thanin its dip, the nature of the stratigraphic horizon is ambiguous. In particular, it is conceivable thatthe E115 model is biased deep if the Wadati-Benioff seismicity upon which it is based resides sig-nificantly below the plate interface. Attempts were made to identify converted phases from thehorizon, generated by deeper events, but no compelling examples were found. The two most plau-sible candidates for the horizon are the slab Moho and the top-of-plate interface. One could arguethat the former might be favoured as more likely to maintain its structural integrity to depths of∼100 km. Alternatively, the interpretation that the rupture plane represents a reactivated outer risefault implies that the fault must cut through the oceanic crust and permits an association betweenevent lineation and the top-of-plate interface. In either case, the apparent control of this horizonon faulting suggests an abrupt change in stress or frictional behaviour. Furthermore, the size of theearthquake (Mw 7.1) and USGS finite fault solution [U.S. Geological Survey, 2017e] both imply atleast some slab mantle involvement.Previous studies in Chile [Fuenzalida et al., 2013; Marot et al., 2012] and Japan [Nakajimaet al., 2011] have documented aftershock sequences following intraslab earthquakes along steeplydipping structures inferred to be reactivated outer rise faults. Fuenzalida et al. [2013] interpretedrupture of the 2007 Mw 6.8 Michilla earthquake at∼50 km below northern Chile to occur primarilywithin the oceanic crust based on a plate boundary well defined by numerous thrust events and areceiver function image of the slab. In contrast, Nakajima et al. [2011] provided evidence fromseismicity and tomographic images that the rupture plane of a 2011 M7.1 intraslab event at 66 kmdepth below NE Japan was confined to the lithospheric mantle. At greater depths, as here and thestudy of Marot et al. [2012] of a Mw 5.7 event at 113 km depth below central Chile in 2003, thelocation of the plate boundary and thus a stratigraphic reference for intraslab event rupture planesbecome more ambiguous.Moment tensors for seven other major intermediate-depth earthquakes (60-300 km depth) insouthern Alaska [U.S. Geological Survey, 2021] were examined, revealing that most exhibit strike-slip focal mechanisms similar to that of the Iniskin earthquake, as shown in Fig. A.7 in Appendix A.The nodal planes dipping to the southeast appear to have strike angles that vary in concert with thoseof the slab over along-strike distances of hundreds of kilometres. If the interpretation is correct, thisbehaviour indicates that the genesis of large intraslab earthquakes in this region is controlled byreactivated outer-rise faults. Note that this control is not universal. At some subduction zones,notably along the Cocos Plate in Central America [Warren et al., 2008] and Tonga-Kermadec northof New Zealand [Warren et al., 2007], subhorizontal faults are more common than subvertical faultsat intermediate depths, which are inconsistent with an outer-rise origin.Stress inversions reflect the stress field in a region sampled by input MTs under the assumptionthat the stress field is homogeneous. Accordingly, results shown in Fig. 2.10 B pertains primarilyto the main fault plane of the Iniskin sequence, especially at its upper edge where the majority of30400 600 800 1000 1200 1400Time since 2014 [day]00.10.20.30.40.50.60.70.80.91Stress Ratio25303540Window SizeFigure 2.12: Left: Stress ratio recovered for different inversion window sizes. The verticalgray line indicates the mainshock date. The light blue line is the fit of an er f(w2√Dt)time dependence. Right: Zoom of the stress ratio increase just after the earthquake.events are located and where the slip (and so stress) orientation is locally less heterogeneous thanthe fault plane as a whole. Thus, the temporal stress evolution before and after the mainshock isconsidered, under the assumption that it reflects conditions on shallow portions of a subducted outerrise fault.The stress history R(t) revealed in Fig. 2.10 B has a relatively simple form. Prior to the main-shock, values of R are high (∼ 0.7−0.8), but display an abrupt drop to R∼ 0.3 at the time of rupture,followed by a more gradual return to high values (∼ 0.6−0.8). The drop is defined by stress inver-sions that contain at least 18 (of a total of 30) aftershocks. This behaviour also occurs for differentsliding window lengths (see Fig. 2.12). To interpret this behavior, a recent study by Warren-Smithet al. [2019] is considered. These authors argued that R can be used as a proxy for pore fluid pressurein environments where porosity is fracture controlled and anisotropic, as, for example, in subductedoceanic lithosphere wherein fractures are dominated by bend-related faulting. In such cases, porepressure influences the principal stresses differently as a result of poroelastic effects [Healy, 2012].The case here is more complex than that considered by Warren-Smith et al. [2019] owing to a stressstate at depth that departs significantly from that at the outer rise responsible for the origin of themainshock fault plane and ancillary structures (see Fig. 2.10 A). If σ2 resides in the average frac-ture plane orientation (as in this case), then as pore pressure increases, R is expected to decrease,and vice versa (see Appendix A.4). Observations made by Martı´nez-Garzo´n et al. [2016] that Rvariation along the San Jacinto Fault Zone (strike-slip/normal faulting environment) in California isconsistent with shear wave anisotropy measured by Li et al. [2015], support this relation betweenfluid pressure and R.Accordingly, the time variation of R might be expected to be modeled by pore pressure diffusionaway from an initially pressurized fault zone into the surrounding wall rock. Assuming the fault31zone permeability is the same as that of the wall rock, the decay of pressure at the center of the faultzone isP(0, t) = P0 erf(w2√Dt)(2.2)where P0 is the pressure in the fault zone immediately after the earthquake, w is the half-width ofthe fault zone, and D is the diffusion constant for a fluid-filled porous medium. This model affords areasonable fit to the data as evidenced by the light blue curve in Fig. 2.12 which assumes w = 2.2 mand D = 1.7x10−6 m2/s (a reasonable value for granite). There is a trade-off between w and D. Eachorder of magnitude increase in D corresponds to about a factor of three increase in w. Considerableuncertainty exist in the assignment of time (i.e. each inversion involves a subset of earthquakesunevenly distributed in time). Therefore windows ranging from 25 to 40 events were investigated.All show a relaxation time of about 100 days which constrains w√D. The source of this equation, thechoice of physical parameters, and the fit to data is discussed further in the Appendix A.4.This modeling implies that fluids at high pore pressures existed along the fault at the time ofor immediately following the mainshock. As discussed earlier, dehydration embrittlement as thefundamental triggering mechanism for intermediate earthquakes has fallen into disfavour in recentyears due to the negative volume change (negative Clapeyron slope) associated with the antigoritedehydration reaction at these conditions (-1% at 3.5 GPa) [Ferrand et al., 2017; Gasc et al., 2011;Jung et al., 2004]. However, mechanisms of dehydration assisted stress transfer [Ferrand et al.,2017] that involve fluid production along serpentinized outer rise faults at time of rupture combinedwith production of high pore pressures through compaction at large confining pressures along thefault zone [Brantut et al., 2017] and/or thermal pressurization [Rice, 2006], are consistent with theobservations.2.6.2 AnchorageThe Anchorage sequence displays considerably greater complexity than that of the Iniskin earth-quake, both in terms of the spatial distribution of aftershock hypocentres and in the diversity of theirMTs. It is likely that this complexity is associated with a mainshock that involved simultaneousrupture on two conjugate planes [Guo et al., 2020; Ruppert et al., 2019], one west-dipping and theother east-dipping. This study hypocentral locations show strong agreement with those of Ruppertet al. [2019], with a northern, west-dipping group displaying a clearly defined planar geometry thatmatches one of the mainshock nodal planes and that contrasts with a more diffuse southern group.Guo et al. [2020] argue that slip occurred on the west-dipping plane and the deeper portions of theeast-dipping plane. This rupture process is consistent with the interpretation that the more diffusenature of the southern group hypocenters is the consequence of a strong stress accumulation dueto up-dip rupture termination along a deeper east-dipping fault. The remarkably uniform faultingdocumented by Ruppert et al. [2019] in their analysis of 18 larger (Mw4.1-5.0) aftershocks is further32emphasized in the dominance (209 events of a total 435) of the cluster #1 (Fig. 2.7) that includesthe mainshock. Nonetheless, smaller numbers of events with differently oriented but dominantlynormal mechanisms suggest complexity in fault structures and/or stress orientation at more localscales.The northern event group contains one foreshock within the west-dipping plane with a recov-ered MT similar to the mainshock, supporting the notion of reactivation of a pre-existing fault. Itis conceivable that the mainshock initiated by rupture on the east-dipping plane and triggered con-comitant rupture on the (pre-existing) west-dipping fault 3 km updip of the mainshock hypocentre assuggested by Guo et al. [2020], helping to explain the poor alignment between seismicity and main-shock hypocentre. However, the lack of aftershock seismicity directly associated with the southern,east-dipping plane remains rather puzzling.The orientation of stress axes as determined through stress inversion of MTs from the Anchor-age sequence is relatively stable, revealing downdip tension and slab normal compression. Howeverthe stress ratio R shows marked oscillations in time, in particular, in the days immediately followingthe mainshock when aftershocks were most plentiful (Fig. 2.11). Unlike the Iniskin earthquake,aftershocks appear to be associated with multiple distinct structures that include the primary west-dipping plane, a diffuse cloud of seismicity immediately updip of rupture termination on the primaryeast-dipping plane, and a sparser network of smaller faults distributed about the primary mainshockrupture zones. Given the spatial extent over which faulting is observed and the complexity of main-shock rupture, it is likely that stress is spatially variable, thus contributing to the strong oscillationsin Fig. 2.11 B, C, D. The comparatively smoother temporal dependence exhibited for independentstress inversions involving subgroups of events from the northern and southern groups individu-ally, provides some support for this conjecture. Nonetheless, the absence of an obvious, systematictime dependence in R(t) and the structural complexity of the hypocentral region preclude a simpleanalysis like that carried out above for the Iniskin earthquake.In conclusion, the Anchorage earthquake and associated aftershocks are unlikely to be fullyrepresentative of intraslab events in general. The location of the earthquake near a plate bend atthe boundary of Yakutat terrane and the Pacific Plate is hardly typical. The earthquake is thought tohave ruptured primarily through or next to thickened (15-30 km) oceanic-plateau crust [Worthingtonet al., 2012] and may have initiated within the mantle [Guo et al., 2020; Ruppert et al., 2019], thoughsmaller events farther downdip appear to be confined to the crust [Abers et al., 2013]. In contrast,the Iniskin earthquake almost certainly ruptured through a large tract of oceanic mantle. Thus it isquite plausible that the mechanisms of Anchorage and Iniskin earthquakes have little in commonbeyond a stress regime that is dominated by downdip tension.332.7 ConclusionThe Iniskin and Anchorage foreshock and aftershock sequences, both associated with Mw 7.1 in-traslab earthquakes, differ in their tectonic setting, hypocentral distribution, moment tensor diversityand stress regime. An analysis of these factors has revealed the following interesting features.Most aftershocks of the Iniskin earthquake are located within or near the upper portions of thesoutheast dipping nodal plane identified with the rupture. The fault plane dips steeply into thePacific lithosphere, as represented by the E115 model, and is defined by a linear concentration ofaftershocks at its upper termination parallel to the intersection of the slab with the rupture plane.These observations suggest that the Iniskin earthquake reactivated an outer rise fault and that updiprupture initiated within the oceanic mantle and terminated at a stratigraphic horizon within the slab,likely either the slab Moho or plate interface. Inversion of moment tensors for time-dependentstress reveals a relatively stable stress orientation marked by downdip tension; however, the stressratio undergoes an abrupt decrease immediately following the mainshock succeeded by a moregradual recovery. Under the assumption that fracturing within the Pacific plate is anisotropic, thisbehaviour is interpreted to manifest diffusion of pore fluids and pore pressures, generated throughmetamorphic dehydration reactions, away from the fault zone. This interpretation is consistent withdehydration-driven stress transfer as a triggering mechanism.The Anchorage sequence took place within a less typical intraslab environment involving shal-low subduction of an oceanic plateau. The aftershock distribution is more diffuse and falls largelyinto two groups: a southern group associated with the east dipping mainshock nodal plane, and anorthern, more tightly aligned group associated with the west dipping nodal plane. Like previousauthors, this distribution is interpreted to manifest a mainshock rupture on dual, conjugate ruptureplanes where the northern aftershocks define one plane, and the southern aftershocks lie updip of thesecond plane that is otherwise undefined by aftershocks. Moment tensor diversity is significantlygreater than that for the Iniskin sequence, though almost 50 percent of aftershocks share a focalmechanism similar to that of the mainshock. Stress inversions indicate a stable orientation markedby downdip tension and slab normal compression. In contrast, stress ratio is variable without ob-vious, systematic trends perhaps manifesting heterogeneous stress adjustment following mainshockrupture on conjugate faults.34Chapter 3Relative Moment Tensor Inversion forMicroseismicity: Application toClustered Earthquakes in the CascadiaForearc3.1 Chapter SummaryThe relative abundance of small earthquakes affords significant opportunities for improved un-derstanding of regional seismotectonics; however, determining moment tensors for such eventsrecorded on regional networks is complicated by low signal-to-noise ratios, sparse station sam-pling, and complex wave propagation at short periods. Built upon previous work, a multiple-event,simultaneous moment tensor inversion scheme for small earthquakes is designed, which employsconstraints from P-wave polarities, relative amplitudes of P- and S-waves recorded at commonstations, and local magnitude estimates. This method does not require a priori knowledge of areference moment tensor. High-fidelity polarity and relative amplitude data are recovered usingprincipal component decomposition of clustered-event waveforms. These data are employed withina multi-stage iterative framework to invert for moment tensors and incorporate local magnitudeinformation. Synthetic examples employing as few as four high-quality and spatially-distributedstations yield accurate moment tensor estimates. The approach is demonstrated on a cluster ofseismicity near San Juan Island, Washington, USA, within the Cascadia forearc. The results areconsistent with previous characterization of the local stress regime and support an interpretation ofswarm behaviour resulting from migration of fluids originating from dehydration of the subductingJuan de Fuca plate.353.2 IntroductionMoment tensors (MTs) characterize the size and radiation pattern of seismic events approximated aspoint sources. Accordingly, they offer fundamental information on source type (e.g., shear, implo-sion/explosion or hybrid mechanisms) and, in particular, the style (normal, reverse, strike-slip) ofshear faulting that describes most earthquakes. When moment tensors from many events are avail-able for a given region, they reveal the complexity of fault systems, magnitude-frequency statisticsand the nature of the tectonic stress regime, and consequently are useful in seismic hazard assess-ment.Moment tensors are routinely and automatically computed for larger events, usually MW ≥ 4[e.g., Ekstro¨m et al., 2012]. The predominant approaches are based on centroid moment tensors(CMTs) that involve the modelling of long-period (≥∼ 20s) waveforms using Green’s functionsthat are typically computed for radially stratified Earth models. Events smaller than MW 4 generallypossess low signal-to-noise ratios in the long-period band, and are poorly modelled using 1-D Earthmodels, a tendency that becomes increasingly pronounced with decreasing magnitude. This cir-cumstance is unfortunate given that smaller events occur with greater frequency, such that accuratemodelling of their moment tensors could provide improved characterization of fault systems andstress variability.The classic approach to source characterization of small events involves the recovery of the fo-cal mechanism, a binary representation of the P-wave radiation pattern assuming a shear-faulting(double-couple) source [e.g., Shearer, 2009]. First-motion P-polarities are identified for a givenevent across an array of stations encompassing diverse azimuth and take-off angles. Thereafter,two orthogonal planes bounding quadrants of positive and negative polarity are determined, oftenemploying grid search methods. Although straightforward in principle, P-wave polarity assignmentcan be ambiguous and time-consuming in practice. Furthermore, the event-station geometry con-trols the uncertainties of the recovered focal mechanism and is frequently insufficient to constraina unique solution. P-polarities are sometimes augmented with S-to-P amplitude ratios to providefurther constraints on focal mechanisms. Hardebeck and Shearer [2003] note only relatively mod-est improvements in focal mechanism recovery, likely due to large measurement uncertainties (seehowever Shelly et al. [2022]), and increased sensitivity to the typically poorly constrained P andS-velocity and attenuation models.Consideration of multiple events in spatial proximity affords additional opportunities for con-straints on moment tensors using relative P and S amplitudes [Dahm, 1996] that are analogous tothose afforded by traveltime double differences in the determination of relative hypocenters [Wald-hauser and Ellsworth, 2000]. In particular, the use of a common mode normalizes propagationeffects such that the relative amplitudes are sensitive primarily to the differences in source strengthand orientation. Plourde and Bostock [2019] have extended Dahm’s original approach by introduc-ing a general treatment for inclusion of S-waves and employing principal component analysis in the36measurement of precise relative amplitudes. Their approach has been applied to yield improved un-derstanding of the mechanics of deep Wadati-Benioff seismicity in Alaska [[Plourde and Bostock,2019] and Chapter 2 of this thesis [Drolet et al., 2022]]. Other recent studies incorporating relativeamplitude measurements include the work of Cheng et al. [2023] and Zhang [2023].This study presents an improved version of the relative moment tensor inversion method ofPlourde and Bostock [2019] tailored for small earthquakes, that incorporates P-wave first motionpolarities and removes the requirement of a known reference MT. In addition, our method solves forrelative or absolute magnitudes depending on availability of local magnitudes. The following sec-tion provides an overview of the general relative moment tensor method and a detailed descriptionof our implementation. Synthetic examples are presented to gauge the performance of our method,and apply it a dense cluster of seismicity in the northern Cascadia forearc.3.3 MethodologyThis section begins by introducing the general assumptions underpinning body-wave relative mo-ment tensor inversion and then follows with a detailed description of the implementation.3.3.1 AssumptionsThe usual assumptions made in developing relative moment tensor inversion using body waves aredepicted in Fig. 3.1 and rationalized as follows: 1) Ray theory governs wave propagation, implyingthat direct arrivals can be traced back to the focal sphere. This information is utilized for take-off,azimuth, and incident angle calculations. 2) The immediate source region representing the focalsphere is assumed to be homogeneous and isotropic, ensuring no perturbation in seismic radiationpatterns due to heterogeneity or anisotropy. This assumption allows the use of standard far-fieldexpressions for the radiation of direct body waves. 3) The Green’s function is the same for pairs ortriples of events, requiring events to be clustered such that the event-station distance is much greaterthan the inter-event distance. This requirement ensures the normalization of propagation effects(Green’s function) in event comparisons. 4) Direct P- and S-waves are assumed to be isolatedfrom other arrivals within short-duration windows, to avoid contamination by multipathing. Thisassumption facilitates precise relative amplitude calculation and, for P-waves, allows for improvedpolarity determination. 5) It is assumed that slip is unidirectional, and that the source time functionis normalized to a common bandlimited signature for all events. This assumption involves low-passfiltering seismograms below their respective source corner frequencies, effectively removing thetemporal source signature in event comparisons.37Figure 3.1: Schematic diagram summarizing the assumptions required to apply the relativemoment tensor inversion.3.3.2 Relative ConstraintsUnder the assumptions laid out above, the relative amplitudes of P or S-waves of multiple eventsmeasured at a given receiver can be directly related to their relative amplitudes on the focal spherethat surrounds their common source location (see Plourde and Bostock [2019], for more details). Inparticular, for P-waves with polarization constrained to a single direction (e.g., parallel to the rayfor isotropic media), the displacement waveforms ui(x, t) measured at a receiver at x in direction ifor a pair of events a, b with differing moment tensors can be related asAabubi (x, t)−uai (x, t) = 0, (3.1)while at the focal sphere, a similar relation holdsAabRabγbi γbj γbk Mbjk− γai γaj γak Majk = 0, (3.2)where γi represents the direction cosine (and so incorporates the take-off angle and azimuth) at thefocal sphere, Rab is a geometric correction to compensate for the difference in source locations(necessarily close to 1) and e.g., Mai j is the moment tensor for source a. Each relative P-wave mea-surement thus provides a constraint on the moment tensors of two events. S-wave polarizations areconstrained to lie within a plane and hence exhibit 2 degrees of freedom (vs 1 for P). Accordingly,38relative amplitude waves across triples of events c,d,e at the receiver are related as:Bcdeudi (x, t)+Bceduei (x, t)−uci (x, t) = 0, (3.3)while at the focal sphere, the relation isBcdeRcd(δi j− γdi γdj )γdk Mdjk +BcedRce(δi j− γei γej )γek Mejk− (δi j− γci γcj )γck Mcjk = 0. (3.4)Note that (3.4) yields 2 independent constraints; hence each S-wave triple places 2 constraints on 3moment tensors.The relative amplitude coefficients Aab, Bcde, Bced defined above can be accurately determinedusing principal component analysis as described in Plourde and Bostock [2019], and the highestfidelity results are expected for 3-component waveforms. Equations (3.2) and (3.4) represent linearconstraints on the MT elements incorporated within a large simultaneous system. In the ideal case,where all possible pairs and triples are formed for N events and K stations, there are N(N−1)K/2possible P-wave constraints and 2(N(N − 1)(N − 2)K/6) possible S-wave constraints on 6N un-known MT elements. However, these equations alone create an under-determined linear systembecause the only information used is relative and so a) the signs of the MTs are not uniquely deter-mined (i.e., the negative of a solution vector of MTs remains a solution vector), and b) the absolutescale (or equivalently, the absolute scalar moments) is not determined. Therefore, extra informationis required to solve the system. This information can be supplied in the form of one or more apriori known, high-quality reference MTs [e.g., [Plourde and Bostock, 2019] and Chapter 2 of thisthesis [Drolet et al., 2022]]. However, accurate MTs are typically not available for small events andpoorly constrained estimates may lead to biased solutions [Dahm, 1993]. Instead, P-polarities andlocal magnitudes are incorporated to constrain the polarity of the radiation pattern and scaling ofthe seismic moments, respectively.3.3.3 The AlgorithmIn order to recover reliable MTs using relative MT inversions, the 9-step framework depicted inFig. 3.2 and summarized below is adopted.Step 1 - Waveform ProcessingThe first step is to collect all displacement waveforms, preferably with instrument response removedin the case that individual stations have experienced changes in equipment over time. Ideally wave-forms should be recorded at 3-component stations for S-waves and at 1- or 3-component stationsfor P-waves. One must then determine the optimal frequency band for waveform comparison overwhich signal-to-noise ratio is maximized. Morever, the chosen band must lie below the event corner39Figure 3.2: Summary of the relative moment tensor algorithm. Steps, inputs, and output arein yellow, green, and red boxes. Detailed descriptions of each step are given in section3.3.3.frequency to remove details of rupture and so validate the point-source approximation (assumption5). Waveforms with a signal-to-noise ratio below 1 are discarded.Next, waveforms must be accurately aligned. This step is crucial since misalignment can sig-nificantly impact the relative amplitude measurement, and, for example, a half-cycle skip will yieldan incorrect sign. For this purpose, a multi-channel alignment technique is employed, extendingthe work of VanDecar and Crosson [1990] to S-waves generated by sources with variable MTs[Bostock et al., 2021] by exploiting principal component analysis. The alignment is systematicallyperformed for all events recorded at a given station (or instrument, if accurate response removal isnot achieved) for a given (P or S) mode. Only waveforms with accurate alignment conforming tothe rank 1 (P) or rank 2 (S) principal component description are accepted for the next stage. Whenthis alignment technique is applied to a large number of events, it is generally straightforward tovisually discern and correct misaligned events.40Step 2 - Relative Amplitude MeasurementOnly pairs or triplets of events in sufficiently close proximity, and so sharing a similar ray path, willhave their relative amplitude computed. A shared path value (ψ) between events a and b is definedas:ψ = 1−2∥xa−xb∥∥xa+xb∥ (3.5)where e.g., xa is the source-to-receiver vector for event a. Event comparisons with ψ values belowa default threshold of 0.6 are discarded.To calculate relative amplitudes, principal components are exploited to mitigate the effect ofnoise by comparing only the common signal between pairs or triplets of waveforms. Waveformsare bandpass-filtered to the highest optimal frequency band found in the previous step within thepair or triplet. Details on how to calculate the relative amplitude are found in Plourde and Bostock[2019]. The misfit between the reconstructed waveform using the relative amplitude and the originalwaveform is calculated and assigned as a quality value. Relative amplitudes associated with highmisfits (default greater than 0.6) are discarded. Furthermore, misfit values are used to weight theconfidence of the relative amplitudes in step 6.Use of 3-component waveforms for determining P-wave relative amplitudes is recommended,as they enable minor rotation adjustment of one waveform a pair into the orientation of the other,thereby improving measurement accuracy. However, when events are extremely close to one anothercompared to the event-station distance, the improvement is negligible; thus, it is possible to usesingle-component stations. Note that 3-component recordings are still more relevant for S-waverelative amplitudes due to the two degrees of freedom in their polarizations.While it is possible to employ either P- or S-wave relative amplitudes alone within relative MTinversion, it is preferable to use both modes when the data permit. In particular, in environmentscharacterized by high levels of random noise affecting both P- and S-wave phases equally, thealgorithm demonstrates superior performance when both phases are included.Step 3 - P Polarity MeasurementAs in traditional focal mechanism determination, the acquisition of P-wave first motion polaritypresents challenges for small regional events. Principal component analysis is incorporated to im-prove polarity estimation. Polarities are measured on the Z-component as it typically exhibits themost pronounced first motion, although, in principle, the E- or N-components can also be employed.For each station, a waveform stack (X) is created by weighting the waveform of each event i by itsfirst principal component expansion coefficient Yi1 as follows, using the singular value decomposi-tion:W = UΛVT (3.6)41Yi1 =Vi1Λ11, (3.7)Xk =∑iYi1Wˆki. (3.8)Here the matrix W in (3.6) contains the filtered, amplitude-normalized and aligned waveform sec-tion, whereas Wˆki in (3.8) is its unfiltered (or filtered to as wide a band as signal levels allow)counterpart. Each column represents a different event i, whereas the rows represent discrete timesamples k. A band-pass (minimum-phase) Bessel filter is employed, reducing phase shift relativeto e.g., the corresponding Butterworth filter. The expansion coefficient Yi1 implicitly includes thesign of the waveform for event i such that the weighted waveforms on the right-hand side of (3.8)constructively interfere, yielding a high SNR stacked trace Xk, whereas the inclusion of unfilteredwaveforms reduces potential ambiguity in polarity assignment introduced by band limitation. Theuser then only has to visually assign the first motion polarity of the stacked trace Xk to uniquelydefine the polarities of all individual events i that follow from their expansion coefficients, i.e.polarity{Wki}= sign{Vi1}×polarity{Xk} (3.9)Step 4 - Event CullAfter these processing steps, some events may not have sufficient constraints to fulfill the qual-ity control criteria and are, therefore, removed. More specifically, each event is required to havemeasured relative amplitudes for at least four different stations that represent independent source-receiver geometries. As a rule of thumb, it is recommended that the minimum number of stationsbe equal to the number of independent MT elements to be resolved (4 for DC, 5 for isotropic, and 6for full MT) when noise levels are low.Step 5 - Geometrical FactorsThe linear equations for relative amplitudes (3.2, 3.4) require the direction cosine and, if desired,geometric corrections for source separation, which are obtained using the event and station locationsassuming a velocity model. The directions cosines specify the ray exit points at the focal sphere fora given velocity model under the assumption 1.Step 6 - Linear System AssemblyA linear system is constructed, incorporating equations (3.2) and (3.4) for all relevant event com-binations. As the number of constraints increases rapidly with the number of events (O(N2) for P,O(N3) for S), but the constraints only involve 2 or 3 MTs per equation, sparse matrices are used tostore information efficiently. The resultant system is represented as Am = b, where matrix A in-corporates misfit weighted equations pertaining to both relative and, if available, MTs for reference42events. Vector b consists primarily of zero entries (with the exception of equations linked to the ref-erence event(s) if available), and vector m encompasses all independent moment tensor elements.The length of vector m depends on whether the recovered moment tensors are to be full (i.e. includeisotropic component) or deviatoric, resulting in dimensions of either 6 or 5 times the number ofevents N, respectively. A comprehensive description of the system construction methodology canbe found in Plourde and Bostock [2019].Step 7 - System conditioningLarge linear systems are difficult to solve accurately if not properly conditioned. Therefore, scalingof A and b is performed to reduce numerical error. Elements are scaled by events and by equations.In one or multiple reference events are included within the system, magnitude scaling is applied tothe reference event equations. The event and, as relevant, reference event scaling factors are retainedfor reversing scaling of the MT elements after solution of the linear system.Step 8 - System SolutionDifferent strategies are employed to solve the linear system dependent on the nature of the MTpolarity (i.e. reference MT vs. P polarities) and model (DC vs. deviatoric / full) constraints. Thereare two cases that require attention for solution. First, the absence or inclusion of reference event(s)will dictate if the linear system is under-determined or over-determined. Secondly, if DC solutionsare desired, a non-linear constraint is added to the system, and a gradient descent optimizationapproach is adopted. The solution algorithm is structured into four distinct stages to accommodatethe different cases. A summary of these stages is presented in Fig. 3.3, while detailed explanationsof each stage are provided in subsequent subsections. For the present, tensile MTs without CLVDcomponent are not solved for, but note that this could be accomplished by modifying stage 2 withthe projected Newton descent technique used by Plourde and Bostock [2019].Step 8 - Stage 1: Model InitializationAt this stage, an initial moment tensor guess solution (m /0) is introduced for the linear system,which plays a pivotal role when employing gradient-descent optimization as required to recoverDC MTs. A good initial guess will significantly improve the convergence of the algorithm and sodiminish computation time. An intuitive way of securing an m /0 would be to solve the linear system(Am /0 = b). In instances where one or more reference events are included, this approach is possiblebecause the system Am /0 = b is over-determined.However, when reference events are not incorporated, an under-determined linear system isencounter. Accordingly, the linear system is temporarily augmented by appending an equation:1000xTm /0 = mk = 1000, where x represents a unit column vector. The new equation forces the43Figure 3.3: Summary of the four solving stages. Each of the stages is described in the maintext.kth element of m /0 to a value of 1 and therefore fixes the sign of all MTs. The new system isAkmk/0 = bk. It is important to note that this procedure does not guarantee the correct polarity of thesolution mk/0 as reflected by the measured P polarities; thus, a check is performed. If less than 50%44of the P polarities agree with those predicted by mk/0; the sign of mk/0 is reversed: (mk/0 = −mk/0). Itis further recommended to discarding mk/0 when the P polarity agreement is less than 60% to avoidpolarity reversals during the iterative solution. In addition, low polarity agreement may signifythat the corresponding element mkk is either near zero or that some polarity information is incorrect.This manner of initial guess selection should force m /0 closer to the eventual desired solution than apurely random selection. Consequently, the linear system Akmk/0 = bk is resolved for each k, alwaysrequiring that mk/0 satisfy at least 60% of the P polarities constraints, thereby generating MN m /0initial solutions (where M = 5,6 for deviatoric and full moment tensors, respectively, and N is thenumber of events) to be used in subsequent stages. Note that in the further steps, the augmentedsystem is not used.Step 8 - Stage 2: Iterative ProcessAt this stage, it is determined whether an iterative approach (gradient descent) is necessary for theinversion process. In particular, if the desired MTs are full or deviatoric, an iterative approach isunnecessary, and m is set equal to m /0. Conversely, if DC MTs are desired, a gradient descent ap-proach is employed to retrieve the optimal m. This requires defining a non-linear objective functionsatisfying the relative amplitude and the DC constraints.DC solution implies that a MT is characterized by only four independent elements, in contrast tothe six independent elements that must be specified for a full MT. One degree of freedom is lost byrequiring an non-isotropic source (no volume change), which is enforced mathematically by settingthe MT trace to 0. This condition is addressed during the construction of the relevant matrices (step6). Another degree of freedom is lost by prohibiting a compensated-linear-vector-dipole (CLVD)component, which is achieved by forcing one of the MT eigenvalues to 0.To effectively apply this latter constraint, the star norm, also referred to as the nuclear norm, isemployed. The star norm is defined as the sum of the singular values, mathematically representedas:∥M∥⋆ =3∑i=1σi = trace(√M∗M) (3.10)Here, σi is the ith singular value of M. For a square symmetric matrix such as M, the singularvalues correspond to the absolute values of the eigenvalues. Since σ1 ≥ σ2 ≥ σ3 and σ1 = σ2 for aDC source, the star norm is minimized when σ3 = 0. In the more general case of the linear system,the goal is to minimize the sum of the star norm of the MTs for all events. Accordingly, a modifiedsolution for the objective function is defined as the sum of a quadratic misfit term and a penalty termthat enforces DC MTs for all sources:f (m) = ∥Am−b∥2︸ ︷︷ ︸relative constraints+α#events∑j∥M j∥⋆︸ ︷︷ ︸DC constraints(3.11)45where α is a tradeoff parameter that weighs the importance of the DC constraint over the relativeconstraints. By default, the value is set to 10.During the gradient descent, m is updated as ml+1 = ml − µ δm at each iteration l with itsgradient δm calculated with respect to the objective function f (m) and where µ quantifies thestep size in the gradient direction. The gradient is calculated using the autograd function fromthe python package PyTorch [Paszke et al., 2019]. To improve convergence, the step size µ variesthrough the iterations. For each iteration, µ is divided by two until the penalty term, computed fromml+1 = ml − µ δm, decreases. If this criterion is met on the first try, µ is increased by 1.5 for thenext iteration.The iteration process stops when one of the following conditions is met: 1) the DC solution isachieved within an error margin, 2) the step size µ becomes too small, 3) a maximum number of it-erations is reached or 4) the improvement for the past 200 iterations is very small. Upon completion,the solution is converted to the closest perfect DC solution. In the case of condition 3), the solutioncan not be trusted and thus is discarded. This situation is generally due to noisy data rendering thesystem ill-conditioned.Note that the star norm is preferred over the Langrangier multiplier method used by Plourdeand Bostock [2019] to set the determinant of each event to 0. This star norm technique is faster andcan be used when no reference events are given (when A is under-determined). As a final remark,when solving for DC solutions using reference events, all reference events should be supplied asDC; otherwise, instabilities may arise.This stage is applied to each mk/0, and is easily performed in parallel.Step 8 - Stage 3: Final SolutionIn the presence of a reference event(s) only one solution, the final solution, makes it to this stage.However, in the absence of a reference event, multiple solutions are evaluated corresponding to thedifferent initial guesses mk/0. Most solutions are local minima that characterize the presence of noiseand, for DC solving, the non-linear optimization function. These solutions can be assessed by plot-ting their percentage of incorrect polarities vs. the relative constraint misfit ||Amk/0||, the two valuesbeing minimized. The optimal solution should lie near the lower left corner of this plot, where boththe number of incorrect polarities and the misfit are near zero. The solution with the lowest relativeconstraint misfit (||Am||) within 95-percentile solutions with the lowest percentage of incorrect po-larities becomes the final solution. Specification of the misfit percentile may vary depending on thenumber of events, so the lowest misfit solution is chosen among a range of solutions.Step 8 - Stage 4: Quality ControlThis stage is only applied in the absence of reference events. Here, the stable event solutions areidentified by selecting a subset of possible solutions from the relative amplitude misfit vs. incorrect46polarity plot. The most likely solution is selected by taking the bottom left solutions delimited by 95-percentile (or chosen-percentile) solutions with the most polarity agreement and the 60-percentilesolutions with the lowest misfit. Subsequently, the Kagan angles between these solutions and thefinal solution are calculated if the solutions are DC. The Kagan angle is a metric that quantifies theminimum angle necessary to rotate one DC focal mechanism into another orientation. Its numericalvalue spans a range from 0◦ for absolute agreement between the two mechanisms, to 120◦ signify-ing complete disagreement [Kagan, 1991]. Typically, when the Kagan angle falls below 60◦, it isregarded as indicative of a good correspondence [Pondrelli et al., 2006]. The Kagan angle standarddeviation for each event serves as an uncertainty proxy, representing a measure of solution stability.Events are considered ‘stable’ if their standard deviation is below 20◦. Solutions with a standarddeviation between 20-30◦ are classified as ‘likely’, since there uncertainties are larger. Events witha standard deviation above 30◦ are considered ‘bad’ and are discarded for further steps which im-prove the magnitude scaling accuracy. For non-DC solutions, the cross-correlation value betweenthe six independent elements of the normalized MTs is used as a proxy agreement between MTradiation patterns, instead of the Kagan angle. A value of 1 would represent the same radiationpattern, whereas a value of -1 would constitute the MT with opposite polarity. Standard deviationsof cross-correlation values below 0.15 are considered ‘stable’, between 0.15 and 0.2 ‘likely’ andabove 0.2 ‘bad’.Step 9 - Solution ScalingThe MTs must be rescaled after solution, by the inverse of the scaling factors applied in step 7. Inthe case of the use of a reference event(s), the recovered MTs are now properly referenced to theirseismic moments M0. In the absence of a reference event(s), all MTs values are relative, but byconsidering magnitude information can be rescaled to their (approximately) correct M0 values.In most cases, only local magnitudes will be provided. Local magnitudes (ML) can be convertedto moment magnitudes (MW ) using empirical relations [e.g., Shearer et al., 2006]:MW ={ [(23)( 10.96)]ML+0.917 if ML < 3.0ML if ML ≥ 3.0.(3.12)Note that this conversion is based on California earthquakes, and more appropriate conversions maybe derived for other regions. Moment magnitudes are then converted to seismic moments using thedefining formula from Hanks and Kanamori [1979]:M0 = 1032 MW+9.1. (3.13)When a catalogue of seismic magnitudes is available, a square system is assembled to recoveroptimal scalar moments for all events. The linear system comprises the constraints on M0i for each47event i from the relative MT inversion along with scalar moment estimates MCAT0i based on cataloguemagnitudes which need not be available for all N events. For example, the system for 4 events withcatalogue magnitudes provided for events 1 and 2 is written as:1 −1 0 01 0 −1 01 0 0 −11 1 0 0log(M01)log(M02)log(M03)log(M04)=log(M01/M02)log(M01/M03)log(M01/M04)log(MCAT01 )+ log(MCAT02 ) . (3.14)This system ensures that relative scalar moments ratios determined from the relative MT solutionare honored whereas local magnitude information supplies the optimal scaling. Since the matrixin (3.14) is square (i.e. not overdetermined), no weighting of the equations is necessary. Althoughother scalar moment reconstructions are possible [e.g., Hayward and Bostock, 2017], synthetic testshave shown that the above reconstruction yields superior accuracy across the different schemesconsidered.3.4 Synthetic ExamplesTo demonstrate the efficacy of the algorithm, its performance is showcased through synthetic ex-amples. This section focuses on scenarios without reference events since scenarios with referenceevents have already been examined in Plourde and Bostock [2019] and Chapter 2 of this thesis[Drolet et al., 2022].The examples comprise a cluster of 20 earthquakes randomly located within an area of 25km2with depths varying between 20 and 25km. The event area is enclosed within a 120km x 120kmarea that includes two station configurations. Four/six stations are randomly distributed within thislarger area and employed to study DC/full MT recovery, respectively. The MTs for each event arerandomly defined with logarithmically distributed magnitudes between Mw 1 and 3. For simplic-ity, a 1D homogeneous velocity model is employed so that relevant angular quantities are easilydefined. Synthetic amplitudes are contaminated by up to ±20% of random noise, and 10% of thepolarities are reversed from their correct values and polarities close to nodal planes are removed.Only 10% of the total possible P- and S-wave relative amplitude information (76/114 P- and 456/684S-comparisons for DC/full MT) is included within the data set to simulate a real-case scenario. Thelargest event magnitude is assumed known and used to scale MTs to their scalar moments.Fig. 3.4 displays the MT inversion results for the setup described above for a typical realization.Since no reference event is employed, many candidate solutions are produced. The final solution isthe solution with the lowest relative amplitude misfit within 80-percentile solutions with the lowestpercentage of incorrect polarities. Note that this figure does not display the best synthetic resultsbut, rather, more typical results demonstrating a typical range in uncertainty. The magnitudes andradiation patterns are in good agreement for both scenarios (DC/full MTs).48Figure 3.4: Synthetic examples as described in the main text. A, B, C) show the event-stationdistribution. Red triangles are stations used for both inversions, whereas green trianglesare only used for the full MT inversion. Moment tensors are plotted in lower and southernhemisphere representation for B and C) respectively, and scale with their magnitude. Thefull MT solutions are displayed; the closest DC solutions are used for the DC inversion.The lower panel shows the DC MT (D) and full MT (E) solutions. The true solution isblue, whereas the red contours show the recovered solution. True magnitude is labelledin black and recovered magnitude in gray, below each MT. Dots are polarities used in theinversion. Filled/empty dots are associated with up/down first motion; dots highlightedin yellow are incorrect polarities used for the inversion. For the DC or full solution, theKagan angle or cross-correlation between the true and the recovered solutions is indicatedbelow the magnitude.49To demonstrate the robustness of the algorithm, the previous scenario is repeated 1,000 timeswith different random number seeds. Fig. 3.5 displays the results. Step 8.4 (Quality Control)was not applied to simplify comparison so that unstable events can be part of the final solution.Generally, a larger error in radiation pattern accompanies the largest error in scalar moment recov-ery. There is no obvious bias toward under/over-estimating magnitude for DC solutions. However,a strong bias toward overestimating the moment magnitude is noticeable for full MT solutions.No explanation for this bias have been found yet. Ninety percent of the median errors are below16.7◦/0.88 for DC/full MT solutions, respectively, indicating good agreement. Almost all medianmagnitude errors lie within ±0.1 for DC and full MTs. Figs B.1 and B.2 are similar to Fig. 3.5,but with 2 extra stations added and with 10 added events, respectively. As expected, increasing thenumber of stations significantly improves the accuracy of the recovered MTs. The spatial coveragealso affects results: improved coverage yields more accurate solutions. Increasing the number ofevents does not appear to significantly affect the quality of the results while the relative amplitudeerrors remain the same.Figure 3.5: Synthetic results for 1,000 realizations DC MTs (left) and full MTs (right). Notethat the normalized cross-correlation axis on the right figure has been flipped for bettervisual comparison. The magnitude difference is calculated as the truth minus the recov-ered moment magnitude.3.5 Real Case Example: San Juan ClusterIn this section, the moment tensor inversion algorithm is applied to one of the largest forearc seis-micity concentrations in northern Cascadia. The earthquake cluster is located near the Canada-U.S. border along the northeast coast of the San Juan Islands Archipelago, Washington State (see50Fig. 3.6). A suite of more than 700 crustal earthquakes (maximum depth 30km) with local mag-nitudes ≤ 3.2 recorded between January 2000 and December 2020 which were well recorded byregional stations is examined. More, the largest event ever recorded within this region was an ML3.4 earthquake on September 24, 1997 [U.S. Geological Survey, 2017f].Figure 3.6: Seismic stations used for MT inversion of local seismicity near San Juan Island.Blue, red and purple triangles are stations that contributed relative amplitudes, polarities,and both, respectively. Black dots are earthquakes from Merrill et al. [2020, 2022]. Thegray dashed line is the Canada-U.S. border. The purple lines are the Moho contour ofthe subducting Juan a Fuca Plate from Bloch et al. [2023]. The red box defines the studyarea.The upper crust of the San Juan Island Archipelago (comprising San Juan, Orcas, Lopez and51numerous smaller islands) is composed of a 5-10 km thick mid-Cretaceous sequence of nappesand thrust faults referred to as the San Juan-Cascade Thrust System. Brown [2012] suggests thatthe nappes were formed through initial accretion and subduction-zone metamorphism south of thecurrent location, followed by uplift, exhumation, orogen-parallel northward transport as a forearsliver component, and ultimately obduction at the current site over the truncated southern end ofWrangellia and the Coast Plutonic Complex. However, the nature of the deeper crust beneath thenappes and within which the majority of the seismicity likely occurs, is unconstrained.Fig. 3.7 displays the event timeline and identifies the most pronounced sequences by their dates.Half of these sequences are associated with an initial mainshock, whereas the other half behavemore like earthquake swarms (a group of events occurring closely in time and space without asingle dominant shock). Swarms have often been interpreted to be caused by fluid [Vidale andShearer, 2006] or magmatic activity [Hill, 1977].Figure 3.7: Temporal distribution of San Juan Island seismicity. Each vertical bar representsa seismic event; a darker region signifies a shorter recurrence interval. Middle panelin red displays events comprising the San Juan cluster; more pronounced sequences areidentified by their starting date. Red dots are events with recovered magnitudes from theMT inversion (see axis at right for MW scale). Top/bottom panels in black/blue are thetimelines for earthquakes/tremors surrounding the San Juan seismicity. The surroundingarea is 9 times larger than the study region and is identified in Fig. 3.13. Backgroundseismicity does not include events from the San Juan cluster. Tremors between 2003and 2007 are high accuracy (waveform correlation) locations and times from Armbrusteret al. [2014] whereas tremors after 2009 are from Pacific Northwest Seismic Network[2009-2021] (PNSN) (envelope correlation) catalogue. Note that the deployment of thePolaris array (2002-02-26 to 2006-07-04) considerably increased the number of stationsin the region and hence decreased earthquake detection thresholds over the period.523.5.1 Data Set and ProcessingThe event catalogue from Merrill et al. [2022] and stations within 200 km of San Juan Island(Fig. 3.6) from diverse networks were utilized. Waveforms were downloaded from data centers ofthe Incorporated Research Institutions for Seismology (IRIS) Center and Natural Resources Canada.Phase picks for stations in Canada were obtained from Merrill et al. [2022] as a mix of manuallyand automatically detected phases. All picks were visually checked and revised as needed. Phasesat U.S. stations were manually picked and represent approximately half of the arrival time cata-logue. First motion P-wave polarities were primarily determined using the technique described instep 3, but some were also manually picked if a clear and consistent polarity was observed across4 frequency bands ([1,6],[2,7],[3,8],[0.5,8]Hz). If inconsistency was noted between manual andautomatic polarities, those measurements were discarded.Earthquake RelocationAfter augmenting the arrival-time phase catalogue, earthquakes were relocated utilizing the Double-Difference Hypocenter Locations (hypoDD) program developed by Waldhauser and Ellsworth [2000].For this relocation, a 1-D velocity model as the average of the Merrill et al. [2022] 3-D velocitymodel centred on the cluster and extending out to a radius of approximately 50 km (see Fig. B.3 forprevious location) was generated. The 549 relocated earthquakes are displayed in Fig. 3.8.HypoDD relocation places most of the earthquakes at depths between 10-23 km within reason-ably tightly constrained subclusters that also tend to be clustered in time. Very few events are locatedabove 10km within the nappe structures described by Brown [2012]. Furthermore, the seismicity inaggregate displays a northward dip which was previously noted by Balfour et al. [2012], whereasthe east-west profile displays some clusters as near-vertical planar structures.Moment Tensor InversionThe steps outlined in section 3.3.3 was followed using the default parameters. Take-off and azimuthangles were obtained from the hypoDD outputs. Missing values from the hypoDD file (especiallyfor P-wave first motion) were calculated using the 1-D velocity model and the python packageobspy.taup [Crotwell et al., 1999].A total of 350 events passed all quality control steps and were admitted to the double-coupleMT inversion. A subset of 20 solutions was selected based on the 95th percentile polarity and 60th-percentile of amplitude constraint misfits, to assess uncertainty in the final solution of 350 momenttensors as described in step 8.4 (Quality Control). Fig. 3.9 displays the solution distribution inrelative residual versus polarity space. There are 171 stable events, 134 likely events and 45 badevents.The final solution agrees with 83% (1624/1948) of the polarity measurements and possesses an53Figure 3.8: Earthquake distribution color-coded by year. A) Top view with a black cross in-dicating the cluster center. B) East-west cross-section centred on the cluster. C) South-north cross-section centred on the cluster.54Figure 3.9: Optimal solution selection. Blue dots represent all recovered solutions. Green dotsare the subset used for the uncertainty assessment based on 95th-polarity and 65th-linearconstraint percentiles. The red dot is the preferred final solution. Note that the x-axis isthe residual calculated using the conditioned matrix A from step 7 and m is not scaled toabsolute value. The inset is a histogram of the Kagan angle standard deviation showingthe distribution of stable (green), likely (yellow) and poor (red) MTs.amplitude constraint misfit of 1.82. Figs B.4 to B.13 show the individual MT solutions along withtheir numbers of constraints and uncertainties.Of the 305 (stable and likely) events retained for analysis, 23 events had a local magnitude (cal-culation based on S-waves) reported by the Geological Survey of Canada. The moment magnitudesMW recovered from the MT inversion range between 1.07 to 3.47 and are shown in Fig. 3.7. Allstable and likely MTs are displayed in Fig. 3.10. Their diversity of orientations highlights the com-plexity of the fault network; although groups of similar mechanisms can be observed in localizedareas (e.g., Fig. 3.10 A and C; see also Fig. B.14 that displays MT solutions in map view clusteredby similarity). As shown in the Kaverian diagram (Fig. 3.11), most of the events are reverse faultsand reverse strike-slip faults, many with a northeasterly strike.3.5.2 Stress InversionThe individual MTs were clustered into 6 geographic groups using a K-means algorithm and per-form stress inversions using the Vavrycˇuk [2014] method. The 305 best (stable and likely) events,55Figure 3.10: Stable and likely MTs colour-coded by depth with zoom-in sections. All MTsare shown in lower hemisphere representation. Map region shown on Fig. 3.6.with their uncertainties, were employed in the inversion. Results are displayed in Fig 3.12.In their study of the regional (southwest B.C. and northwest Washington) stress regime basedon focal mechanism data, Balfour et al. [2011] showed that orientation of the maximum horizontalcompressive stress (σHmax) varies with distance to the trench: along the Vancouver Island coast,σHmax is approximately margin-normal due to stress accumulation along the locked subduction zoneinterface; whereas 100-150 km inland of the coast, σHmax is margin parallel due to the northward56Figure 3.11: Kaverian diagram illustrating representation of fault styles determined usingFMC software [A`lvarez Go´mez, 2019].Figure 3.12: Stress inversion results. Maximal horizontal compressive stress directions areplotted as colored arrows corresponding to the contributing event epicenters (coloreddots) displayed on map of epicenters. Note that the green events represent the shal-low group (< 12 km). The white arrows represent the average horizontal compressivestress using stable and likely MT solutions, and compare closely with that reportedby Balfour et al. [2011] plotted in black. Black epicenters represent earthquakes forwhich MTs were not retrieved. Purple lines are the 52 and 56 km oceanic Moho con-tours from Bloch et al. [2023]. The stereonet shows the direction of all stress axes(σ1:compression; σ2: neutral; σ3:tension) with the same colour code as map.push of the clockwise rotation of Oregon block [McCaffrey et al., 2000]. However, the local stressorientation recovered for the San Juan cluster [Balfour et al., 2011] deviates from the regionalpattern in a counter-clockwise direction, with σHmax oriented northwest-southeast. In their study,Balfour et al. [2011] used 38 earthquakes from the San Juan Island cluster with depths rangingfrom 0.03 to 24.5 km that occurred between August 1980 to October 2009. Their data set is almostindependent of that employed here (317 earthquakes), with only one shared event (2009/06/25).57Note that the close agreement (∼ 3◦) between their recovered stress orientation and the average (inwhite) solution for the entire MT suite.3.5.3 DiscussionThe San Juan earthquake cluster represents a dense concentration of seismic activity apparently as-sociated with a local perturbation to the regional stress regime [Balfour et al., 2011]. It comprisesspatiotemporally segregated subclusters that, in aggregate, outline a broadly north dipping fault net-work (Fig. 3.8 C). The temporal sequences associated with the individual subclusters are typicallynot preceded by large magnitude foreshocks, but rather comprise persistent low-magnitude events.Moreover, several of the temporally confined subclusters appear to form near-vertical structures (seeeast-west profile in Fig. 3.8 B). These observations are consistent with swarm behaviour, suggestingfluid migration as the underlying cause for seismicity [e.g., Vidale and Shearer, 2006]. The sourceof fluids is interpreted to be the underlying, dehydrating Juan de Fuca plate.A range of independent observations supports this interpretation. The San Juan cluster repre-sents one of the largest concentrations of forearc seismicity in the region, a part of Cascadia that hasbeen imaged in a number of geophysical studies. The tomographic models of Savard et al. [2018]and Merrill et al. [2022] locate the San Juan cluster in a zone of low Vp/Vs whereas magnetotelluricwork [Calvert et al., 2020; Rippe et al., 2013; Wannamaker et al., 2014] characterizes the crust hereas highly conductive. Both observations favour an interpretation involving the presence of fluids.A correlation between seismicity and low Vp/Vs has been observed throughout the Puget Sound -Strait of Juan de Fuca region [Merrill et al., 2022; Savard et al., 2018] and may manifest, at leastpartially, the presence of silica precipitated by fluids rising from the slab [Audet and Bu¨rgmann,2014; Bostock et al., 2019; Ramachandran et al., 2006]. Furthermore, the MOCHA magnetotel-luric experiment that imaged the forearc/arc between southern Washington and northern California[Egbert et al., 2022] has revealed conductivity anomalies interpreted to represent a fluid concentra-tion peak between 17.5 and 30 km depth controlled by crustal structure and metamorphic processes.This depth range matches that of the San Juan cluster (majority of the events between 13 and 23 kmdepth).Another interesting correspondence involves the location of the San Juan cluster landward of apronounced deep tectonic tremor concentration, as mapped by Armbruster et al. [2014] (see Fig.3.13), marking the downdip limit of tremor beneath the eastern Haro Strait. As such, it representsa local example of the more general tendency for forearc seismicity and tremor to be anticorrelatedin their respective epicentral (and depth) distributions in Cascadia [e.g., Boyarko and Brudzinski,2010; Kao et al., 2009]. Tremor is thought to occur where fluids are trapped at near-lithostaticpressures in the vicinity of the plate boundary [Audet et al., 2009; Rubinstein et al., 2008]. Bostocket al. [2019] have proposed that the complementary epicentral distributions of tremor and forearccrustal seismicity results from a landward breach of the permeability seal responsible for maintain-58ing high fluid pressures in the tremor zone, plausibly due to i) eclogitization of oceanic crust whichinvolves a ∼ 10% volume reduction, and/or ii) high strain associated with along-strike curvature inthe subducting slab together with increased slab dip, that allows entry of fluids into the overlyingcrust and the promotion of fluid-driven swarm seismicity like that proposed here below San JuanIsland. This interpretation is supported by the resistivity dip profiles across the Cascadia margincompiled by Wannamaker et al. [2014], indicating high conductivity zones emanating from 35-40km depth near the downdip termination of tremor and evolving landward, in some instances as here,into crustal depths. It is noted, however, that according to Fig. 3.7 there appears to be no obvioustemporal correlation between the occurrence of tremor and seismicity below San Juan Island. If theFigure 3.13: Seismicity Map. The red dots are the earthquake part of the San Juan cluster.Blue dots are tremors between 2003 and 2007 relocated by Armbruster et al. [2014].Gray dots are other earthquakes detected in the area. The orange box identifies the areaused for Fig. 3.7. The oceanic Moho from Bloch et al. [2023] is displayed in purplelines. The bold black lines are the Devil’s Mountain Fault Zone (DMFZ) [Barrie andGreene, 2018], Skipjack Island Fault Zone (SJIFZ), the South Lopez–San Juan Fault(SLSJF) and the Leech River Fault (LRF). Solid black lines are other mapped faults inthe region Greene and Barrie [2022]; Greene et al. [2018].general north-dipping geometry of the seismicity clusters evident between∼ 13-23 km depth in Fig.3.8 C is extrapolated updip, it would surface in proximity to the trace of the Devils Mountain FaultZone (DMFZ) or South Lopez–San Juan (SL-SJF). The DMFZ and its westward extension into theLeech River Fault (LRF) on southern Vancouver Island are believed to be currently active [Johnsonet al., 2001; Morell et al., 2017], and potentially capable of generating M6-7 earthquakes. At thesurface, the DFMZ is reported to dip at between 45◦ to 75◦ to the north [Johnson et al., 2001]. It is59therefore conceivable that the deeper reaches of seismicity present below San Juan Island (Fig. 3.8)are associated with this structure. However, the more limited, shallower (< 13 km depth) events donot follow this extrapolation but instead occur directly above the deeper seismicity suite and mayindicate intersection of the fault network with steeper secondary structures extending to the surface.Note that Lidar and bathymetry data [Greene and Barrie, 2022] reveal that the San Juan IslandArchipelago, including San Juan Island in particular, are crisscrossed by multiple fault strands atthe surface. Helium isotope ratios measured at hotsprings to the south in Puget Sound, where theSeattle Fault crosses Bainbridge Island, provide evidence that slab-derived fluids do reach the sur-face along major crustal faults [McCrory et al., 2016]. Additionally, on Vancouver Island near thetransition from the DMFZ to the LRF, mineralization is interpreted as resulting from the distillationof volatiles and solutes from the Juan de Fuca plate [Fyfe et al., 1987]. Moreover, Greene and Bar-rie [2022] have identified pockmarks and a mud volcano in the proximity of the DMFZ, providingevidence of fluid seepage in the region (see their Fig. 10). However, Johnson et al. [2022] haveanalyzed methane emissions in gas plumes issuing from major fault traces on the southern PugetSound seafloor and do not record an obvious thermal or chemical signature that would favour adeep-seated slab origin.3.6 ConclusionsThe relative moment tensor inversion approach of Plourde and Bostock [2019] has been modifiedto incorporate constraints from P-wave first motion polarities, thereby resolving the moment-tensorpolarity ambiguity that exists in the absence of available reference moment tensors. This adapta-tion is particularly important for studying small earthquakes, where conventional moment tensoralgorithms based on waveform modelling may fail due to the complexities of high-frequency wavepropagation. Moreover, the challenge of moment tensor scaling or, equivalently, scalar momentrecovery, is addressed through a post-inversion step utilizing available local magnitude information.The efficacy of the algorithm is demonstrated through both synthetic examples and application toa dense concentration of seismicity beneath San Juan Island in Washington State, in the Cascadiaforearc.The principal stress directions obtained from the moment tensors recovered for the San JuanIsland dataset exhibit little variability across subclusters and are consistent with findings from aprevious, regional focal-mechanism study conducted by Balfour et al. [2011], using an independentdataset. This agreement provides increased confidence in the observation that the stress regime localto the San Juan Island region is significantly perturbed from the regional trends defined in Balfouret al. [2011]. Moreover, other characteristics of the seismicity, namely the temporal dependenceof moment magnitudes, hypocentral distribution, and association with low Poisson’s ratio and highconductivity imaged in previous structural studies are consistent with an interpretation involvingseismic swarms. More specifically, the San Juan seismicity is attributed to the migration of fluids60originating from the dehydration of the underlying subducting Juan de Fuca plate. These findingscollectively enhance the understanding of earthquake processes in the Cascadia forearc.61Chapter 4Moment Tensor Analysis in the SouthernHaida Gwaii: Insights intoUnderthrusting Processes4.1 Chapter SummaryOblique subduction zones often result in slip partitioning, where the thrust component of slip (nor-mal to the trench) is restricted to the subduction interface, and the strike-slip component (parallel tothe trench) occurs on a near-vertical strike-slip fault that defines the landward boundary of a crustalsliver. This scenario is observed in southern Haida Gwaii, an offshore island archipelago off thewest coast of British Columbia, Canada. In this study, several earthquake clusters in the vicinityof the crustal sliver are analyzed to gain a better understanding of slip partitioning. The region’ssparse seismic network poses challenges for studying microseismicity. To address this shortcoming,a novel relative moment tensor inversion technique is employed to recover mechanisms from oversixty earthquakes. Our results indicate that deeper events exhibit predominantly reverse faultingmechanisms, consistent with the underthrusting of the Pacific Plate beneath the North AmericanPlate.4.2 IntroductionDynamic interactions at plate boundaries offer valuable insights into Earth’s geological processes.A particularly interesting example is afforded by the case of oblique subduction, where one platesubducts beneath another at an oblique angle, resulting in slip partitioning into strike-slip (trenchparallel ) and thrust (trench normal) components [Jarrard, 1986]. In most cases, a forearc sliver isformed between the trench/deformation front and a strike-slip fault as illustrated in Fig. 4.1. Exam-62Figure 4.1: Schamatic showing the tectonic setting of oblique convergence adapted for SouthHaida Gwaii with the expected moment tensors shown by the beachballs. Modified fromCassidy et al. [2014].ples include the Chile trench and Liquine-Ofqui fault [Stanton-Yonge et al., 2016], Indo-Australianplate subduction along Indonesia and the Great Sumatran fault [Ishii et al., 2013], the PuysegurTrench and Puysegur Fault [Shuck et al., 2021], Philippine sea plate subduction and the PhilippineFault [Barrier et al., 1991], and Pacific plate convergence at Haida Gwaii and the Queen Charlottefault [Cassidy et al., 2014; Oliva et al., 2024].This paper focuses on Haida Gwaii, a high-obliquity endmember, where the Pacific plate movesat between 48-53 mm/yr with respect to the North American plate [Brothers et al., 2020; DeMetsand Merkouriev, 2016] at an angle between 6-25◦ clockwise from the offshore deformation front[Brothers et al., 2020; DeMets and Merkouriev, 2016; DeMets et al., 2010; Tre´hu et al., 2015] (seeFig. 4.2 for geographic context). This is among the youngest examples of incipient subduction onEarth having initiated between 12-11 to 6-5 Ma [Austermann et al., 2011; DeMets and Merkouriev,2016]. Due to its young age and high obliquity, no volcanic arc and no clear Wadati–Benioff zoneare observed in the region. However, several geologic observations support the occurrence of under-thrusting [Hyndman and Hamilton, 1993]: (1) the presence of an accretionary sedimentary wedgenamed the Queen Charlotte Terrace (QCT), (2) an off-shore depression and bathymetric high 100km seaward of the depression known as the Oshawa rise, (3) characteristic parallel bands of low-high gravity anomaly, (4) the heat flow pattern [Hyndman et al., 1982], (5) thrust focal mechanisms63in the region [Bird, 1997; Be´rube´ et al., 1989; Kao et al., 2012; Ristau et al., 2007], and (6) receiverfunctions that indicate the presence of a landward dipping low-velocity zone [Bustin et al., 2007;Gosselin et al., 2015; Smith et al., 2003]. Furthermore, receiver functions from [Bustin et al., 2007;Gosselin et al., 2015] and gravity models [Hyndman, 2015] are consistent with underthrusting to adepth of ∼ 40 km and extending 65-70 km inland from the Queen Charlotte Fault (QCF).The QCF is a∼900 km long, right-lateral strike-slip fault that merges with the Fairweather Faultin mainland Alaska to the north. Near the southern termination of the QCF (i.e. south of latitude∼53.4◦), it overlaps with the Haida Gwaii Thrust (HGT), forming the QCT that together with un-derlying basement, plays the role of forearc sliver [Riedel et al., 2020]. The sliver is approximately25 km wide and reaches a depth up to ∼20 km based on seismicity distribution from Oliva et al.[2024].The Haida Gwaii region has hosted a majority of the largest Canadian earthquakes recordedwith modern instrumentation. The 1929 M7.0 earthquake (no mechanism recovered, but strike-slipmotion suspected), the 1949 M8.1 (largest instrumentally recorded earthquake in Canada), and theM7.4 1970 event, all occurred along the transform Queen Charlotte fault (QCF) (see Fig. 4.2 forlocation). In contrast, the M7.8 2012 earthquake occurred on the HGT and exhibited reverse slip.The occurrence of major strike-slip and thrust earthquakes in the same region is an expression ofslip partitioning.A new regional earthquake catalogue [Oliva et al., 2024] further characterizes the nature offaulting in the region through the recovery of more than 11,000 well-located hypocenters of mi-croearthquakes that occurred largely as aftershocks to the 2012 earthquake. However, the paucity ofavailable moment tensors for the region renders the driving stresses of these events uncertain. Seis-mic stations in the region are sparsely distributed. Only events with a local magnitude (ML) greaterthan 3.5 are routinely processed for Centroid Moment Tensors (CMTs), and no event with ML < 3.9is represented within the regional CMT catalogue [Kao et al., 2012]. Approximately twenty focalmechanisms for events between 1984 and 1995 with ML between 2.0 and 4.3 along the Queen Char-lotte transform boundary were recovered by Bird [1997] using first motion P-wave polarities andup to 9 stations located on Haida Gwaii and the adjacent Mainland. Her mechanisms indicate boththrust and strike-slip faulting. A microseismic study by Be´rube´ et al. [1989] also documented thrustmechanisms in the vicinity of the HGT during a 1983 ocean-bottom seismometer deployment.The goal of this study is to recover moment tensors near the expected sliver location to gain adeeper understanding of how slip-partitioning is accommodated in the region. Notably, the earth-quake locations near the Queen Charlotte Fault (QCF) form a distinct linear structure that deviatesfrom the seafloor trace of the QCF. This discrepancy raises important questions about the role of theQCF in the seismotectonic framework of the region. Despite the potential utility of focal mecha-nisms and moment tensors from microseismic (MW ≤∼ 3) events for unravelling the distribution ofstress across the complex plate boundary at Haida Gwaii, the generally sparse distribution of per-64Figure 4.2: Overview of the study region. The Haida Gwaii region is an archipelago with∼ 150 islands. Previously reported moment tensors (lower hemisphere) with Mw ≥ 6.0are shown [Kao et al., 2012; Ristau et al., 2007; Rogers, 1983; U.S. Geological Survey,2017c]. Black dots are ∼ 11,000 well-located earthquakes from Oliva et al. 2024. Redtriangles are stations used in this study. The three coloured squares are the center loca-tions of the earthquake groups analyzed here. The surface traces of the Queen CharlotteFault (QCF), the Revere-Dellwood Fault (RDF) and the Haida Gwaii Thrust (HGT) areshown as dark brown lines [Brothers et al., 2020; Oliva et al., 2024]. The bathymetryand topography are from GEBCO [2023].65manent seismometers located primarily landward of the seismicity poses major challenges. Here,the recovery of over sixty moment tensors (MTs) for microseismic events in the region of the HGT,using a newly developed relative moment tensor inversion method [Plourde and Bostock, 2019],described in more detail in Chapter 3 [Drolet et al., 2024] is reported and discussed. This approachincorporates relative amplitude measurements in addition to P-wave polarities to mitigate the effectsof the sparse and irregular station distributions.4.3 Data and MethodologyThe approach used for recovering moment tensors from microseismic events requires that theseevents occur in dense clusters such that relative amplitude measurements of multiple events madeat an individual station are independent of propagation effects and reflect dominantly differences inmoment tensor strength and radiation pattern [[Plourde and Bostock, 2019], Chapter 3 of this thesis].Eleven separate clusters (henceforth referred to as “groups”) of seismicity from the well-located11,000 event catalogue of Oliva et al. [2024] were originally targeted. These groups were selectedbecause of their high earthquake density and geographic location. As shown in Fig. 4.3, the majorityof these groups lie along an offshore trend of seismicity that approximately follows, albeit at a smallangle to, the QCF trace with deeper events occurring close to the inferred intersection with the HGT.A few groups lie landward of this trend, including one intriguing inland group L1, analyzed in theAppendix C.3. All selected events occurred after the M7.8 October 28, 2012 earthquake coincidingwith an increase in seismic station coverage effected to improve aftershock characterization. Eventsrecorded during and prior to the 2012 earthquake were not recorded by sufficient numbers of stationsto permit analysis. Moreover, for reasons explained below, eight of the selected groups failed toproduce reliable MT solutions and so will not be analyzed. Table 4.1 provides details on the finalselection of event groups with those clusters along the offshore trend labelled as O1-O3.Group Date RangeGood / Recovered /Total EventsDepth Range(km)Area Estimate(km2)O12013-01-022020-01-2123 / 33 / 36 12 - 18 30O22013-01-032016-01-1614 / 34 / 35 15 - 21 20O32013-01-032015-01-2927 / 45 / 47 14 - 19 15Table 4.1: General information about each earthquake group. Good/Recovered/Total Eventsrefers to the number of MTs selected for further analysis, the number of recovered MTs,and the total number of events that are part of the inversion, respectively.Notwithstanding the deployment of additional stations following the main 2012 event, seismic66Figure 4.3: Overview of the earthquake groups analyzed in this study. Same data as figure 4.2.The square outlines delimit the different earthquake groups characterized by respectivecolour codes. Grey groups were investigated, but no MTs were recovered with a highdegree of confidence. Group O1 includes four reference events, as shown in the bottomleft corner. Those MT solutions from Kao et al. [2012] have been converted to theirdouble-couple equivalents , and the black lines identify their locations as employed forthe CMT inversion. Note the Queen Charlotte Fault (QCF) splits into two segments nearlatitude of 52.25◦ [Brothers et al., 2020]. HGT: Haida Gwaii Thrust faults. Brown linesare surface faults [Cui et al., 2017; Thompson, 1990]. White lines are the cross-sectiondirection used in Fig. 4.5.The black line 1-1’ and 2-2’ are the direction of the cross-section on Fig. 4.4.67station coverage remained relatively sparse, with poor azimuthal coverage, in particular for off-shore earthquakes. These conditions present significant challenges for focal mechanism recovery inmicroseismic studies. To mitigate these deficiencies, a new MT inversion method based on relativeamplitudes (see Chapter 3) is employed. This approach simultaneously recovers moment tensors forclustered events by utilizing constraints from the relative amplitudes of P- and S-waves recorded atcommon stations, P-wave polarities, and local magnitude estimates. Relative amplitudes constrainthe radiation pattern, while a known (or “reference”) MT, if available, and/or P-wave first motionpolarities determine the sign of the MT. In addition, either a reference MT or local magnitude infor-mation is used to establish seismic moments for all events.Waveforms were downloaded from the Incorporated Research Institutions for Seismology (IRIS)and Natural Resources Canada (NRCan) data centers. Phase picks were obtained from Oliva et al.[2024] comprising a mix of analyst and automatic phase picks, all of which were visually checkedand revised as needed. P-wave first motion polarities were determined using the principal compo-nent analysis (PCA) approach developed by Drolet et al. [2024] (see section 3.3.3 of Chapter 3),which yields more accurate and consistent polarities. Cross-correlation delay times were used toalign the waveforms, allowing accurate relative amplitude measurements based on principal com-ponents analysis [Plourde and Bostock, 2019]. Local magnitudes were obtained from the NRCancatalogue [Natural Resources Canada, 2021a] for all groups. Furthermore, only group O1 possessesknown reference moment tensors (4 total), as supplied by the Geological Survey of Canada (GSC)[Kao et al., 2012] and displayed in Fig. 4.3. Note that centroid locations determined by GSC, whoemploy 3 regional 1D velocity models, differ from those determined by Oliva et al. [2024], whocomputed hypocenters using a simultaneous inversion for earthquake locations and a regional 3Dvelocity model. Oliva et al. [2024] hypocenter locations were used for this analysis.Like first-motion P-polarity focal-mechanism determination, relative MT inversion is highlysensitive to the assumed velocity model since the latter controls the ray azimuth and take-off anglesand, so, the point on the focal sphere to which P-polarities and relative amplitudes are mapped.Ideally, the full 3D velocity model of the region from Oliva et al. [2024] would be used to extractthese angles. However, this approach is computationally expensive and complicated by the presenceof several low-velocity zones. Therefore, the inversions were run using five different suites of 1Dvelocity models: (1) 1D velocity models extracted from Oliva et al. [2024] for each event-station-phase combination as the average velocity between the event and station; (2) 1D velocity modelsdetermined for each earthquake group by averaging velocities within a 50 km radius circle centredon the group, based on Oliva et al. [2024]’s 3D velocity model (for the offshore groups, only thelandward half of the circle was used in the averaging); (3) the standard 1D velocity model usedto locate earthquakes in the Haida Gwaii region by the Geological Survey of Canada [Farahbodand Cassidy., 2022]; (4) a 1D velocity model employed in receiver function work performed at thestation MOBC [Bustin et al., 2007]; and (5) three station-dependent 1D velocity models used by68the Geological Survey of Canada for centroid moment tensor inversion in the region [Kao et al.,2012]. From now on, velocity model 1 will be referred to as 1D Oliva et al. [2024] phase averagemodel and model 2 as 1D Oliva et al. [2024] group average model. Since all velocity models are1D, azimuth is defined by the event-to-station direction. Calculation of the take-off and incidenceangles were performed using the python package obspy.taup [Crotwell et al., 1999]. Low-velocityzones were removed as needed to ensure stable ray paths to stations.For MT inversion, only events with phases recorded at a minimum of four seismic stationswith an azimuth gap smaller than 300° were considered. Only clearly defined phases with realisticarrival times and waveforms with high signal-to-noise ratio were selected. Due to the small numberof stations, every moment tensor was forced to be double-couple, thereby effectively reducing thenumber of unknown parameters from 6 to 4 per event. For each group, the inversion results fromthe five different velocity models were compared, and the velocity model yielding the lowest misfit(highest polarity agreement and lowest relative amplitude residual) was employed to define the finalsolution. Only events with a median Kagan angle (the minimum angle required to rotate one DCMT into another, [Kagan, 1991]) deviation smaller than 30° were selected. That is, for each event,these deviations were computed between its best solution (solution with the lowest misfit) and thesolutions for the other velocity models. Eight of the groups did not pass this last quality control; seeFig 4.3 for their location. Note that these groups tended to comprise shallow (< 12 km) events likelywithin the low-velocity QCT, suggesting difficulties in achieving stable ray paths in this region ofstrong velocity contrasts. Note that MTs for the inland lineation L1 are presented in the appendix,as the results have high uncertainty, although they appear to be reasonable.As a final processing step, a stress inversion for each group is performed using the Vavrycˇuk[2014] method with assigned errors to focal mechanisms of 20°. This technique is an iterative jointinversion for stress and fault orientation. For each iteration, the principal stress orientation andstress ratio (based on the relative magnitudes of the principal axes) are found using the Michael[1984] method. Then, the ambiguity of identifying the fault orientations is resolved by applyingfault instability constraints using the stress ratio result. The most unstable fault planes are then usedas rupture planes in the following inversion. Five hundred noisy realizations were performed toestimate the principal stress uncertainty.4.4 Results4.4.1 Moment Tensor SolutionsTable 4.2 presents the chosen velocity model for the final solutions (see Appendix C.2.1 Table C.1for detailed misfit values). For all three remaining groups, the Bustin et al. [2007] model yieldsthe lowest misfit. One might have expected that the velocity model based on the average velocitybetween the event and station from a 3D velocity model (model 1) would yield the best results.69However, the large velocity contrasts occurring between earthquakes located offshore and stationson land likely render averaged velocity models less than ideal for characterizing ray paths. Note thatgroup O1 displays the highest relative amplitude residuals group, which is explained by the largermagnitude events that form this group.Group ScalingBest VelocityModelPolarityAgreementRelative AmplitudeResidual (Nm)O1 4 MW + 4 ML Bustin et al. (2007) 95/122 (78%) 9.61e+30O2 1 ML Bustin et al. (2007) 81/95 (85%) 1.01e+23O3 14 ML Bustin et al. (2007) 106/128 (83%) 2.72e+24Table 4.2: Properties of the final solutions for each group. The polarity agreement and relativeamplitude residual represent all events included in the inversion.Here, only the best-constrained events (those with a median Kagan angle under 30°) from therelative MT solution with the lowest misfit are presented. However, Figs C.2 to C.5 in AppendixC.2.2 display the MT solutions for each velocity model and all inverted events in each group. Fig.4.5 shows the best recovered double-couple MT solutions. The estimated top of the underthrustPacific Plate (or HGT) based on a constant angle of 27° (for group O1) and 30° (for group O2 andO3) from the deformation front is also demarcated in Fig 4.5 A. For Fig 4.5 C and E, the simplifiedmodelled slab-top occurs 1-2 km above the upper right corner and thus is not visible in the figure.Overall, Figs 4.5 displays many reverse faulting mechanisms oriented in the expected dip direc-tion of the HGT. Fig. 4.6 A displays a Kaverian diagram for all off-shore earthquakes and revealsthat most recovered earthquakes have reverse mechanisms.4.4.2 Stress Inversion ResultsFig. 4.7 displays the stress inversion results in terms of principal stress orientations (σ1 (compres-sion)≥ σ2 ≥ σ3 (tension)), for all groups along with their associated error margins determined fromthe noisy realizations. Note that group O2 reveals the most variability in stress orientation due tothe low number of events in the inversion (14 events). According to Martı´nez-Garzo´n et al. [2016],for focal mechanisms with noise levels up to 20°, a minimum of 20 events is necessary to retrievethe stress field orientation with an accuracy of 10°.Including the error margins, all groups display σ3 axis close to vertical and so are consistent withdominantly reverse faulting, albeit with significant σ2 girdling for O2. Also, note that σ3 generallypoints to the east, the dip direction of the subducting slab. These results are consistent with Fig. 4.6that summarizes dominantly reverse faulting.Another useful means of characterizing stress regimes is the stress ratio, R [Etchecopar et al.,1981], also used by Vavrycˇuk [2014] to identify the fault planes in stress inversion. The stressratio (also called relative stress magnitude, eigenvalue ratio or shape ratio) quantifies the relative70Figure 4.4: Cross-section showing the cluster location compared to the possible position ofthe sliver as defined by an assumed subduction interface (blue dotted line) and the QCFsurface traces(s) (brown dotted lines). B displays two possible locations of the QCF asthe cross-section is near a bifurcation of the QCF surface trace. Dots are earthquakeswithin ± 10 km on A and ± 15 km on B of the cross-sections. A) 1-1’ segment from 4.3showing group O1 in pink. B) 2-2’ segment from 4.3 showing group O2 in yellow andO3 in orange.71Figure 4.5: Moment tensor solutions colour-coded by depth for groups O1 (A,B), O2 (C,D)and O3 (E, F). Plots A, C and E are map views (lower hemisphere MTs) with pink linesrevealing the orientation of SW-NE cross-sections in plots B, D and F . Note that B, Dand F show the MTs in side view as indicated by the blue arrows. The brown line on Arepresents the surface trace of the Queen Charlotte Fault. The blue line on B, D and F isan estimate of the top of Pacific Plate (see main text for description).72Figure 4.6: Kaverian diagram, showing the different types of earthquake mechanisms forgroups O1 (pink rim), O2 (yellow rim) and O3 (orange rim). This diagram was pre-pared using Focal Mechanisms Classification (FMC) python code from A`lvarez Go´mez[2019].Figure 4.7: Stress inversion results for the three earthquake groups, with same colour codeas Fig. 4.3. Gray lines on the stereonet are plotted in 15° increments. Diamonds: σ1(compression); circles: σ2 (neutral); squares: σ3 (tension). The left stereonet showsthe preferred solutions for all groups, and on the right are the solutions for 500 noisyrealizations. The corresponding stress ratio value, R, is included under each individualgroup stereonet.magnitudes of the principal stress axes and is defined as:R =σ1−σ2σ1−σ3 , (4.1)where R can take values between 0 (uniaxial tension) and 1 (uniaxial compression). The stress ratiosare at the high end of values characteristic of compressional regimes.734.5 Discussion4.5.1 Main ResultsUnfortunately, I was unable to recover reliable MT solutions for offshore clusters of seismicity atshallow depths where seismicity is thought to occur within the QCT (see Oliva et al. [2024], theirFig 8). However, I was able to retrieve 64 MTs for groups O1-3 at greater depths in the regionwhere the QCF is potentially intersecting the underthrusting Pacific Plate. My results reveal thatthe majority of microseismic events at these depths represent reverse faulting. Group O1 presentsmore variety in MTs that could be explained by its location near the intersection between the thrustinterface and the respected QCF. The results indicate that the governing stress regime in this regiondelivers dip-parallel, partitioned slip at the frictional contact between the Pacific plate and QCT. Inother words, the least compressional axis (σ3) is pointing steeply downdip. The precise location ofthe top of the Pacific plate here is not well constrained by the seismicity. However, it is anticipatedthat the thrust events occur within a possibly distributed plate boundary zone.4.5.2 Slip PartitioningThe recovered mechanisms are predominantly pure thrust pointing downdip along the slab, sug-gesting they occur on the thrust interface of the sliver, in line with the traditional slip-partitioninghypothesis (refer to Fig. 4.1). However, the three studied groups are located downdip of the inferredsliver landward boundary (i.e. the QCF), as shown in cross-sections 1-1’ and 2-2’ in Fig. 4.4.Since the general seismicity pattern does not align with a near-vertical Queen Charlotte Fault(QCF), the simplest interpretation is that the subsurface geometry of the QCF exhibits curvature orinclination. This adjustment would position groups O1, O2, O3 near the intersection of the sliverlower and landward boundaries. In this case, these earthquakes, as aftershocks, reflect the releaseof residual compressive stress following the 2012 M7.8 event. They are located near the downdiprupture terminus of the 2012 earthquake, estimated to extend to a depth of 19 km by Hobbs et al.[2015]. Additionally, the overlying upper terrace earthquakes are then attributed to diffuse stressrelease in the overriding crust.4.5.3 Methodology limitationsAs shown in Fig. 4.3, about half of the seismicity clusters targetted for our investigation did notpass our final quality control criterion. The scarcity of stations makes this a particularly challengingregion for studying microseismicity. Most rejected groups were shallow, and largely concentratedin the sediments of the QCT, rendering it more difficult to accurately trace polarities and relativeamplitude measurements to the focal sphere. When using a local seismic networks, deeper eventstend to have straighter ray trajectories than shallower ones, resulting in better-constrained take-off74angles. In addition, large and poorly modelled lateral variations in shallow velocities will furtherexacerbate the accuracy of traced rays.It is important to note that the traditional P-wave first motion technique would not have beeneffective in this scenario even for groups O1−3, because most events comprised only 0-4 reliablepolarity measurements. Using the relative P- and S-wave amplitudes significantly increased thenumber of constraints, allowing the recovery of reliable MTs.4.6 ConclusionIn this study, I was able to recover reliable MTs from only 3 of the 11 seismicity clusters that wereoriginally targeted to investigate the stress regime associated with slip partitioning at the highlyoblique convergent margin of Haida Gwaii. These 3 clusters occur in the deeper reaches (≥ 12 km)of the QCT forearc sliver and likely reflect partitioned dip slip within or near a possibly distributedplate boundary zone. As aftershocks to the Mw 7.8 Haida Gwaii thrust event, they likely representslip due to release of residual stresses. The remaining clusters occurred at shallower levels withina region of strong lateral velocity gradients. I believe that this factor, along with the sparse stationcoverage restricted to eastern azimuths, resulted in a failure to generate reproducible moment tensorsolutions for different velocity models.75Chapter 5ConclusionThis thesis aimed to explore two main research questions related to the insights that small earth-quakes can provide into regional stress regimes and methods by which these events can be moreeffectively analyzed. Through a series of detailed studies, a significant understanding of both topicswas gained. In the following paragraphs, I review how each question was addressed and the keyfindings of this thesis.5.1 Question 1: What information can we extract from smallearthquakes?In Chapter 2, I conducted a detailed analysis of the aftershock sequence of two significant intraslabearthquakes in Alaska, revealing key insights into the rupture processes and fault dynamics asso-ciated with these earthquakes. These aftershock sequences comprised hundreds of well-recordedsmall events, rendering them an ideal data set for the task. For the Anchorage sequence, the af-tershock analysis and moment tensor recovery identified two distinct fault zones: a west-dipping(northern) and an east-dipping (southern) group, both of which correspond to the mainshock mo-ment tensor nodal planes. This observation suggests the possibility of simultaneous rupture onboth nodal planes. It is plausible that the mainshock initiated on the east-dipping plane and subse-quently triggered rupture on the pre-existing west-dipping fault, located 3 km updip of the main-shock hypocenter, as suggested by Guo et al. [2020]. This scenario could explain the misalignmentbetween the seismicity and the mainshock hypocenter. The occurrence of a foreshock on the west-dipping fault indicates that this fault already existed prior to the mainshock. Lastly, the diversity inmechanisms recoverd in this study highlights the complexity of the fault structure activated by themajor, damaging intraplate event.For the Iniskin sequence, the aftershock distribution delineated the extent of the mainshock rup-ture. The observations suggest that the Iniskin earthquake reactivated a pre-existing outer rise fault,with updip rupture initiating within the oceanic mantle and terminating at a stratigraphic horizon76within the slab, likely at the slab Moho or plate interface. The inversion of moment tensors fortime-dependent stress revealed a stable stress orientation characterized by downdip tension consis-tent with the downward pull of the negatively buoyant subducting plate. However, the stress ratiounderwent a sharp decrease immediately following the mainshock, followed by a gradual recovery.Assuming anisotropic fracturing within the Pacific plate, this stress evolution is interpreted as evi-dence of pore fluid and pore pressure diffusion away from the principal fault zone, likely driven bymetamorphic dehydration reactions. This interpretation supports the notion that dehydration-drivenstress transfer played a crucial role in triggering the earthquake.In Chapter 3, the analysis of small earthquakes within the San Juan cluster, an anomalouslydense group of earthquakes located in the Cascadia forearc, revealed multiple localized bursts ofevents with comparable magnitudes. Moment tensors displayed systematic variability as a functionof geographic location. These observations characterize earthquake swarms. The region displaysa relatively homogeneous stress regime, distinct from the larger regional-scale pattern previouslydescribed by Balfour et al. [2011]. Additionally, reverse faulting emerged as the predominantmechanism within the cluster. Since earthquake swarms are believed to be driven by migratingfluids, I infer that the source of these fluids is metamorphic dehydration reactions taking place in theunderlying Juan de Fuca plate and consistent with other geophysical observations from the region.In Chapter 4, the study of microseismicity in the Haida Gwaii region has provided additionalinsights into slip partitioning in oblique subduction zones. The analysis of small earthquakes deep(> 12 km) within the Queen Charlotte Terrance sliver allowed us to identify reverse faulting mech-anisms at depth likely associated with the plate boundary zone formed by the overriding NorthAmerican plate and the subducting Pacific slab. The orientation of these mechanisms is consis-tent with slip perpendicular to strike occurring on the lower boundary of the silver, questioning thesubsurface geometry of the Queen Charlotte Fault.In summary, small earthquakes, when thoroughly analyzed, can yield rich datasets that providecritical insights into fault networks, rupture processes, stress regimes, and the role of fluids in seis-mic activity. Their frequent occurrence, compared to large earthquakes, allows for a more detailedunderstanding of the spatio-temporal evolution of these processes.5.2 Question 2: How can we improve or adapt pre-existing methods torecover mechanisms for microseismic events given theirchallenges?In Chapter 2, the relative moment tensor inversion method developed by Plourde and Bostock [2019]is utilized which proved effective for the two aftershock sequences studied, largely because the se-quences included large aftershocks with reliable moment tensors that contributed to the inversionprocess. However, in many studies of microseismicity, a known moment tensor is not available. To77address this challenge, Chapter 3 introduced additional constraints to the inversion process, includ-ing first-motion P-wave polarities and local magnitude estimates. I also improved the computationof double-couple moment tensor solutions from Plourde and Bostock [2019] by incorporating thestar norm, which increased the inversion speed. Furthermore, I parallelized the main computationalbottleneck of the code (step 8 - System Solution) to improve its efficiency. Through synthetic ex-amples and the application to the San Juan cluster, I demonstrated the robustness of my method.In Chapter 4, I applied this refined method to earthquakes in southern Haida Gwaii, a regionparticularly challenging for moment tensor recovery. During this study, I encountered some limita-tions due to the sparsity of the local seismic network. In particular, for shallow crustal earthquakesand large, poorly modelled velocity variations, it can be difficult to recover reliable moment tensorresults when using only a small number of stations located to one side of the sources.By incorporating additional constraints and optimizing computational efficiency, I have madesignificant progress in improving moment tensor recovery, even in challenging environments. How-ever, the limitations encountered highlight the ongoing need for further refinement, particularly inregions with sparse seismic networks or high levels of noise.5.3 Future WorkWhile this thesis has significantly advanced our understanding of small earthquakes and the meth-ods used to analyze them, several promising avenues for future research remain. Below are somedirections that would be particularly interesting to explore:1. Fig A.7 in Appendix A for Chapter 2 highlights that most intermediate-depth events in south-ern Alaska have a similar mechanism relative to the Pacific plate curvature, likely caused bythe reactivation of outer-rise faulting. Investigating the aftershock sequence of other largeevents in this region would reveal whether similar time dependency in stress ratio signaturesis present. Such a study would deepen our understanding of intermediate-depth earthquakesby supporting (or placing into question) the dehydration-driven stress transfer mechanism asdiscussed in this thesis.2. The relative moment tensor inversion method presented in Chapter 3 holds potential for fur-ther optimization and improvements. Some code components have not yet been fully op-timized, but once done it could significantly speed up the inversions process. It would beinteresting to see this method further automated and applied in an industrial setting, whereparameters such as rock velocity are well-known, and the seismic station coverage is dense.This could provide more insights into the method’s accuracy and limitations.3. If additional seismic stations are deployed in southern Haida Gwaii, or if an improved ve-locity model of the region becomes available, it would be highly beneficial to examine the78seismicity clusters at shallower levels. In particular, it would be useful to ascertain whethershallow seismicity occurring near the trace of the Queen Charlotte Fault is dominated byright-lateral strike-slip faulting as anticipated for the slip-partitioning model and supportedby our observations of dip-slip mechanisms for deeper events.In conclusion, this thesis has significantly improved our understanding of the insights that smallearthquakes can offer in elucidating fault dynamics and the role of fluids in seismic processes.Moreover, it has contributed to the refinement of analytical methods for microseismic events thatcan be utilized in future seismological research.79BibliographyAbers, G. A., J. Nakajima, P. E. van Keken, S. Kita, and B. R. Hacker, Thermal–PetrologicalControls on the Location of Earthquakes within Subducting Plates, Earth and Planetary ScienceLetters, 369-370, 178–187, doi:https://doi.org/10.1016/j.epsl.2013.03.022, 2013. → pages11, 33Albuquerque Seismological Laboratory (ASL)/USGS, United States National Seismic Network,doi:10.7914/SN/US, 1990. → page 124A`lvarez Go´mez, J. A., FMC—Earthquake Focal Mechanisms Data Management, Cluster andClassification, SoftwareX, 9, 299–307, doi:https://doi.org/10.1016/j.softx.2019.03.008, 2019. →pages 57, 73, 136Armbruster, J. G., W.-Y. Kim, and A. M. Rubin, Accurate Tremor Locations from Coherent S andP waves, Journal of Geophysical Research: Solid Earth, 119(6), 5000–5013,doi:https://doi.org/10.1002/2014JB011133, 2014. → pages 52, 58, 59Audet, P., and R. Bu¨rgmann, Possible Control of Subduction Zone Slow-Earthquake Periodicity bySilica Enrichment, Nature, 510, doi:https://doi.org/10.1038/nature13391, 2014. → page 58Audet, P., M. G. Bostock, N. I. Christensen, and S. M. Peacock, Seismic Evidence forOverpressured Subducted Oceanic Crust and Megathrust Fault Sealing, Nature, 457,doi:https://doi.org/10.1038/nature07650, 2009. → page 58Austermann, J., Z. Ben-Avraham, P. Bird, O. Heidbach, G. Schubert, and J. M. Stock, Quantifyingthe Forces Needed for the Rapid Change of Pacific Plate Motion at 6Ma, Earth and PlanetaryScience Letters, 307(3), 289–297, doi:https://doi.org/10.1016/j.epsl.2011.04.043, 2011. → page63Balfour, N. J., J. F. Cassidy, S. E. Dosso, and S. Mazzotti, Mapping Crustal Stress and Strain inSouthwest British Columbia, Journal of Geophysical Research: Solid Earth, 116(B3),doi:https://doi.org/10.1029/2010JB008003, 2011. → pages 56, 57, 58, 60, 77Balfour, N. J., J. F. Cassidy, and S. E. Dosso, Identifying Active Structures UsingDouble-Difference Earthquake Relocations in Southwest British Columbia and the San JuanIslands, Washington, Bulletin of the Seismological Society of America, 102(2), 639–649,doi:10.1785/0120110056, 2012. → page 5380Barrie, J. V., and H. G. Greene, The Devils Mountain Fault Zone: An Active Cascadia Upper PlateZone of Deformation, Pacific Northwest of North America, Sedimentary Geology, 364,228–241, doi:https://doi.org/10.1016/j.sedgeo.2017.12.018, 2018. → page 59Barrier, E., P. Huchon, and M. Aurelio, Philippine Fault: A Key for Philippine Kinematics,Geology, 19(1), 32–35, doi:10.1130/0091-7613(1991)019⟨0032:PFAKFP⟩2.3.CO;2, 1991. →page 63Bird, A. L., Earthquakes in the Queen Charlotte Islands Region, 1982-1996, Ph.D. thesis,University of Victoria, 1997. → page 64Bloch, W., M. G. Bostock, and P. Audet, A Cascadia Slab Model From Receiver Functions,Geochemistry, Geophysics, Geosystems, 24(10), e2023GC011,088,doi:https://doi.org/10.1029/2023GC011088, e2023GC011088 2023GC011088, 2023. → pages51, 57, 59Bostock, M., A. Plourde, D. Drolet, and G. Littel, Multichannel Alignment of S Waves, Bulletin ofthe Seismological Society of America, 112(1), 133–142, doi:10.1785/0120210076, 2021. →pages 16, 40Bostock, M. G., N. I. Christensen, and S. M. Peacock, Seismicity in Cascadia, Lithos, 332-333,55–66, doi:https://doi.org/10.1016/j.lithos.2019.02.019, 2019. → page 58Boyarko, D. C., and M. R. Brudzinski, Spatial and Temporal Patterns of Nonvolcanic TremorAlong the Southern Cascadia Subduction Zone, Journal of Geophysical Research: Solid Earth,115(B8), doi:https://doi.org/10.1029/2008JB006064, 2010. → page 58Brantut, N., I. Stefanou, and J. Sulem, Dehydration-Induced Instabilities at Intermediate Depths inSubduction Zones, Journal of Geophysical Research: Solid Earth, 122(8), 6087–6107,doi:https://doi.org/10.1002/2017JB014357, 2017. → page 32Brothers, D. S., N. C. Miller, J. V. Barrie, P. J. Haeussler, H. G. Greene, B. D. Andrews, O. Zielke,J. Watt, and P. Dartnell, Plate Boundary Localization, Slip-Rates and Rupture Segmentation ofthe Queen Charlotte Fault Based on Submarine Tectonic Geomorphology, Earth and PlanetaryScience Letters, 530, 115,882, doi:https://doi.org/10.1016/j.epsl.2019.115882, 2020. → pages63, 65, 67Brown, E., Obducted Nappe Sequence in the San Juan Islands – Northwest Cascades ThrustSystem, Washington and British Columbia, Canadian Journal of Earth Sciences, 49(7),796–817, doi:10.1139/e2012-026, 2012. → pages 52, 53Bustin, A. M. M., R. D. Hyndman, H. Kao, and J. F. Cassidy, Evidence for Underthrusting beneaththe Queen Charlotte Margin, British Columbia, from Teleseismic Receiver Function Analysis,Geophysical Journal International, 171(3), 1198–1211,doi:10.1111/j.1365-246X.2007.03583.x, 2007. → pages 64, 68, 69, 127, 133, 135, 136Be´rube´, J., G. C. Rogers, R. M. Ellis, and E. O. Hasselgren, A Microseismicity Study of the QueenCharlotte Islands Region, Canadian Journal of Earth Sciences, 26(12), 2556–2566,doi:10.1139/e89-218, 1989. → page 6481Cakir, R., Monitoring Active Faults for Tectonic Mapping Efforts in Washington State [Data set].International Federation of Digital Seismograph Networks, doi:10.7914/SN/X4\ 2016, 2016. →page 110Calvert, A. J., M. G. Bostock, G. Savard, and M. J. Unsworth, Cascadia Low FrequencyEarthquakes at the Base of an Overpressured Subduction Shear Zone, Nature Communications,11(3874), doi:https://doi.org/10.1038/s41467-020-17609-3, 2020. → page 58Cassidy, J. F., G. C. Rogers, and R. D. Hyndman, An Overview of the 28 October 2012 Mw 7.7Earthquake in Haida Gwaii, Canada: A Tsunamigenic Thrust Event Along a PredominantlyStrike-Slip Margin, Pure and Applied Geophysics, 171, 3457–3465,doi:10.1007/s00024-014-0775-1, 2014. → page 63Cheng, Y., R. M. Allen, and T. Taira, A New Focal Mechanism Calculation Algorithm (REFOC)Using Inter-Event Relative Radiation Patterns: Application to the Earthquakes in the ParkfieldArea, Journal of Geophysical Research: Solid Earth, 128(3), e2022JB025,006,doi:https://doi.org/10.1029/2022JB025006, e2022JB025006 2022JB025006, 2023. → pages7, 37Chernak, L. J., and G. Hirth, Syndeformational Antigorite Dehydration Produces Stable Fault Slip,Geology, 39(9), 847–850, doi:10.1130/G31919.1, 2011. → page 11Choi, J., J. Byun, S. J. Seol, and S. K. Lee, Convolutional Neural Network-Based Moment TensorInversion Using Domain Adaptation for Microseismicity Monitoring, Exploration Geophysics,54(2), 133–143, doi:10.1080/08123985.2022.2086798, 2023. → page 8Choy, G. L., J. L. Boatwright, and S. H. Kirby, The Radiated Seismic Energy and Apparent Stressof Interplate and Intraplate Earthquakes at Subduction Zone Environments; Implications forSeismic Hazard Estimation, Tech. rep., U.S. Dept. of the Interior, U.S. Geological Survey,doi:10.3133/ofr015, 2001. → page 11Creager, K., J. Vidale, and S. Malone, Cascadia Array of Arrays [Data set]. InternationalFederation of Digital Seismograph Networks, doi:10.7914/SN/XG\ 2009, 2009. → page 110Crotwell, H. P., T. J. Owens, and J. Ritsema, The TauP Toolkit: Flexible Seismic Travel-Time andRay-Path Utilities, Seismological Research Letters, 70(2), 154–160, doi:10.1785/gssrl.70.2.154,1999. → pages 53, 69Cui, Y., D. Miller, P. Schiarizza, and L. Diakow, British Columbia Digital Geology, BritishColumbia Ministry of Energy, Mines and Petroleum Resources, British Columbia GeologicalSurvey Open File 2017-8, 9p., 2017. → pages 67, 131, 134Dahm, T., Relative Moment Tensor Inversion to Determine the Radiation Pattern of SeismicSources, Ph.D. thesis, Geophysikalisches Institut, Universitat Karlsruhe, Germany, (in German),1993. → page 39Dahm, T., Relative Moment Tensor Inversion Based on Ray Theory: Theory and Synthetic Tests,Geophysical Journal International, 124(1), 245–257, doi:10.1111/j.1365-246X.1996.tb06368.x,1996. → pages 8, 16, 3682DeMets, C., and S. Merkouriev, High-Resolution Reconstructions of Pacific–North America PlateMotion: 20 Ma to Present, Geophysical Journal International, 207(2), 741–773,doi:10.1093/gji/ggw305, 2016. → page 63DeMets, C., R. G. Gordon, and D. F. Argus, Geologically Current Plate Motions, GeophysicalJournal International, 181(1), 1–80, doi:https://doi.org/10.1111/j.1365-246X.2009.04491.x,2010. → page 63Drolet, D., M. G. Bostock, A. P. Plourde, and C. G. Sammis, Aftershock Distributions, MomentTensors and Stress Evolution of the 2016 Iniskin and 2018 Anchorage Mw 7.1 Alaskan IntraslabEarthquakes, Geophysical Journal International, 231(1), 199–214, doi:10.1093/gji/ggac165,2022. → pages 20, 21, 37, 39, 48, 95Drolet, D., M. Bostock, and S. Peacock, Relative Moment Tensor Inversion for Microseismicity:Application to Clustered Earthquakes in the Cascadia Forearc, Seismica, 2(4),doi:10.26443/seismica.v2i4.1311, 2024. → pages 66, 68Dziewonski, A. M., T.-A. Chou, and J. H. Woodhouse, Determination of Earthquake SourceParameters from Waveform Data for Studies of Global and Regional Seismicity, Journal ofGeophysical Research: Solid Earth, 86(B4), 2825–2852,doi:https://doi.org/10.1029/JB086iB04p02825, 1981. → page 17D’Amico, S., B. Orecchio, D. Presti, A. Gervasi, L. Zhui, I. Guerra, G. Neri, and R. Herrmann,Testing the Stability of Moment Tensor Solutions for Small and Moderate Earthquakes in theCalabrian-Peloritan Arc Region (South Italy), Bollettino di Geofisica Teorica ed Applicata,52(283-298), doi:10.4430/bgta0009, 2011. → page 17Egbert, G. D., B. Yang, P. A. Bedrosian, K. Key, D. W. Livelybrooks, A. Schultz, A. Kelbert, andB. Parris, Fluid Transport and Storage in the Cascadia Forearc Influenced by Overriding PlateLithology, Nature Geoscience, 15, 2022. → page 58Ekstro¨m, G., M. Nettles, and A. Dziewon´ski, The Global CMT Project 2004–2010:Centroid-Moment Tensors for 13,017 Earthquakes, Physics of the Earth and Planetary Interiors,200-201, 1–9, doi:https://doi.org/10.1016/j.pepi.2012.04.002, 2012. → pages 7, 17, 36Elliott, J. L., C. F. Larsen, J. T. Freymueller, and R. J. Motyka, Tectonic Block Motion and GlacialIsostatic Adjustment in Southeast Alaska and Adjacent Canada Constrained by GPSMeasurements, Journal of Geophysical Research: Solid Earth, 115(B9),doi:https://doi.org/10.1029/2009JB007139, 2010. → page 15Etchecopar, A., G. Vasseur, and M. Daignieres, An Inverse Problem in Microtectonics for theDetermination of Stress Tensors from Fault Striation Analysis, Journal of Structural Geology,3(1), 51–65, doi:10.1016/0191-8141(81)90056-0, 1981. → page 70Farahbod, A. M., and J. Cassidy., Temporal Variations in Coda Q before and after the 2017 BarrowStrait earthquake (Mw 5.9) in Nunavut and the 2012 Haida Gwaii Earthquake (Mw 7.8) inBritish Columbia, Tech. rep., Geological Survey of Canada, Open file, 2816-7155 ; 8930, 2022.→ pages 68, 13383Farahbod, A. M., and H. Kao, Spatiotemporal Distribution of Events during the First Week of the2012 Haida Gwaii Aftershock Sequence, Bulletin of the Seismological Society of America,105(2B), 1231–1240, doi:10.1785/0120140173, 2015. → page 137Ferrand, T. P., et al., Dehydration-Driven Stress Transfer Triggers Intermediate-Depth Earthquakes,Nature Communications, 8(1), 15,247, doi:10.1038/ncomms15247, 2017. → pages 11, 32Fuenzalida, A., B. Schurr, M. Lancieri, M. Sobiesiak, and R. Madariaga, High-ResolutionRelocation and Mechanism of Aftershocks of the 2007 Tocopilla (Chile) Earthquake,Geophysical Journal International, 194(2), 1216–1228, doi:10.1093/gji/ggt163, 2013. → page30Fyfe, W. S., R. Kerrich, S. R. Hicock, and F. C. Colloza, LITHOPROBE: Faults and Fluids,Geoscience Canada, 14(2), 1987. → page 60Gasc, J., A. Schubnel, F. Brunet, S. Guillon, H.-J. Mueller, and C. Lathe, Simultaneous AcousticEmissions Monitoring and Synchrotron X-Ray Diffraction at High Pressure and Temperature:Calibration and Application to Serpentinite Dehydration, Physics of the Earth and PlanetaryInteriors, 189(3), 121–133, doi:https://doi.org/10.1016/j.pepi.2011.08.003, 2011. → pages11, 32GEBCO, GEBCO 2023 Grid, GEBCO Compilation Group,doi:10.5285/f98b053b-0cbc-6c23-e053-6c86abc0af7b, 2023. → page 65Geological Survey of Canada, Canadian Seismic Research Network, doi:10.7914/601S-AK68,2002. → page 124Gosselin, J. M., J. F. Cassidy, and S. E. Dosso, Shear-Wave Velocity Structure in the Vicinity of the2012 Mw 7.8 Haida Gwaii Earthquake from Receiver Function Inversion, Bulletin of theSeismological Society of America, 105(2B), 1106–1113, doi:10.1785/0120140171, 2015. →pages 64, 136Greene, H. G., and J. V. Barrie, Faulting within the San Juan–southern Gulf Islands Archipelagos,Upper Plate Deformation of the Cascadia Subduction Complex, in From Continental Shelf toSlope: Mapping the Oceanic Realm, Geological Society of London,doi:10.1144/SP505-2019-125, 2022. → pages 59, 60Greene, H. G., J. V. Barrie, and B. J. Todd, The Skipjack Island Fault Zone: An ActiveTranscurrent Structure within the Upper Plate of the Cascadia Subduction Complex,Sedimentary Geology, 378, 61–79, doi:https://doi.org/10.1016/j.sedgeo.2018.05.005, 2018. →page 59Guo, Y., K. Miyakoshi, and M. Tsurugi, Simultaneous Rupture on Conjugate Faults during the2018 Anchorage, Alaska, Intraslab Earthquake (MW 7.1) Inverted from Strong-MotionWaveforms, Earth, Planets and Space, 72(1), 176, doi:10.1186/s40623-020-01315-x, 2020. →pages 15, 32, 33, 76Hacker, B. R., S. M. Peacock, G. A. Abers, and S. D. Holloway, Subduction Factory 2. AreIntermediate-Depth Earthquakes in Subducting Slabs Linked to Metamorphic Dehydration84Reactions?, Journal of Geophysical Research: Solid Earth, 108(B1),doi:https://doi.org/10.1029/2001JB001129, 2003. → page 11Hanks, T. C., and H. Kanamori, A Moment Magnitude Scale, Journal of Geophysical Research:Solid Earth, 84(B5), 2348–2350, doi:https://doi.org/10.1029/JB084iB05p02348, 1979. → page47Hardebeck, J. L., and E. Hauksson, Crustal Stress Field in Southern California and its Implicationsfor Fault Mechanics, Journal of Geophysical Research: Solid Earth, 106(B10), 21,859–21,882,doi:https://doi.org/10.1029/2001JB000292, 2001. → pages 18, 19Hardebeck, J. L., and P. M. Shearer, A New Method for Determining First-Motion FocalMechanisms, Bulletin of the Seismological Society of America, 92(6), 2264–2276,doi:10.1785/0120010200, 2002. → page 16Hardebeck, J. L., and P. M. Shearer, Using S/P Amplitude Ratios to Constrain the FocalMechanisms of Small Earthquakes, Bulletin of the Seismological Society of America, 93(6),2434–2444, doi:10.1785/0120020236, 2003. → pages 7, 36Hasegawa, A., and J. Nakajima, Seismic Imaging of Slab Metamorphism and Genesis ofIntermediate-Depth Intraslab Earthquakes, Progress in Earth and Planetary Science, 4(1), 12,doi:10.1186/s40645-017-0126-9, 2017. → page 11Hayes, G., Slab2 - A Comprehensive Subduction Zone Geometry Model, U.S. Geological Surveydata release, doi:https://doi.org/10.5066/F7PV6JNV, 2018. → pages 19, 102Hayward, T. W., and M. G. Bostock, Slip Behavior of the Queen Charlotte Plate Boundary beforeand after the 2012, MW 7.8 Haida Gwaii Earthquake: Evidence from Repeating Earthquakes,Journal of Geophysical Research: Solid Earth, 122(11), 8990–9011,doi:https://doi.org/10.1002/2017JB014248, 2017. → page 48Healy, D., Anisotropic Poroelasticity and the Response of Faulted Rock to Changes in Pore-FluidPressure, Geological Society, London, Special Publications, 367(1), 201–214,doi:10.1144/SP367.14, 2012. → pages 31, 105, 108Helffrich, G. R., How Good are Routinely Determined Focal Mechanisms? Empirical StatisticsBased on a Comparison of Harvard, USGS and ERI Moment Tensors, Geophysical JournalInternational, 131(3), 741–750, doi:10.1111/j.1365-246X.1997.tb06609.x, 1997. → page 16Hill, D. P., A Model for Earthquake Swarms, Journal of Geophysical Research (1896-1977), 82(8),1347–1352, doi:https://doi.org/10.1029/JB082i008p01347, 1977. → pages 52, 135Hobbs, T. E., J. F. Cassidy, and S. E. Dosso, Rupture process of the 2012 mw 7.8 haida gwaiiearthquake from an empirical green’s function method, Bulletin of the Seismological Society ofAmerica, 105(2B), 1219–1230, doi:10.1785/0120140175, 2015. → page 74Hyndman, R. D., Tectonics and Structure of the Queen Charlotte Fault Zone, Haida Gwaii, andLarge Thrust Earthquakes, Bulletin of the Seismological Society of America, 105(2B),1058–1075, doi:10.1785/0120140181, 2015. → page 6485Hyndman, R. D., and T. S. Hamilton, Queen Charlotte Area Cenozoic Tectonics and Volcanismand their Association with Relative Plate Motions along the Northeastern Pacific Margin,Journal of Geophysical Research: Solid Earth, 98(B8), 14,257–14,277,doi:https://doi.org/10.1029/93JB00777, 1993. → page 63Hyndman, R. D., T. J. Lewis, J. A. Wright, M. Burgess, D. S. Chapman, and M. Yamano, QueenCharlotte Fault Zone: Heat Flow Measurements, Canadian Journal of Earth Sciences, 19(8),1657–1669, doi:10.1139/e82-141, 1982. → page 63Ingersoll, L. R., O. J. Zobel, and A. C. Ingersoll, Heat Conduction: with Engineering andGeological Applications, 1st ed., McGraw-Hill Book Co, New York, 1948. → pages 104, 105IRIS Transportable Array, USArray Transportable Array [Data Set], International Federation ofDigital Seismograph Networks, doi:10.7914/SN/TA, 2003. → page 110Ishii, M., E. Kiser, and E. L. Geist, Mw 8.6 Sumatran Earthquake of 11 April 2012: Rare SeawardExpression of Oblique Subduction, Geology, 41(3), 319–322, doi:10.1130/G33783.1, 2013. →page 63Izutani, Y., and H. Kanamori, Scale-Dependence of Seismic Energy-to-Moment Ratio forStrike-Slip Earthquakes in Japan, Geophysical Research Letters, 28(20), 4007–4010,doi:https://doi.org/10.1029/2001GL013402, 2001. → page 11Jadamec, M. A., and M. I. Billen, Reconciling Surface Plate Motions with RapidThree-Dimensional Mantle Flow around a Slab Edge, Nature, 465(7296), 338–341,doi:10.1038/nature09053, 2010. → pages 13, 14, 19, 20, 102Jarrard, R. D., Relations among Subduction Parameters, Reviews of Geophysics, 24(2), 217–284,doi:https://doi.org/10.1029/RG024i002p00217, 1986. → page 62John, T., S. Medvedev, L. H. Ru¨pke, T. B. Andersen, Y. Y. Podladchikov, and H. Austrheim,Generation of Intermediate-Depth Earthquakes by Self-Localizing Thermal Runaway, NatureGeoscience, 2(2), 137–140, doi:10.1038/ngeo419, 2009. → page 12Johnson, H. P., S. G. Merle, T. A. Bjorklund, S. L. Hautala, T. Baumberger, S. L. Walker, J. Liu,N. D. Ward, and C. Wang, Methane Plume Emissions Associated With Puget Sound Faults inthe Cascadia Forearc, Geochemistry, Geophysics, Geosystems, 23(1),doi:https://doi.org/10.1029/2021GC010211, 2022. → page 60Johnson, S., S. Dadisman, D. Mosher, R. Blakely, and J. Childs, Active Tectonics of the Devil’sMountain Fault and Related Structures, Northern Puget Lowland and Eastern Strait of Juan deFuca Region, Pacific Northwest, Tech. rep., U.S. Geological Survey Professional Paper 1643,2001. → page 59Jung, H., H. W. Green II, and L. F. Dobrzhinetskaya, Intermediate-Depth Earthquake Faulting byDehydration Embrittlement with Negative Volume Change, Nature, 428(6982), 545–549,doi:10.1038/nature02412, 2004. → pages 11, 3286Kagan, Y. Y., 3-D Rotation of Double-Couple Earthquake Sources, Geophysical JournalInternational, 106(3), 709–716, doi:10.1111/j.1365-246X.1991.tb06343.x, 1991. → pages17, 47, 69Kao, H., S.-J. Shan, H. Dragert, and G. Rogers, Northern Cascadia Episodic Tremor and Slip: ADecade of Tremor Observations from 1997 to 2007, Journal of Geophysical Research: SolidEarth, 114(B11), doi:https://doi.org/10.1029/2008JB006046, 2009. → page 58Kao, H., S.-J. Shan, A. Bent, C. Woodgold, G. Rogers, J. F. Cassidy, and J. Ristau, RegionalCentroid-Moment-Tensor Analysis for Earthquakes in Canada and Adjacent Regions: AnUpdate, Seismological Research Letters, 83(3), 505–515, doi:10.1785/gssrl.83.3.505, 2012. →pages 7, 64, 65, 67, 68, 69, 133Kao, H., S. Shan, and A. M. Farahbod, Source Characteristics of the 2012 Haida Gwaii EarthquakeSequence, Bulletin of the Seismological Society of America, 105(2B), 1206–1218,doi:10.1785/0120140165, 2015. → page 137Kelemen, P. B., and G. Hirth, A Periodic Shear-Heating Mechanism for Intermediate-Depthearthquakes in the mantle, Nature, 446(7137), 787–790, doi:10.1038/nature05717, 2007. →page 12Kim, Y., G. A. Abers, J. Li, D. Christensen, J. Calkins, and S. Rondenay, Alaska Megathrust 2:Imaging the Megathrust Zone and Yakutat/Pacific Plate Interface in the Alaska SubductionZone, Journal of Geophysical Research: Solid Earth, 119(3), 1924–1941,doi:https://doi.org/10.1002/2013JB010581, 2014. → page 15Kirby, S., Interslab Earthquakes and Phase Changes in Subducting Lithosphere, Reviews ofGeophysics, 33(S1), 287–297, doi:https://doi.org/10.1029/95RG00353, 1995. → page 11Kirby, S. H., K. Wang, and S. Dunlop, The Cascadia Subduction Zone and Related SubductionSystems: Seismic structure, Intraslab Earthquakes and Processes, and Earthquake Hazards,Tech. rep., U.S. Geological Suvey, Reston, VA, doi:10.3133/ofr02328, 2002. → page 11Klein, F., User’s Guide to HYPOINVERSE-2000, a Fortran Program to Solve for EarthquakeLocations and Magnitudes, Version 1.40, June 2014, Tech. rep., U. S. Geological Survey,doi:10.13140/2.1.4859.3602, 2014. → page 15Li, Z., Z. Peng, Y. Ben-Zion, and F. L. Vernon, Spatial Variations of Shear Wave Anisotropy nearthe San Jacinto Fault Zone in Southern California, Journal of Geophysical Research: SolidEarth, 120(12), 8334–8347, doi:https://doi.org/10.1002/2015JB012483, 2015. → page 31Liu, C., T. Lay, Z. Xie, and X. Xiong, Intraslab Deformation in the 30 November 2018 Anchorage,Alaska, MW 7.1 Earthquake, Geophysical Research Letters, 46(5), 2449–2457,doi:https://doi.org/10.1029/2019GL082041, 2019. → page 15Lund, B., and R. Slunga, Stress Tensor Inversion Using Detailed Microearthquake Information andStability Constraints: Application to O¨lfus in Southwest Iceland, Journal of GeophysicalResearch: Solid Earth, 104(B7), 14,947–14,964, doi:https://doi.org/10.1029/1999JB900111,1999. → page 1887Malone, S., Seismic Array Studies of Cascadia Deep Tremor [Data set]. International Federation ofDigital Seismograph Networks, doi:10.7914/SN/XH\ 2004, 2004. → page 110Malone, S., K. Creager, S. Rondenay, T. Melbourne, and G. Abers, Collaborative Research:Earthscope Integrated Investigations of Cascadia Subduction Zone Tremor, Structure andProcess [Data Set]. International Federation of Digital Seismograph Networks,doi:10.7914/SN/XU\ 2006, 2006. → page 110Marot, M., T. Monfret, M. Pardo, G. Ranalli, and G. Nolet, An Intermediate-Depth TensionalEarthquake (MW 5.7) and its Aftershocks within the Nazca Slab, Central Chile: A ReactivatedOuter Rise Fault?, Earth and Planetary Science Letters, 327-328, 9–16,doi:https://doi.org/10.1016/j.epsl.2012.02.003, 2012. → page 30Martı´nez-Garzo´n, P., Y. Ben-Zion, N. Abolfathian, G. Kwiatek, and M. Bohnhoff, A RefinedMethodology for Stress Inversions of Earthquake Focal Mechanisms, Journal of GeophysicalResearch: Solid Earth, 121(12), 8666–8687, doi:https://doi.org/10.1002/2016JB013493, 2016.→ pages 18, 31, 70McCaffrey, R., M. D. Long, C. Goldfinger, P. C. Zwick, J. L. Nabelek, C. K. Johnson, andC. Smith, Rotation and Plate Locking at the Southern Cascadia Subduction Zone, GeophysicalResearch Letters, 27(19), 3117–3120, doi:https://doi.org/10.1029/2000GL011768, 2000. →page 57McCrory, P. A., J. E. Constantz, A. G. Hunt, and J. L. Blair, Helium as a Tracer for Fluids Releasedfrom Juan de Fuca Lithosphere Beneath the Cascadia Forearc, Geochemistry, Geophysics,Geosystems, 17(6), 2434–2449, doi:https://doi.org/10.1002/2015GC006198, 2016. → page 60Merrill, R., M. G. Bostock, S. M. Peacock, A. J. Calvert, and N. I. Christensen, A DoubleDifference Tomography Study of the Washington Forearc: Does Siletzia Control CrustalSeismicity?, Journal of Geophysical Research: Solid Earth, 125(10), e2020JB019,750,doi:https://doi.org/10.1029/2020JB019750, e2020JB019750 2020JB019750, 2020. → page 51Merrill, R. J., M. G. Bostock, S. M. Peacock, A. J. Schaeffer, and S. W. Roecker, ComplexStructure in the Nootka Fault Zone Revealed by Double-Difference Tomography and a NewEarthquake Catalog, Geochemistry, Geophysics, Geosystems, 23(2), e2021GC010,205,doi:https://doi.org/10.1029/2021GC010205, e2021GC010205 2021GC010205, 2022. → pages51, 53, 58Michael, A. J., Determination of Stress from Slip Data: Faults and Folds, Journal of GeophysicalResearch: Solid Earth, 89(B13), 11,517–11,526, doi:https://doi.org/10.1029/JB089iB13p11517,1984. → pages 8, 69Michael, A. J., Spatial Variations in Stress within the 1987 Whittier Narrows, California,Aftershock Sequence: New Techniques and Results, Journal of Geophysical Research: SolidEarth, 96(B4), 6303–6319, doi:https://doi.org/10.1029/91JB00195, 1991. → page 19Morell, K. D., C. Regalla, L. J. Leonard, C. Amos, and V. Levson, Quaternary Rupture of a CrustalFault beneath Victoria, British Columbia, Canada, GSA Today, 7, doi:10.1130/GSATG291A.1,2017. → page 5988Mu¨ller, R. D., W. R. Roest, J.-Y. Royer, L. M. Gahagan, and J. G. Sclater, Digital Isochrons of theWorld’s Ocean Floor, Journal of Geophysical Research: Solid Earth, 102(B2), 3211–3214,doi:https://doi.org/10.1029/96JB01781, 1997. → page 13Nakajima, J., A. Hasegawa, and S. Kita, Seismic Evidence for Reactivation of a Buried HydratedFault in the Pacific Slab by the 2011 M9.0 Tohoku Earthquake, Geophysical Research Letters,38(7), doi:https://doi.org/10.1029/2011GL048432, 2011. → page 30Natural Resources Canada, Canadian National Seismograph Network [Data set]., Tech. rep.,International Federation of Digital Seismograph Networks, doi:10.7914/SN/CN, 1975. → pages110, 124Natural Resources Canada, Search the Earthquake Database, Tech. rep., NRCan,https://www.earthquakescanada.nrcan.gc.ca/stndon/NEDB-BNDS/bulletin-en.php, accessedJuly 9, 2024, 2021a. → page 68Natural Resources Canada, What Is Shown on a Seismogram?,https://www.earthquakescanada.nrcan.gc.ca/info-gen/smeters-smetres/seismogram-en.phpaccessed on August 13, 2024, 2021b. → page 4Nayak, A., D. Eberhart-Phillips, N. A. Ruppert, H. Fang, M. M. Moore, C. Tape, D. H.Christensen, G. A. Abers, and C. H. Thurber, 3D Seismic Velocity Models for Alaska from JointTomographic Inversion of Body-Wave and Surface-Wave Data, Seismological Research Letters,91(6), 3106–3119, doi:10.1785/0220200214, 2020. → page 16NOAA National Oceanic and Atmospheric Administration (USA), National Tsunami WarningCenter Alaska Seismic Network, doi:10.7914/SN/AT, 1967. → page 124Oliva, S. J., et al., Incipient Subduction and Slip Partitioning at High Obliquity: The Haida GwaiiPlate Boundary, Journal of Geophysical Research: Solid Earth, 129(5), e2024JB028,752,doi:https://doi.org/10.1029/2024JB028752, e2024JB028752 2024JB028752, 2024. → pages63, 64, 65, 66, 68, 69, 74, 126, 132, 133, 135, 136, 137Pacific Northwest Seismic Network, Tremor Map, Tech. rep., Pacific Northwest Seismic Network,https://pnsn.org/tremor, accessed November 26, 2023, 2009-2021. → page 52Paszke, A., et al., PyTorch: An Imperative Style, High-Performance Deep Learning Library, inNeural Information Processing Systems, 2019. → page 46Peacock, S. M., Are the Lower Planes of Double Seismic Zones Caused by SerpentineDehydration in Subducting Oceanic Mantle?, Geology, 29(4), 299–302,doi:10.1130/0091-7613(2001)029⟨0299:ATLPOD⟩2.0.CO;2, 2001. → pages 11, 12, 29, 137Peacock, S. M., and K. Wang, Seismic Consequences of Warm Versus Cool SubductionMetamorphism: Examples from Southwest and Northeast Japan, Science, 286(5441), 937–939,doi:10.1126/science.286.5441.937, 1999. → page 11Plourde, A. P., and M. G. Bostock, Relative Moment Tensors and Deep Yakutat Seismicity,Geophysical Journal International, 219(2), 1447–1462, doi:10.1093/gji/ggz375, 2019. → pages8, 9, 16, 36, 37, 38, 39, 41, 43, 46, 48, 60, 66, 68, 77, 7889Pondrelli, S., S. Salimbeni, G. Ekstro¨m, A. Morelli, P. Gasperini, and G. Vannucci, The ItalianCMT Dataset from 1977 to the Present, Physics of the Earth and Planetary Interiors, 159(3),286–303, doi:https://doi.org/10.1016/j.pepi.2006.07.008, 2006. → page 47Preston, L. A., K. C. Creager, R. S. Crosson, T. M. Brocher, and A. M. Trehu, IntraslabEarthquakes: Dehydration of the Cascadia Slab, Science, 302(5648), 1197–1200,doi:10.1126/science.1090751, 2003. → page 11Prieto, G. A., M. Florez, S. A. Barrett, G. C. Beroza, P. Pedraza, J. F. Blanco, and E. Poveda,Seismic Evidence for Thermal Runaway during Intermediate-Depth Earthquake Rupture,Geophysical Research Letters, 40(23), 6064–6068, doi:https://doi.org/10.1002/2013GL058109,2013. → page 12Proctor, B., and G. Hirth, Role of Pore Fluid Pressure on Transient Strength Changes and FabricDevelopment during Serpentine Dehydration at Mantle Conditions: Implications forSubduction-Zone Seismicity, Earth and Planetary Science Letters, 421, 1–12,doi:https://doi.org/10.1016/j.epsl.2015.03.040, 2015. → page 11Ramachandran, K., R. D. Hyndman, and T. M. Brocher, Regional P Wave Velocity Structure of theNorthern Cascadia Subduction Zone, Journal of Geophysical Research: Solid Earth, 111(B12),doi:https://doi.org/10.1029/2005JB004108, 2006. → page 58Rice, J. R., Heating and Weakening of Faults during Earthquake Slip, Journal of GeophysicalResearch: Solid Earth, 111(B5), doi:https://doi.org/10.1029/2005JB004006, 2006. → page 32Riedel, M., S. Yelisetti, C. Papenberg, K. Rohr, M. Coˆte´, G. Spence, R. Hyndman, and T. James,Seismic Velocity Structure of the Queen Charlotte Terrace off Western Canada in the Region ofthe 2012 Haida Gwaii Mw 7.8 Thrust Earthquake, Geosphere, 17(1), 23–38,doi:10.1130/GES02258.1, 2020. → page 64Rippe, D., M. J. Unsworth, and C. A. Currie, Magnetotelluric Constraints on the Fluid Content inthe Upper Mantle beneath the Southern Canadian Cordillera: Implications for Rheology,Journal of Geophysical Research: Solid Earth, 118(10), 5601–5624,doi:https://doi.org/10.1002/jgrb.50255, 2013. → page 58Ristau, J., G. C. Rogers, and J. F. Cassidy, Stress in Western Canada from Regional MomentTensor Analysis, Canadian Journal of Earth Sciences, 44(2), 127–148, doi:10.1139/e06-057,2007. → pages 64, 65Rogers, G. C., Seismotectonics of British Columbia, Ph.D. thesis, University of British Columbia,doi:http://dx.doi.org/10.14288/1.0052938, 1983. → page 65Rowan, C., 5 Focal Mechanisms,https://all-geo.org/highlyallochthonous/2009/12/5-focal-mechanisms/ accessed on August 13,2024, 2009. → page 6Rubinstein, J. L., M. L. Rocca, J. E. Vidale, K. C. Creager, and A. G. Wech, Tidal Modulation ofNonvolcanic Tremor, Science, 319(5860), 186–189, doi:10.1126/science.1150558, 2008. →page 5890Ruppert, N. A., A. Nayak, C. Thurber, and C. Richards, Aftershock Analysis of the 2018 Mw 7.1Anchorage, Alaska, Earthquake: Relocations and Regional Moment Tensors, SeismologicalResearch Letters, 91(1), 114–125, doi:10.1785/0220190199, 2019. → pages 15, 20, 32, 33, 96Savard, G., M. G. Bostock, and N. I. Christensen, Seismicity, Metamorphism, and Fluid EvolutionAcross the Northern Cascadia Fore Arc, Geochemistry, Geophysics, Geosystems, 19(6),1881–1897, doi:https://doi.org/10.1029/2017GC007417, 2018. → page 58Shearer, P. M., Introduction to Seismology, chap. 9, Cambridge University Press,doi:https://doi.org/10.1017/CBO9780511841552, 2009. → pages 5, 36Shearer, P. M., G. A. Prieto, and E. Hauksson, Comprehensive Analysis of Earthquake SourceSpectra in Southern California, Journal of Geophysical Research: Solid Earth, 111(B6),doi:https://doi.org/10.1029/2005JB003979, 2006. → page 47Shelly, D. R., R. J. Skoumal, and J. L. Hardebeck, S/P Amplitude Ratios Derived fromSingle-Component Seismograms and Their Potential Use in Constraining Focal Mechanisms forMicroearthquake Sequences, The Seismic Record, 2(2), 118–126, doi:10.1785/0320220002,2022. → page 36Shiina, T., J. Nakajima, and T. Matsuzawa, Seismic Evidence for High Pore Pressures in theOceanic Crust: Implications for Fluid-Related Embrittlement, Geophysical Research Letters,40(10), 2006–2010, doi:https://doi.org/10.1002/grl.50468, 2013. → page 11Shuck, B., et al., Strike-Slip Enables Subduction Initiation Beneath a Failed Rift: New SeismicConstraints From Puysegur Margin, New Zealand, Tectonics, 40(5), e2020TC006,436,doi:https://doi.org/10.1029/2020TC006436, e2020TC006436 2020TC006436, 2021. → page 63Silver, P., Active Seismic Study of Episodic Tremor and Slip, Detecting Structural changes Duringan ETS Event: Proof of Concept [Data set]. International Federation of Digital SeismographNetworks, doi:10.7914/SN/Z3\ 2009, 2009. → page 110Smith, A. J., R. D. Hyndman, J. F. Cassidy, and K. Wang, Structure, Seismicity, and ThermalRegime of the Queen Charlotte Transform Margin, Journal of Geophysical Research: SolidEarth, 108(B11), doi:https://doi.org/10.1029/2002JB002247, 2003. → page 64Souther, J., and A. Jessop, Dyke swarms in the Queen Charlotte Islands, British Columbia, andImplications for Hydrocarbon Exploration, Evolution and Hydrocarbon Potential of the QueenCharlotte Basin, British Columbia, Geological Survey of Canada, Paper 90-10, 1991. → page136Stanton-Yonge, A., W. A. Griffith, J. Cembrano, R. St. Julien, and P. Iturrieta, Tectonic Role ofMargin-Parallel and Margin-Transverse Faults during Oblique Subduction in the SouthernVolcanic Zone of the Andes: Insights from Boundary Element Modeling, Tectonics, 35(9),1990–2013, doi:https://doi.org/10.1002/2016TC004226, 2016. → page 63Sutherland, B. A., Geology of the Queen Charlotte Islands, British Columbia, vol. 54, BritishColumbia Dept. of Mines and Petroleum Resources, Victoria, 1968. → pages 136, 13791Thompson, R. I., Geology, Cumshewa Inlet, British Columbia, Geological Survey of Canada,doi:https://doi.org/10.4095/131495, Preliminary Map, 1990. → pages 67, 131, 136, 137Tre´hu, A. M., M. Scheidhauer, K. M. M. Rohr, B. Tikoff, M. A. L. Walton, S. P. S. Gulick, andE. C. Roland, An Abrupt Transition in the Mechanical Response of the Upper Crust toTranspression along the Queen Charlotte Fault, Bulletin of the Seismological Society of America,105(2B), 1114–1128, doi:10.1785/0120140159, 2015. → page 63University of Washington, Pacific Northwest Seismic Network - University of Washington [Dataset]., International Federation of Digital Seismograph Networks, doi:10.7914/SN/UW, 1963. →page 110U.S. Geological Survey, United States National Strong-Motion Network [Data set]. , InternationalFederation of Digital Seismograph Networks, doi:10.7914/SN/NP, 1931. → page 110U.S. Geological Survey, Subduction Fault Zone Diagram, Communications and Publishing,https://www.usgs.gov/media/images/subduction-fault-zone-diagram accessed August 7, 2024,2016. → page 4U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/eventpage/ak20421672/executive accessed July 13,2021, 2017a. → pages 21, 95U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, Advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/map/ accessed October 16, 2021, 2017b. → page 101U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, Advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/eventpage/ak0138esnzr/moment-tensor accessedSeptember 16, 2024, 2017c. → page 65U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, Advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/eventpage/ak018fcnsk91/executive accessed June 1,2021, 2017d. → pages 15, 23, 95U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, Advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/eventpage/ak01613v15nv/executive accessed June 1,2021, 2017e. → pages 12, 19, 29, 30, 9592U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, Advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/eventpage/uw10427908/executive accessed October16, 2023, 2017f. → page 51U.S. Geological Survey, Earthquake Hazards Program, Tech. rep., USGS,doi:https://doi.org/10.5066/F7MS3QZH, Advanced National Seismic System (ANSS)Comprehensive Catalog of Earthquake Events and Products: Various,https://earthquake.usgs.gov/earthquakes/map/ accessed October 16, 2021, 2021. → pages30, 95Utsu, T., A Statistical Study on the Occurrence of Aftershocks, The Geophysical Magazine, 30,1961. → page 106VanDecar, J. C., and R. S. Crosson, Determination of Teleseismic Relative Phase Arrival TimesUsing Multi-Channel Cross-Correlation and Least Squares, Bulletin of the Seismological Societyof America, 80(1), 150–169, doi:10.1785/BSSA0800010150, 1990. → page 40Vavrycˇuk, V., Iterative Joint Inversion for Stress and Fault Orientations from Focal Mechanisms,Geophysical Journal International, 199(1), 69–77, doi:10.1093/gji/ggu224, 2014. → pages8, 18, 55, 69, 70Vidale, J. E., and P. M. Shearer, A Survey of 71 Earthquake Bursts Across Southern California:Exploring the Role of Pore Fluid Pressure Fluctuations and Aseismic Slip as Drivers, Journal ofGeophysical Research: Solid Earth, 111(B5), doi:https://doi.org/10.1029/2005JB004034, 2006.→ pages 52, 58, 135Waldhauser, F., and W. L. Ellsworth, A Double-Difference Earthquake Location Algorithm:Method and Application to the Northern Hayward Fault, California, Bulletin of theSeismological Society of America, 90(6), 1353–1368, doi:10.1785/0120000006, 2000. → pages15, 36, 53, 132Wang, K., Unbending Combined with Dehydration Embrittlement as a Cause for Double andTriple Seismic Zones, Geophysical Research Letters, 29(18), 36–1–36–4,doi:https://doi.org/10.1029/2002GL015441, 2002. → page 11Wang, K., J. He, F. Schulzeck, R. D. Hyndman, and M. Riedel, Thermal Condition of the 27October 2012 Mw 7.8 Haida Gwaii Subduction Earthquake at the Obliquely Convergent QueenCharlotte Margin, Bulletin of the Seismological Society of America, 105(2B), 1290–1300,doi:10.1785/0120140183, 2015. → page 136Wannamaker, P. E., R. L. Evans, P. A. Bedrosian, M. J. Unsworth, V. Maris, and R. S. McGary,Segmentation of Plate Coupling, Fate of Subduction Fluids, and Modes of Arc Magmatism inCascadia, Inferred from Magnetotelluric Resistivity, Geochemistry, Geophysics, Geosystems,15(11), 4230–4253, doi:https://doi.org/10.1002/2014GC005509, 2014. → pages 58, 59Warren, L. M., A. N. Hughes, and P. G. Silver, Earthquake Mechanics and Deformation in theTonga-Kermadec Subduction Zone from Fault Plane Orientations of Intermediate- and93Deep-Focus Earthquakes, Journal of Geophysical Research: Solid Earth, 112(B5),doi:https://doi.org/10.1029/2006JB004677, 2007. → page 30Warren, L. M., M. A. Langstaff, and P. G. Silver, Fault Plane Orientations of Intermediate-DepthEarthquakes in the Middle America Trench, Journal of Geophysical Research: Solid Earth,113(B1), doi:https://doi.org/10.1029/2007JB005028, 2008. → page 30Warren-Smith, E., et al., Episodic Dtress and Fluid Pressure Cycling in Subducting Oceanic Crustduring Slow Slip, Nature Geoscience, 12(6), 475–481, doi:10.1038/s41561-019-0367-x, 2019.→ page 31West, M. E., et al., The 30 November 2018 Mw 7.1 Anchorage Earthquake, SeismologicalResearch Letters, 91(1), 66–84, doi:10.1785/0220190176, 2019. → page 15Worthington, L. L., H. J. A. Van Avendonk, S. P. S. Gulick, G. L. Christeson, and T. L. Pavlis,Crustal Structure of the Yakutat Terrane and the Evolution of Subduction and Collision inSouthern Alaska, Journal of Geophysical Research: Solid Earth, 117(B1),doi:https://doi.org/10.1029/2011JB008493, 2012. → page 33Yamasaki, T., and T. Seno, Double Seismic Zone and Dehydration Embrittlement of theSubducting Slab, Journal of Geophysical Research: Solid Earth, 108(B4),doi:https://doi.org/10.1029/2002JB001918, 2003. → page 11Zhang, M., FocMecDR: A Cross-Correlation-Based Double-Ratio Earthquake Focal MechanismInversion, Poster presented at: Seismological Society of America (SSA) Annual Meeting, April2023, San Juan, PR, USA, 2023. → pages 7, 3794Appendix ASupporting Materials for Chapter 2A.1 Data AvailabilityWaveforms used for this study are available form the Incorporated Research Institutions for Seis-mology (IRIS). A few moment tensors from USGS [U.S. Geological Survey, 2017a, d, e, 2021]were employed. Text files with the relocations and recovered MTs for this study are included in theSupplementary Material of the published paper [Drolet et al., 2022].95A.2 FiguresFigure A.1: Anchorage sequence relocation comparison between this study and Ruppert et al.[2019]. Blue circles show the Ruppert et al. [2019] relocations and orange circles showthis study relocations. Only the 337 common events are displayed. Note that Ruppertet al. [2019] included smaller magnitude events and used a smaller time window.96Figure A.2: Choice of HASH solution weights for Iniskin (blue) and Anchorage (orange) se-quences. Solid lines: mean Kagan angle per moment tensor between the solution withoutweight and weighted solution. Dashed lines: same as solid lines but for the subset ofevent with HASH solutions. As weight increases, solutions diverge from the no-weightsolution and but events with HASH solutions are less affected. The preferred solutionwas generated with a weight of 100 near the points of maximum curvature. Note thatHASH solutions were based on a 5° grid search.97Figure A.3: Quality test on moment tensors of the Iniskin sequence. Standard deviation of theKagan angle for each event from bootstrapping. Plot is based on the 70 (of a total 100)iterations that converged to double couple solutions. The plot is saturated at 15◦; redbands are events with a standard deviation of the Kagan angle above 5°. The events aresorted with time; low ids represent earlier origin times. The mainshock id is 44 and wasnot part of the moment tensor inversion.98Figure A.4: Quality test on moment tensors of the Anchorage sequence. Standard deviationof the Kagan angle for each event from bootstrapping. Plot is based on the 70 (of a total120) iterations that converged to double couple solutions. The plot is saturated at 15°;red bands are events with a standard deviation of the Kagan above 5°. The events aresorted with time; low ids represent earlier origin times. The mainshock id is 30 and wasnot part of the moment tensor inversion.99Figure A.5: Reference events used for the Iniskin moment tensor inversion.Figure A.6: Reference events used for the Anchorage moment tensor inversion.100Figure A.7: Intermediate deep events of Mw≥6 with moment tensor solutions (lower hemi-sphere representation) and depths from the USGS catalogue [U.S. Geological Survey,2017b]. Green contour lines: depths to Pacific plate interface between 50 and 250 km at30 km intervals for the E115 model. Purple star: Iniskin mainshock location.101A.3 Slab E115 vs Slab 2 modelsIn the analysis previously presented analysis used the E115 Slab model [Jadamec and Billen, 2010],since it correlates better with the earthquake relocations. However, note that the use of the Slab 2model [Hayes, 2018], would yield similar conclusions. Fig. A.7 displays maps and cross sectionsshowing the difference between the models.In the Iniskin region, the strike angle of the Slab 2 model is∼5° greater than the slab E115 modeland the former is ∼22 km shallower. The angle between the earthquake lineation and intersectionof rupture plane and slab 2 is less than 8° which does not change the interpretation that the lineationis located at the intersection of the slab and the mainshock rupture plane. Slab 2 model should bedeeper to accommodate the relocation and that’s why the E115 slab model is used.There are important differences between the two models in the Anchorage area. However, inthis region, it is hard to define the slab interface due to the poorly defined boundary between thePacific and the Yakutat Plates. Note that the slab model was not used to draw any conclusions aboutthe Anchorage sequence.102Figure A.8: A: Overview of the difference between the E115 Slab (green) and Slab 2 (purple)models. B: Map view of slab models around the 2016 Iniskin sequence. C: Map view ofslab models around the 2018 Anchorage sequence. D: Cross-section showing the depthvariation of the models along I1-I2. E: Cross-section showing the depth variation of themodels along A1-A2.103A.4 Decay of the pressure in a pressurized fault zoneLet’s envision a fault zone of width 2w, which has been pressurized by frictional heating and/orwater released by chemical reactions during an earthquake. It is hypothesized that this pressurizationcauses the rapid drop in the stress ratio R = σ1−σ2σ1−σ3 , as documented in Fig. 2.12 in Chapter 2. Thestress ratio has been argued to be inversely related to the pressure. The interpretation of Fig. 2.12is that the pressure rises rapidly at the time of the earthquake and subsequently falls as a functionof time due to diffusion into the wall rock. This process is modeled as a 1D diffusion problem, asshown in Fig. A.9Figure A.9: Fault zone geometry, where w is the haft with of the fault zone, D f r is diffusionconstant of the fault zone the and Dwr is the diffusion constant of the rock wall.A.4.1 Analytic modelIngersoll et al. [1948] give an analytic solution of the thermal diffusion equation for the transientcooling of dike having the geometry in Fig. A.9, where the initial temperature of the slab−w < x <w is T degrees above that in the surrounding semi-infinite half-spaces. The diffusion equation forfluid flow in a porous medium has the same form as that for heat transferdPdt= Dd2Pdx2(A.1)where the temperature perturbation T is replaced by the pressure perturbation P and the thermaldiffusivity D is replaced by the dynamic viscosity defined asD =kφµct(m2/s). (A.2)104In equation A.2, k is permeability, φ is porosity, µ is fluid viscosity, and ct is the compressibilityof the saturated porous rock. The compressibility ct depends on the compressibility of the fluid andthat of the porous matric as ct = c f + cφ . Typical values of these constants for granite are given inTable A.1. When used in (2) they yield D = 1.7x10−6 m2/s. There is considerable uncertainty inthese values, particularly at depths on the order of 100 km depth. The resultant value of D shouldbe viewed as an order of magnitude estimate.Symbol Description Valuek permeability of granite 1.5 x 10−17 (m2)φ porosity of granite 0.004 – 0.015µ viscosity of water 3 x 10−4 (Pa·s at 100°C)c f compressibility of water 4.4 x 10−10 (Pa−1)cφ compressibility of rock matrix 10−8 - 10−9 (Pa−1 for jointed rock)ct total compressibility of rock + water (c f + cφ ) 10−9 - 10−11 (Pa−1 for solid rock)Table A.1: Value of parameters in the diffusion equation.The initial conditions are P(x) = P0 for−w≤ x≤w and P(x) = 0 for−w> x>w. The transientsolution from Ingersoll et al. [1948] isP(x, t) =P02[−er f(−w− x2√Dt)+ er f(w− x2√Dt)]. (A.3)At the center of the fault zone (x = 0),P(0, t) = P0 er f(w2√Dt). (A.4)A.4.2 Relation between R and PTo fit the data in Fig. 2.12, P(0, t) must be converted to R(0, t). Although it has been argued that R isinversely correlated with P, the functional form of this relation is uncertain. It is assumed that R andP are linearly related as R = AP+B where the constants A and B are determined by P(Rmin) = P0and P(Rmax) = 0. These constraints give B = Rmax and A = (Rmin−Rmax)/P0, soR(0, t) = Rmax− (Rmax−Rmin)er f(w2√Dt). (A.5)It is apparent that R(0, t) depends on both w and D.The linear relation between R and P that was assumed in the derivation of equation A.5 can bejustified using arguments developed in Healy [2012]. For a strike-slip fault, the symmetry axis of thefractures comprising the damage zone surrounding the core lies in the plane of the fault. The normalvectors to the fractures are all normal to the symmetry axis with a random azimuthal distribution.105This is what Healy calls a “uniaxial girdle”.If the fractures are all connected, and thus have the same fluid pressure Pf , the effective stressesare σ ′1 = σ1−Pf , σ′2 = σ2, and σ′3 = σ3−Pf so the stress ratio becomesR′=σ ′1−σ′2σ ′1−σ ′3=σ1−Pf −σ2σ1−σ3 = R−Pfσ1−σ3 (A.6)Note that R′= R when Pf = 0, and that R′= 0 when Pf = σ1−σ2. In this case, the stress ratiodecreases linearly with Pf as assumed in the derivation of equation A.5.The scenario envisaged for the Iniskin earthquake sequence is slightly more complicated thanthat outlined above. In particular, the geometrical context of the principal stresses with the respect tothe dominant fracture set where the latter is generated (the outer rise, seaward of the Alaska trench),differs from that in the region (∼100 km depth) of the Iniskin sequence. This reorganization resultsin an obliquity that variably affects the effective stresses in each of the principal axis directionsand leads generally to a non-linear dependence between R′and P. However, under conditions ofCoulomb failure, this dependence is expected to remain negative (i.e. decreasing R′manifests in-creasing P), and for the sake of simplicity it is assumed that the relation is adequately approximatedas linear.The fit shown in Fig. 2.12 of Chapter 2 gives w2√D= 2.84. Converting D = 1.7x10−6 m2/s toD= 0.147 m2/day gives a fault half-width of w= 2.2 m. The trade-off between w and D means thatan order of magnitude increase in D produces a factor of√10 ≈ 3 increase in w for the same fit.Explicit finite difference modeling was used to explore the effect of increasing the permeability ofthe fault zone rock relative to that of the wall rock. Fig. A.10 shows that making the permeabilityin the fault zone a factor of 2 larger than that of the wall rock results in a faster decay of pressureat the center of the fault zone. This means that a wider fault zone is required to fit the data for agiven wall rock permeability. The scatter in the data in Fig. 2.12 precludes an attempt to quantifythis permeability contrast.A.4.3 Relation between the Omori-Utsu law and pressureIt was hypothesized that the observed decrease in R directly following the mainshock and its subse-quent increase over the following 150 days (see Fig. 2.12) is caused by a rapid pressurization of thefault zone associate with the mainshock followed by a slow decrease in pressure as fluid diffusesaway. This hypothesis is supported by the observed temporal distribution of aftershocks.The temporal distribution of aftershocks is usually described by the Omori-Utsu law [Utsu,1961]dNdt=k(t+ c)p(A.7)106Figure A.10: Comparison of analytic and numerical solutions when the diffusion constant forthe fault zone rocks is twice that of the wall rock. Note that the larger diffusion constantin the fault zone leads to a more rapid drop in fault zone pressure. Fitting the R decaydata in this case requires a wider fault zone.where dN/dt is the rate of aftershocks (without regard to magnitude) and t is the time followingthe mainshock. The parameter k specifies the overall level of activity, c is a small time-delay thataccounts for missing events at the onset of the process, and p specifies the decay rate. For mostaftershock sequences 0.5 < p < 1.5 with an average near p = 1.Integration of equation A.4.3 givesN(t+ c) = k(t+ c)1−p (A.8)Figure A.11 shows the cumulative number of observed aftershocks and the best fit of the Omori-Utsu law to the data for 200 < t < 750 days. The sudden increase in activity at about t = 750 days107coincides with the largest event in the aftershock sequence (Mw = 5.1) and its aftershocks. Sincethis event is so large and occurs so long after the mainshock, it is not considered as an aftershock andterminate the aftershock analysis at 750 days. The discrepancy for t < 200 days can be reduced byincreasing c to a value near c = 10, but this is an unusually large delay and degrades the fit between200 and 750 days.Figure A.11: Cumulative number of aftershocks as a function of time after the mainshock.The dashed curve is the Omori-Utsu Law (equation A.4.3) with parameters k = 23.2,c = 0.01, and p = 0.74.Figure A.12 shows the difference between the data and the Omori-Utsu fit in Fig. A.11. Notethat the excess of observed aftershocks is maximum during the first 30 days following the mainshockand then decays to zero within about 150 days. This is close to the observed decay time of R inFig. 2.12. The low value of R directly following the mainshock has been ascribed to the role ofincreased fluid pressure in the activation of many faults that are not optimally oriented for Coulombslip [Healy, 2012]. As the fault-zone pressure drops, the number of misoriented faults decreasescausing the observed rise in R. Similarly, the activation of non-optimum faults directly following108the mainshock results in an excess of aftershocks above those predicted by the Omori-Utsu law thatfits the later part of the sequence where the excess pressure has decayed away.Figure A.12: Difference between the cumulative number of aftershocks and the Omori-Utsulaw (equation A.4.3) as a function of time after the mainshock from Fig. A.11. Notethat the difference is largest within about 30 days after the mainshock and decays tozero after about 150 days. This is roughly the same decay time as that observed for thestress ratio R in Fig. 2.12. The sudden increase in events at t ∼ 750 days is interpretedas being caused by a large event and its aftershocks that are not part of the aftershocksequence.109Appendix BSupporting Materials for Chapter 3B.1 Data AvailabilityThe following data are provided on Zenodo: phase catalogue (phases.txt), polarity catalogue (po-larity.txt), relocation file (relocation.txt) and moment tensor catalogue (MT solutions.txt). Fur-thermore, the synthetic and full relative moment tensor inversion codes are available on Zenodo(DOI: 10.5281/zenodo.10850729). Future development of the algorithm presented here will beconducted on GitHub as part of the Marie-Sklodowska-Curie Action funded by the European Union(https://github.com/wasjabloch/relMT).The following is a list of all seismic networks used in this study: Canadian Seismic ResearchNetwork (C8), Canadian National Seismograph Network (CN) [Natural Resources Canada, 1975],United States National Strong-Motion Network (NP) [U.S. Geological Survey, 1931], Plate Bound-ary Observatory Borehole Seismic Network (PB), Portable Observatories for Lithospheric Analy-sis and Research Investigating Seismicity (PO), USArray Transportable Array (TA) [IRIS Trans-portable Array, 2003], Pacific Northwest Seismic Network - University of Washington (UW) [Uni-versity of Washington, 1963], Monitoring Active Faults for Tectonic Mapping Efforts in Washing-ton State (X4) [Cakir, 2016] Cascadia Array of Arrays (XG) [Creager et al., 2009], Seismic ArrayStudies of Cascadia Deep Tremor (XH) [Malone, 2004], Collaborative Research: Earthscope inte-grated investigations of Cascadia subduction zone tremor, structure and process (XU) [Malone et al.,2006], Active Seismic Study of Episodic Tremor and Slip, Detecting Structural changes During anETS Event: Proof of Concept (Z3) [Silver, 2009].110B.2 FiguresFigure B.1: Same as Fig. 3.5 from Chapter 3, but two extra stations were added, so there aresix stations for DC and eight for full MTs.Figure B.2: Same as Fig. 3.5 from Chapter 3, but 30 events are now used.111Figure B.3: Earthquake locations color-coded by year before relocation. A) Top view witha black cross indicating the cluster center. B) East-west cross-section centred on thecluster. C) South-north cross-section centred on the cluster.112B.3 Moment Tensor Results for 350 EarthquakesFigure B.4: Colour code used to plot the moment tensors on Figs B.5 to B.13.In the following Figure (Figs B.5 to B.13), MTs are colour-coded with their uncertainties (seeFig. B.4 for colour bar information). MTs are in the lower hemisphere. Purple +/blue x identifythe station with P/S relative information. Circles identify stations with polarity information (emptywith pink rim: down; filled black: up), and the pattern matches the projection. Black planes are thesolution part of the uncertainties analysis. Under each mechanism is the date, event ID and numberof relative P amplitudes in purple/relative S amplitudes in blue / P polarities in gray. See othersupplements for the numerical values.113Figure B.5: Moment tensor results from 2001/11/29 to 2004/8/13.114Figure B.6: MT results from 2004/8/18 to 2007/4/16.115Figure B.7: MT results from 2007/4/17 to 2011/1/3.116Figure B.8: MT results from 2011/2/16 to 2012/7/4.117Figure B.9: MT results from 2012/7/4 to 2013/9/3.118Figure B.10: MT results from 2013/9/11 to 2017/1/30.119Figure B.11: MT results from 2017/1/30 to 2018/8/20.120Figure B.12: MT results from 2018/9/23 to 2019/7/23.121Figure B.13: MT results from 2019/8/21 to 2020/7/11.122Figure B.14: Moment tensor solutions for stable events clustered by source type using the Un-weighted Pair Group Method with arithmetic mean (UPGMA), a hierarchical clusteringmethod. Kagan angles between pair of focal mechanism as distance metric. A mini-mum of 10 events per cluster with a maximum difference of 25° Kagan are required toform a cluster. Other refer to the uncluster moment tensors.123Appendix CSupporting Materials for Chapter 4C.1 Data AvailabilityThe following is a list of all seismic networks used in this study: National Tsunami Warning CenterAlaska Seismic Network (AT) [NOAA National Oceanic and Atmospheric Administration (USA),1967], Canadian Seismic Research Network (C8) [Geological Survey of Canada, 2002], CanadianNational Seismograph Network (CN) [Natural Resources Canada, 1975], United States NationalSeismic Network (US) [Albuquerque Seismological Laboratory (ASL)/USGS, 1990].124C.2 Supplementary Figures and TablesC.2.1 Misfit for the different velocity modelsTable C.1: Each group’s residual and polarity agreement with the five different velocity mod-els used. The green text highlights the lowest relative amplitude and highest polarityagreement for each group. The green box is the final solution used for each group.Note that the choice of the velocity model is a bit arbitrary, but the results should be similar, asonly the events with a median Kagan angle difference smaller than 30°.125C.2.2 MT solutions for the different velocity modelsFigure C.1: Legend for Figs C.2 to C.5. A) The colour bar shows the agreement between thefinal solution and the other velocity model. Only events with median Kagan angle below30° (green & yellow) are used in the analysis. B) Colour code used for the nodal planes.ID corresponds to the ID defined by Oliva et al. [2024]. The filled beachball correspondsto the final solution and is filled with the associated colour from the median Kagan angle(A). All solutions are in lower hemisphere representation.126Figure C.2: MTs results for group O1. Legend shown in Fig. C.1. Note that the MT solutionfor the event #61846 was not recovered by the Bustin et al. [2007] velocity model.127Figure C.3: MTs results for group O2. Legend shown in Fig. C.1.128Figure C.4: MTs results for group O3. Legend shown in Fig. C.1.129Figure C.5: MTs results for group O3 (continuation). Legend shown in Fig. C.1.130C.3 Analysis of the earthquake group L1C.3.1 ResultsThe moment tensor inversion for L1 (see Fig. C.6 for location) did not give stable results for allvelocity models. However, the stable inversion yields reasonable results, which are presented anddiscussed in this Appendix section.Figure C.6: Overview of the earthquake group L1 location with the mapped land faults. Thepink dots are earthquakes part of L1. Same data as Fig. 4.2. CCF: Copper Creek Fault;DCF: Dawson Cove Fault; HGT: Haida Gwaii Thrust faults; MPF: MT. Poole Fault;QCF: Queen Charlotte Fault; RSFZ: Rennell Sound Fault Zone; SCF: Sachs Creek Fault.Brown lines are surface faults [Cui et al., 2017; Thompson, 1990]. The white lines is thecross-section direction used for Fig. C.8.The earthquake group L1, forming a quasi-linear structure, is particularly interesting as most of131its seismicity occurred in a four-month interval (August 2013 to mid-November 2013), ten monthsafter M7.8 2012 earthquake (see Fig. C.9) and are around 20 km deep. These events have only beendetected by Oliva et al. [2024]. Table C.2 summarizes the main information about the cluster.Group Date RangeRecovered /Total EventsDepth Range(km)Structure LenghEstimate (km)L12013-04-262019-12-03135 / 154 10 - 24 35Table C.2: General information about earthquake group L1. Recovered/Total Events refers tothe number of recovered MTs, and the total number of events that are part of the inversion,respectively.Earthquake group L1 was relocated rather than using the Oliva et al. [2024] hypocenters be-cause automatic S picks on a noisy vertical component at station MOBC (refer to Fig. C.6 forstation location) were incorrect and should have been labelled as P. Relocation was performed usinghypoDD [Waldhauser and Ellsworth, 2000] employing a 1D velocity model and both catalogue andcross-correlation delay time data. The 1D velocity model was determined by averaging velocitiesfrom the Oliva et al. [2024] 3D-velocity model within a 50 km area surrounding the geographiccenter of the group. This relocation was important as it corrected a modest in-depth curvature ofgroup L1 toward MOBC resulting from the error. The constituent earthquakes now delineate analmost vertical, planar structure (see Fig. C.7 for location comparison).Figure C.7: Relocation comparison. A) Before relocation. B) After relocation. Earthquakesare dot coloured by depth. The red triangle is the station MOBC.132Table C.3 displays the misfit for the different velocity models. Group L1 contained no eventswith known magnitudes for absolute amplitude scaling, so the magnitude of the smallest recoveredevent was fixed to 1.0 for the relative amplitude residual calculation. The 1D Oliva et al. [2024]phase average velocity model (model 5) is by far the best. The inversions using the velocity modelsfrom Bustin et al. [2007], Kao et al. [2012], and Oliva et al. [2024] group average produced verylarge relative amplitude differences, indicating significant differences in order of magnitude betweenthe maximum and minimum magnitudes. The studied events were only detected by Oliva et al.[2024] and not by NRCan, suggesting that the event local magnitudes (ML) are likely below 3.Given the sparse regional seismic network, the smallest detectable event should not be much lowerthan magnitude 1, as events smaller than this would likely go undetected. Therefore, the expectedmagnitude range should span only about two orders of magnitude. Thus, these three velocity modelssuggest instability in the inversion results, making them unreliable for consideration. Farahbod andCassidy. [2022] velocity model gives reasonable relative amplitude residual; however, its polarityagreement is low (19% below 1D Oliva et al. [2024] phase average). Hence, only the MT solutionsfrom the 1D Oliva et al. [2024] phase average velocity model are considered. Lastly, in this scenario,compared to the off-shore cluster, no large velocity contrast disturbing the averaging techniqueis expected since the ray path stays on the continental plate. Therefore, 1D Oliva et al. [2024]phase average velocity model should perform well. The results of individual earthquakes with theirestimated confidence are shown in the section C.3.4, but note that all MTs are presented here.Table C.3: Residual and polarity agreement for the five different velocity models. The greentext highlights each group’s lowest relative amplitude and highest polarity agreement.The green box is the final solution used for each group. Red values are unreasonablevalues.Group L1 forms an almost linear structure in map view and exhibits many nearly pure vertical(or horizontal) dip-slip mechanisms as shown in Figs C.8 via maps and C.10 A via the Kaveriandiagram. Due to the geometry of the relocated earthquakes, the true fault plane(s) is (are) vertical,with uplift predominantly to the northeast. The Kaverian diagram for L1, Fig C.10 A, displays alarge concentration of vertical dip-slip faults leaning toward normal fault, with a few other mecha-nisms, mainly toward the edge of the structure (see Fig. C.8). Fig. C.10 B shows the stress inversionresults for 500 noisy realizations. The principal stress axis is well-defined with small uncertainty.133The stress ratio, R, is 0.65.Figure C.8: Moment tensor solutions colour-coded by depth for group L1. Plot A displaysmap view (lower hemisphere) MTs with pink lines showing the orientation of the NW-SE cross-section in plot B. The blue arrows indicate the direction viewed in the cross-sections. Brown lines on A are surface faults from Cui et al. [2017].134If the smallest event magnitude is fixed to 1, the maximum magnitude for the inverted eventwould be 3, within the expected magnitude range for these earthquakes. The magnitude distributionversus time shows no obvious pattern (see Fig. C.9). As said previously, event activity for this groupwas concentrated within a short time interval and came in short bursts. Additionally, as reported byOliva et al. [2024], there is no apparent spatiotemporal correlation in the occurrence of the events.These behaviours are observed in earthquake swarms. Swarms have often been interpreted to becaused by fluid [Vidale and Shearer, 2006] or magmatic activity [Hill, 1977]. In the present case, afluid cause is more likely as no volcanic activity has been observed in the region recently.Figure C.9: Time evolution of the earthquakes with their estimated magnitudes. Note this plotfocuses on the most active time period and that the number of earthquakes on the right y-axis is based on all events composing L1 from 2013-04-26 to 2019-12-03. The smallestevent magnitude was fixed to 1, and the other magnitudes were scaled according to theirrelative amplitude.C.3.2 DiscussionThis earthquake group is notable for its lateral extent and depth (∼ 20 km) along a SE-NW strikingapproximately linear fault structure, with most recovered MTs exhibiting vertical dip-slip faulting.Other mechanisms are mainly located toward the edges of the structure and may be associated withstress adjustments at the fault termini. Given that the fault has an estimated length of approximately35 km and a width of around 10 km, a complete rupture of the fault could potentially generate anearthquake with a magnitude of approximately 6.5, which would cause important shaking in theregion.The seismicity of L1 can be placed within a firmer tectonic context by incorporating constraintsfrom previous geophysical studies. The significant depth of the events (many events are below20 km deep, and the deepest is close to 24 km) positions them close to the continental Moho, asestimated by Bustin et al. [2007] using receiver functions, which is dipping landward, from∼25 kmdeep west to∼28 km to the east of station MOBC. Moreover, a landward dipping low-velocity zone135Figure C.10: A) Kaverian diagrams showing the different types of earthquake mechanismsfor groups L1. This diagram was prepared using Focal Mechanisms Classification(FMC) python code from A`lvarez Go´mez [2019]. B) Stress inversion results for 500noisy realizations. The stress ratio value, R, is under the stereonet. Grey lines on thestereonet are by 15° increments.(LVZ) has been documented under MOBC at a depth of ∼ 40 km and has been interpreted to bethe subducting oceanic crust [Bustin et al., 2007; Gosselin et al., 2015]. The simple top of PacificPlate model presented by Oliva et al. [2024] (their Fig. 5 transect C−C′) locates at approximately30-35 km depth below group L1, or 6-10 km below the deepest earthquakes and is consistent withbut partly based on the receiver function estimates. According to the thermal models presented byWang et al. [2015] and based on heat flow data, seismicity in L1 occurs where temperatures approach350-400°C, near the brittle-ductile transition. Due to its young age, a mantle wedge separating thebottom of the group from the top of the Pacific plate is not expected to be significantly serpentinized.Interestingly, L1 seismicity lies parallel or sub-parallel to multiple surface faults including theSandspit Fault, the Copper Creek Fault, and the Rennell Sound Fault zone [Thompson, 1990] (seeFig. 4.3). Dyke swarms that are thought to have formed between Middle Eocene and MiddleMiocene in the vicinity display similar trends [Souther and Jessop, 1991]. Sutherland [1968] notedthat the Rennel Sound Faults Zone east block has been vertically displaced 2,000 to 10,000 feetdownward with a large horizontal offset of 12-58 miles. Thus, near-surface geology suggests a pre-existing fabric in the region favouring the generation of northwest-southeast faults. We also note,however, that L1 seismicity is oriented close to but subparallel (∼5-10° degrees counterclockwise)from the HGT and the offshore seismicity situated along the Oshawa rise and near the QCF.Perhaps the simplest explanation is that L1 seismicity results from a crustal fault activation trig-gered by stress changes induced by the M7.8 2012 earthquake. It aligns approximately with trendsdisplayed by the Sandspit Fault, Mt. Poole Fault, Sachs Creek Fault, Dawson Cove Fault, Copper136Creek Fault, and the Rennell Sound Fault zone, all of which are steeply dipping normal faults. Someexhibit a northeast downthrown side, whereas others a southwest downthrown [Thompson, 1990].The age of these faults is not well constrained, but they existed before the Terirary dyke swarmssince they would have acted as conduits for heat and plutons [Sutherland, 1968].However, the surface faults, likely only a few kilometres deep, may not reflect the deeper crustalfabric at 20 km. Additionally, this scenario does not explain why the seismicity is concentratednear the Moho/brittle-ductile-transition depths, nor does it account for the lack of a spatiotemporalcorrelation in seismic activity.A second explanation would imply a stronger link to the subduction. The vertical fault, L1,is also nearly parallel to the outer-rise seismicity (Oshawa rise) observed by Farahbod and Kao[2015] and Oliva et al. [2024], which exhibits normal faulting according to Kao et al. [2015]. Theflexure of the Pacific Plate, as it bends to subduct under the North American Plate, causes theseouter-rise faults. Outer-rise hydrates the slab by allowing fluid penetration into the crust and, insome cases, in the mantle before it subducts [Peacock, 2001]. Two outer-rise aftershocks from theM7.8 2012 have M>6 with a nominal fault dimension is ∼10km, suggesting that water penetrationreaches the Pacific mantle. Previous outer-rise faults probably have similar water depth penetration.The captured fluid can later be released through dehydration reactions [Peacock, 2001] or pressure-temperature changes, which force fluid migration.The stress change induced by the M7.8 2012 earthquake could have triggered fluid release froman outer-rise fault. This fluid may have migrated upward in the mantle and upon reaching the crusta few months later. The build-up of fluid pressure at the base of the Moho could have causedhydraulic fracturing and fluid migration into a pre-existing fault or geologic interface. Hydraulicfracturing occurs when the pressure gradient exceeds the rock strength, leading to the creation orexpansion of fractures. In this scenario, the fluid was generated by the oceanic plate, where hydro-static pressure (the pressure exerted by the fluid) and lithostatic pressure (the pressure exerted bythe weight of overlying rocks) are in equilibrium. As the fluid moves upward, the rate of changeof lithostatic and hydrostatic pressure with depth differs, generating a pressure gradient that couldexceed the rock’s strength and result in fracturing. Note that natural hydraulic fracturing does notproduce large earthquakes, but when occurring on existing faults, they can be detected, supportingthat those earthquakes have small magnitudes. Furthermore, hydraulic fracturing is likely to openvertical faults, the orientation with the largest pressure gradient. Lastly, the lack of spatio-temporalcorrelation supports this explanation.C.3.3 ConclusionEven though the results from the different velocity models did not all converge, the MTs resultsobtained using Oliva et al. [2024] phase average seem reasonable. The ∼35 km long inland lin-ear structure at depths between 10 and 24 km shows predominantly vertical faulting with similar137orientation to surface faults, the Haida Gwaii Thrust and out-rise faults. Additionally, this crustalstructure terminates at the continental Moho. The cause of those earthquakes could be explained byfluids reactivating a pre-existing fault.C.3.4 Moment Tensors solutionsHere are all MTs results for L1. The final solution is the coloured solution. The colour, as scaledin Fig. C.11, is based on the median Kagan angle difference between the final solution and the bestsubset of solutions. The best subset of solution was selected as described in Chapter 3 section 3.3.3Step 8 - Stage 4: Quality Control.Figure C.11: Colour-code legend for the best solution.In the following Figures (Figs C.12 to C.15) the MTs are in the lower hemisphere. Purple +/ blue x identify the station with P/S relative information. Circles identify stations with polarityinformation (empty with pink rim: down; filled black: up), and the pattern matches the projection.Black planes are the solution part of the confidence analysis. Under each mechanism is the date,event ID and number of relative P amplitudes in purple / relative S amplitudes in blue / P polaritiesin grey.138Figure C.12: MTs results from 2013/4/26 to 2013/8/23.139Figure C.13: MTs results from 2013/8/26 to 2013/9/13.140Figure C.14: MTs results from 2013/9/14 to 2013/10/30.141Figure C.15: MTs results from 2013/11/6 to 2013/11/19.142
Thesis/Dissertation
eng
Geophysics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Bostock, Michael G.
Advancements in seismotectonic analysis through relative moment tensor inversion
Text