Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Seismic velocity structure under Vancouver Island from travel time inversion : insight from Low Frequency… Savard, Geneviève 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2018_september_savard_genevieve.pdf [ 69.61MB ]
JSON: 24-1.0371609.json
JSON-LD: 24-1.0371609-ld.json
RDF/XML (Pretty): 24-1.0371609-rdf.xml
RDF/JSON: 24-1.0371609-rdf.json
Turtle: 24-1.0371609-turtle.txt
N-Triples: 24-1.0371609-rdf-ntriples.txt
Original Record: 24-1.0371609-source.json
Full Text

Full Text

Seismic velocity structure under Vancouver Island from travel timeinversion: Insight from Low Frequency EarthquakesbyGeneviève SavardB.Sc., Physics, Université de Montréal, 2012a 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)August 2018© Geneviève Savard, 2018The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled:Seismic velocity structure under Vancouver Island from travel time inversion: insight from Low Frequency Earthquakessubmitted by Geneviève Savard in partial fulfillment of the requirements forthe degree of Doctor of Philosophyin GeophysicsExamining Committee:Michael G. Bostock, Earth, Ocean and Atmospheric SciencesSupervisor Andrew J. Calvert, Earth Sciences, Simon Fraser UniversitySupervisory Committee Member Catherine L. Johnson, Earth, Ocean and Atmospheric SciencesUniversity ExaminerSimon M. Peacock, Earth, Ocean and Atmospheric SciencesUniversity ExaminerAdditional Supervisory Committee Members:Nikolas I. Christensen,  Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberMark Jellinek, Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberiiAbstractThe Cascadia subduction zone is a warm subduction zone where Episodic Tremor and Slow Slip(ETS), a phenomenon with debated governing physical mechanisms, is observed. The eventualdevelopment of an accurate physical model of ETS relies on a comprehensive set of seismic andgeodetic observations. Low Frequency Earthquakes (LFE) are one of the seismic signaturesthat accompany slow slip and possess discernible phase arrivals that can be used for traveltimeinversion methods. I first address the issue of detecting LFEs by developing an automatedalgorithm that exploits cross-station waveform similarity. The algorithm detects thousands ofevents per ETS episode in southern Vancouver Island that reveal pronounced spatiotemporalclustering and complex propagation patterns that compare favourably with independent ob-servations in Cascadia. This method allows us to build a significantly improved catalogue ofLFE templates that we next employed for local double-difference earthquake tomography withthe aim of resolving structures in the vicinity of ETS and metamorphic controls on its genera-tion. In southern Vancouver Island, tomographic images reveal high Poisson’s ratios associatedwith a dipping low‐velocity zone (LVZ) inferred to be overpressured, upper oceanic crust ofthe Juan de Fuca plate where LFEs and other slow‐slip phenomena occur. We also observe alow Poisson’s ratio anomaly (∼0.225) in the forearc continental crust above the mantle wedgecoinciding with high Vp and high levels of clustered microseismicity. We develop a conceptualmodel where quartz concentration is produced by metasomatism in the forearc crust and cat-alyzed by a focussed ingress of slab-derived fluids at high pore pressure. In northern VancouverIsland, we also observe the LVZ and a low Poisson’s ratio anomaly in a similar position, but noassociated microseismic activity. Distinct, parallel lineaments that trend slightly oblique to thestrike of the megathrust are resolved within seismicity concentrations of the Nootka Fault Zoneimmediately offshore Nootka Island. We speculate that these features may represent faultingdeveloped as a result of high strains associated with slab bending, amplified in the vicinity ofthe newly formed plate boundary. We also find possible evidence for a deep-seated, seismogenicfault near the Explorer Plate edge.iiiLay SummaryThe time of arrival of a seismic wave generated by an earthquake at a network of seismographsat the surface depends upon the location and origin time of the earthquake and the structure ofthe Earth through which it travels. In this thesis, I develop a detection algorithm that can beused to improve the clarity and accuracy of arrival time observations for a recently discoveredtype of earthquake called Low Frequency Earthquake (LFE). Using the arrival times of LFEsand other sources of seismicity, I am able to calculate improved estimates of location for theseevents along with the structure of the Earth under the tectonically active region of VancouverIsland. The models I developed provide insights into the seismotectonics of Vancouver Island.ivPrefaceThis thesis presents research work carried out in its entirety at the University of BritishColumbia by me with the guidance, advice and editorial authorship of my supervisor, MichaelBostock. The three principal research projects presented in this thesis were initially identifiedand suggested to me by Michael. My supervisory committee was composed of Mark Jellinek(UBC), Andrew Calvert (Simon Fraser University) and Nikolas Christensen (UBC, Emeritusprofessor). All members of the committee provided constructive comments and suggestions dur-ing committee meetings and in personal discussions. Nik Christensen was particularly involvedin the interpretation of the tomographic results presented in Chapter 3 and is a co-author ofthis publication.Chapter 2 has been published in Bulletin of the Seismological Society of America as:Savard, G., and Bostock, M. G. (2015). Detection and Location of Low-Frequency EarthquakesUsing Cross-Station Correlation. Bulletin of the Seismological Society of America, 105(4). is reproduced unmodified except for formatting changes. The electronic supplement tothis article is included in Appendix A. I developed the algorithm with guidance from MichaelBostock. I was the lead author and responsible for the analysis of the data and composition ofthe manuscript. Michael Bostock provided feedback, advice and editorial supervision through-out the development of this paper. The results of this research were also presented at differentinternational conferences.Chapter 3 is an accepted manuscript in Geochemistry, Geophysics, Geosystems and cur-rently available as a preprint online:Savard, G. , Bostock, M. G. and Christensen, N. I. (2018), Seismicity, Metamorphism, and FluidEvolution Across the Northern Cascadia Forearc. Geochem. Geophys. Geosyst.. AcceptedAuthor Manuscript. doi:10.1029/2017GC007417It is reproduced unmodified from its preprint version except for formatting changes. TheSupplementary Material of this paper is included in Appendix B. I am the lead author and wasresponsible for the construction of the dataset and tomographic inversions. Takeshi Akuharaprovided help with the use of the double-difference tomography code. Michael Bostock andvNikolas Christensen were highly involved in the interpretation of the models. The analysis ofthe results and proposed conceptual model is a team effort and the result of multiple discussionsbetween myself, Michael Bostock and Nikolas Christensen. We also received critical inputand benefited from helpful discussion from Simon Peacock, Ken Hickey, Andrew Calvert andAlexandre Plourde.Chapter 4 is a manuscript currently in preparation for future publication. Part of theLow Frequency Earthquake template catalogue is the result of unpublished work by LindsayYuling Chuang as a graduate student at UBC. Michael Bostock provided help with the tedioustask of culling false detections and obtained the STA-LTA detections. I completed all arrivaltime picks, constructed the dataset and completed the traveltime inversions. Michael Bostockprovided advice, help with interpretation and editorial supervision of the Chapter.Funding for my research came from a UBC Affiliated Fellowship, NSERC Alexander GrahamBell Canada graduate scholarship (Master’s), NSERC postgraduate scholarship (Doctoral) andUBC Four Year Fellowship awarded to me as well as an NSERC Discovery Grant awarded tomy advisor Michael Bostock.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Supplementary Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Seismicity along the Cascadia subduction zone . . . . . . . . . . . . . . . . . . . 11.2 Fundamentals of traveltime inversion . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.1 Seismic body waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.2 The earthquake location problem . . . . . . . . . . . . . . . . . . . . . . . 51.2.3 The double-difference technique . . . . . . . . . . . . . . . . . . . . . . . . 61.2.4 Basics of local earthquake traveltime tomography . . . . . . . . . . . . . . 71.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Detection and Location of Low Frequency Earthquakes Using Cross-station Correlation 92.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.1 Search Window Parameters and Station Network . . . . . . . . . . . . . . 112.2.2 Search Window Admission/Rejection . . . . . . . . . . . . . . . . . . . . . 122.2.3 LFE Detection via Multi-channel Cross-correlation . . . . . . . . . . . . . 142.2.4 S-wave Arrival Time Estimation and Epicenter Location . . . . . . . . . . 162.2.5 Hypocenter Determination . . . . . . . . . . . . . . . . . . . . . . . . . . 16vii2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3.1 Epicenters for Southern Vancouver Island and Comparison with OtherCatalogues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3.2 Spatio-temporal Clustering of LFE Epicenters . . . . . . . . . . . . . . . 192.3.3 Hypocenters and Depth Distribution . . . . . . . . . . . . . . . . . . . . . 262.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313 Seismicity, Metamorphism, and Fluid Evolution Across the Northern Cascadia Forearc 343.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.2 Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.4.1 Slab Fluid Evolution From Trench to Tremor . . . . . . . . . . . . . . . . 433.4.2 Fluid Expulsion Into Forearc Crust . . . . . . . . . . . . . . . . . . . . . . 473.4.3 Forearc Crustal Composition . . . . . . . . . . . . . . . . . . . . . . . . . 483.4.4 Quartz Concentration Mechanisms . . . . . . . . . . . . . . . . . . . . . . 483.4.5 Crustal Forearcs and Orogenic Gold Deposits . . . . . . . . . . . . . . . . 513.4.6 Leech River Fault and Texada Island Seismicity . . . . . . . . . . . . . . . 523.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534 Double-difference Tomography of Northern Vancouver Island . . . . . . . . . . . . . . 554.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.2 Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.2.1 Data collection and selection . . . . . . . . . . . . . . . . . . . . . . . . . 644.2.2 Minimum 1-D Vp model inversion . . . . . . . . . . . . . . . . . . . . . . 674.2.3 3D graded inversion using TomoDD . . . . . . . . . . . . . . . . . . . . . 734.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.3.1 Checkerboard resolution tests . . . . . . . . . . . . . . . . . . . . . . . . . 744.3.2 Relocated Seismicity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.3.3 Vp, Vs and Poisson’s ratio models . . . . . . . . . . . . . . . . . . . . . . 864.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 945 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.1 Detection and Location of Low Frequency Earthquakes Using Cross-station Cor-relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.2 Seismicity, Metamorphism, and Fluid Evolution Across the Northern CascadiaForearc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.3 Double-difference Tomography of Northern Vancouver Island . . . . . . . . . . . 97Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99viiiA Supplementary material for Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 113A.1 Effect of Station Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113A.2 Epicenter Locations for Northern Washington . . . . . . . . . . . . . . . . . . . . 114A.3 Time-distance Plots for the 2003 and 2005 Episodes . . . . . . . . . . . . . . . . 115B Supplementary material for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 120ixList of FiguresFigure 1.1 Simplified cartoon of the Cascadia Subduction Zone. . . . . . . . . . . . . . . 2Figure 1.2 Simplified cartoon of Episodic Tremor and Slow Slip. . . . . . . . . . . . . . . 3Figure 2.1 Cross-station detection work flow for LFE epicenter and hypocenter detection. 12Figure 2.2 Station density and maximum S-wave delay times for 2004. . . . . . . . . . . 13Figure 2.3 Particle-velocity waveforms and envelope functions for an example of detectedLFE epicenter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Figure 2.4 Particle-velocity waveforms and envelope functions for an example of detectedLFE hypocenter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18Figure 2.5 Detected LFE epicenters in southern Vancouver Island for the ETS episodesof 2003, 2004 and 2005. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21Figure 2.6 Spatio-temporal evolution of LFE epicenters between January 2003 and June2006. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22Figure 2.7 LFE epicenter density for the ETS episode of 2003, 2004 and 2005 and theinter-ETS periods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 2.8 Spatio-temporal evolution of LFE epicenters for the ETS episode of 2004. . . 25Figure 2.9 Examples of LFE migration patterns in southern Vancouver Island. . . . . . 26Figure 2.10 Detected LFE hypocenters for southern Vancouver Island. . . . . . . . . . . . 28Figure 2.11 Example of detection from our catalogue that is used as an initial templateto produce a high SNR template using the network cross-correlation andstacking procedure of Bostock et al. (2012) . . . . . . . . . . . . . . . . . . . 30Figure 2.12 Final template obtained from the initial template of figure 2.11. . . . . . . . . 33Figure 3.1 Large-scale geologic map of study region. . . . . . . . . . . . . . . . . . . . . 36Figure 3.2 Model Poisson’s ratio versus depth along the sections AA’, BB’ and CC’shown in figure 3.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 3.3 Model Vp along the sections AA’, BB’, and CC’ shown in Figure 3.4. . . . . 41Figure 3.4 Relocated epicenters after inversion of the full travel time dataset. . . . . . . 42Figure 3.5 Earthquake clusters linked by high waveform cross-correlation coefficients. . . 44xFigure 3.6 Vp model profiles and earthquake relocations along the eastern section of theLeech River Fault and western end of the Darrington-Devils Mountain FaultZone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 3.7 Schematic model of fluid evolution A) across northern Cascadia forearc, andB) within dehydrating oceanic crust. . . . . . . . . . . . . . . . . . . . . . . . 46Figure 3.8 Plot of forearc crustal anomaly Poisson’s ratio and Vp (red ellipse) on corre-sponding laboratory measurements for various crustal lithologies as measuredby Christensen (1996). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 4.1 Tectonic context of the study region. . . . . . . . . . . . . . . . . . . . . . . . 56Figure 4.2 Regional seismo-tectonics of northern Vancouver Island. . . . . . . . . . . . . 57Figure 4.3 Regional faults and major terranes of the study region. . . . . . . . . . . . . 59Figure 4.4 Major geological groups of northern Vancouver Island. . . . . . . . . . . . . . 60Figure 4.5 Relief and station map for northern Vancouver Island. . . . . . . . . . . . . . 62Figure 4.6 Event catalogue obtained from GSC and LFE hypocenters. . . . . . . . . . . 65Figure 4.7 Minimum 1D Vp model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 4.8 Station corrections computed by VELEST for the minimum 1D model. . . . 71Figure 4.9 Tomographic dataset and relocations within the minimum 1D model. . . . . 72Figure 4.10 RMS traveltime residual distributions after relocation within minimum 1Dmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.11 Map view of earthquakes and LFEs relocated with tomoDD. . . . . . . . . . 75Figure 4.12 Depth slices at 10 and 20 km depth for the checkerboard test with anomaliessized 40x40x20 km3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 4.13 Depth slices at 30 and 40 km depth for the checkerboard test with anomaliessized 40x40x20 km3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Figure 4.14 Depth slices at 5 and 10 km depth for the checkerboard test with anomaliessized 20x20x10 km3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78Figure 4.15 Depth slices at 15 and 20 km depth for the checkerboard test with anomaliessized 20x20x10 km3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 4.16 Depth slices at 25 and 30 km depth for the checkerboard test with anomaliessized 20x20x10 km3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80Figure 4.17 Depth slices at 35 (top) and 40 (bottom row) km depth for the checkerboardtest with anomalies sized 20x20x10 km3. . . . . . . . . . . . . . . . . . . . . . 81Figure 4.18 Along dip depth profiles at y=-40,-20,0,20,40 km for the checkerboard testwith anomalies sized 40x40x20 km3. . . . . . . . . . . . . . . . . . . . . . . . 82Figure 4.19 Along dip depth profiles at y=-40,-20,0,20,40 km for the checkerboard testwith anomalies sized 20x20x10 km3. . . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.20 Along strike depth profiles at x=-80,-60,-40,-20,0 km for the checkerboardtest with anomalies sized 40x40x20 km3. . . . . . . . . . . . . . . . . . . . . . 84xiFigure 4.21 Along strike depth profiles at x=-80,-60,-40,-20,0 km for the checkerboardtest with anomalies sized 20x20x10 km3. . . . . . . . . . . . . . . . . . . . . . 85Figure 4.22 Along dip depth profiles at y=-65,-50,-35 km, near the Nootka Fault Zone. . 89Figure 4.23 Same as 4.22, but for along dip depth profiles at y=-20,-5,10 km, near LFEseismicity and the middle of the Explorer plate. . . . . . . . . . . . . . . . . . 90Figure 4.24 Same as 4.22, but for along dip depth profiles at y=25,40,55 km, near theBrooks Peninsula Fault Zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 4.25 Same as 4.22, but for along strike depth profiles at x=-60,-45,-30,-15,0 km. . 92Figure 4.26 Close-up of the area around the Holberg Fault. . . . . . . . . . . . . . . . . . 93Figure 4.27 Velocity structure near the Holberg Fault. . . . . . . . . . . . . . . . . . . . . 94Figure A.1 Station density and maximum S-wave delay times for 2003. . . . . . . . . . . 113Figure A.2 Station density and maximum S-wave delay times for 2005. . . . . . . . . . . 114Figure A.3 Station density and maximum S-wave delay times for 2013. . . . . . . . . . . 114Figure A.4 LFE epicenters for the 2010 and 2011 ETS episodes in northern Washington. 116Figure A.5 LFE epicenter density for the ETS episodes of 2010 and 2011 in northernWashington. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117Figure A.6 Spatio-temporal evolution of LFE epicenters for the 2003 ETS episode. . . . 118Figure A.7 Spatio-temporal evolution of LFE epicenters for the 2005 ETS episode. . . . 119Figure B.1 Vs model for profiles AA’, BB’ and CC’. . . . . . . . . . . . . . . . . . . . . . 121Figure B.2 Synthetic checkerboard model. Note that profiles cross checkerboard gridobliquely leading to some fading for certain sections. . . . . . . . . . . . . . . 122Figure B.3 Recovered checkerboard model. Relocated seismicity in the checkerboardmodel is only shown to illustrate data coverage. . . . . . . . . . . . . . . . . . 123Figure B.4 Synthetic high Vp/Vs LVZ model. . . . . . . . . . . . . . . . . . . . . . . . . 124Figure B.5 Recovered high Vp/Vs dipping layer. Relocated seismicity in the syntheticmodel is shown to illustrate data coverage. . . . . . . . . . . . . . . . . . . . 125Figure B.6 Average Poisson’s ratio and P-wave velocity along profiles AA’, BB’ and CC’. 126xiiList of Supplementary MaterialsResults from Chapters 3 and 4 are made available on cIRcle. This online repository containsthe 3D velocity models for Vp and Vs and the hypocenter relocations from the travel timetomography studies presented in this thesis.The folder “SVI” contains data for southern Vancouver Island and “NVI”, for northernVancouver Island. Each folder contains 7 files:• tomoDDrelocations.csv: There are 12 columns: event ID, year, month, day, hour, minutes,seconds, latitude, longitude, depth (km), cluster (optional), type. The first column, “eventID” contains individual event identification numbers that we assigned to each event andcan be ignored. “Type” is either “catalog” for events present in the GSC earthquakecatalog (or, for some events for SVI, obtained from the previous study of Thomas Brocher),“lfe” for our LFE catalog and “detection” for events we detected with our modified cross-station method and used for tomography. For the SVI dataset, “cluster” is between 1and 27 for the clusters identified through cross-correlation coefficients linkage (and henceconcerns events between 2002 and 2006 for which we had waveform data available). Thevalue is 0 otherwise.• grdinfo.dat: Information on the origin of the coordinate system (X-Y-Z) and inversiongrid nodes in km from the origin.• Vpmodel.dat: Vp value for each X-Y-Z point of the coordinate system in the followingstructure: 266666666664Vp(1, 1, 1) Vp(2, 1, 1) . . . Vp(nx, 1, 1)Vp(1, 2, 1) Vp(2, 2, 1) . . . Vp(nx, 2, 1)... ... . . . ...Vp(1, ny, 1) Vp(2, ny, 1) . . . Vp(nx, ny, 1)... ... . . . ...Vp(1, ny, nz) Vp(2, ny, nz) . . . Vp(nx, ny, nx)377777777775(1)where nx, ny and nz are the number of grid nodes in the X, Y and Z direction respectively.• Vsmodel.dat: same as Vpmodel.dat but for S-waves.xiii• inversiongrid.dat: corresponding latitude and longitude for each line of the velocity modelfile, separated by “>”.• P.dws: Derivative Weight Sum value for each data point of the Vp model.• S.dws: same as P.dws but for S-waves.Additional data requests should be sent by email to Geneviève Savard at:, I would like to acknowledge the people who kindly provided the datasets that werecritical to this research. I thank Kumar Ramachandran and Thomas M. Brocher for sharingtheir SHIPS travel time catalogue for P-wave first arrivals and metadata. I thank Cliff Thurberfor his permission to use the TomoDD code, Geological Survey Canada staff (namely AllisonBird and Taimi Mulder) for providing me with their traveltime catalogue between 1992 and2014 and Honn Kao for sharing the SEAJADE dataset and his tremor catalogue for southernVancouver Island. I also thank my co-author Nik Christensen for his geological expertise (andthis sample of eclogite I will always keep as a precious souvenir for this Ph.D.!), and mycolleagues Alexandra Royer, Takeshi Akuhara, Lindsay Yuling Chuang and Alexandre Plourdefor their professional help. Thank you to the members of my supervisory committee, Andrew,Nik, and Mark, for their advice and helpful discussion.Five years five months ago, as I was boarding a plane for Vancouver, I was excited at theprospect of a new intellectual challenge and of living close to the Coast Mountains, but alsoterrified at the thought of being completely immersed in an English environment. After all, thelast time I had given an oral presentation in English, I had passed out in front of the class fromnervousness. Would I be able to read lips in English? Would I end up completely isolated dueto communication struggles? What was I getting myself into? My dad wisely said, worst casescenario, I would just have to come back home, but I could say that I tried. Thankfully thisdidn’t happen, and I owe a great deal to the wonderful people I met during the course of thisPh.D., for their presence during the good and the hard moments.First, Michael, thank you for being a fantastic supervisor, for being so welcoming, for yourcontagious enthusiasm for our research and optimism, for your understanding and for keeping meon track throughout this Ph.D. I am also immensely grateful for the opportunities I had to sharemy work and travel to conferences in Japan, Alaska, San Francisco and Mexico. These wereprecious experiences to me, professionally and personally. Thank you for your encouragementsand help particularly during these last few months.I thank Alexandra Royer, Ben Postlethwaite, Lindsay Chuang, Alexandre Plourde and TimHayward for being more than just fellow seismology grad students, but also true friends tome. I also thank the medley of amazing people, the “MJCJ crew”, that I shared an office withthrough the years, there are too many to list, but I’d like to give special thanks to Kathi Unglertand Anna Mittelholz for always lending an ear and providing advice in the tough moments.xvI thank the people of Accurate Realtime who provided real-time captioning during mycourses and group meetings, allowing me to not get lost in fast-paced group discussions. Itwas of tremendous help to me and I am amazed by your skills that I consider akin to virtuosopianists.Life away from home wasn’t too bad thanks to all these weekend adventures in the mountainsthat I’ve shared with wonderful people. Thank you to everyone with whom I shared a rope onthe cliffs, a ski run in amazing powder (or slushy coastal cement), a kayak trip amongst theorcas, a slog up scree slopes, glaciers or logging roads, a delicious backcountry meal, a dip inglacial lakes, a triumphant summit moment or epic mountain views.Thanks to my father Daniel, my mother Veronique and my sister Marie-Eve for your constantpresence despite the distance between us. And finally, thanks to everyone who encouraged meto embark on this adventure and board that plane in December 2012. It was one of the bestdecisions I ever made.xviDedicationJe dédie cette thèse à mon père, Daniel Savard.xviiChapter 1Introduction1.1 Seismicity along the Cascadia subduction zoneIn order to build accurate earthquake hazard models for a given region, a comprehensive por-trayal of all relevant tectonic structures and associated seismogenic potential needs to be made.The southwestern coast of Canada is located along the Cascadia Subduction Zone, one of manyactive subduction zones located around the Pacific Ring of Fire. As shown in figure 1.1, itis formed by the subduction of the young (<12 Ma) and warm oceanic Juan de Fuca plateunder the continental North America plate in a northeastward direction at an average rate of40 mm/yr (e.g., Wilson, 2002). There are three main types of earthquakes with destructionpotential occurring in this region: megathrust (or interplate) earthquakes, deep intraslab (orWadati-Benioff) earthquakes and crustal earthquakes. Megathrust earthquakes are producedby the sudden release of accumulated stress on a locked area of the plate interface. Becausesubduction zone megathrusts are the largest structures susceptible to seismic slip, they producethe world’s largest earthquakes, such as the 1960 Mw 9.5 Great Chilean earthquake, 2004 Mw9.1-9.3 Sumatra-Andaman earthquake and the 2011 Mw 9.1 Tohoku earthquake. Relative toother subduction zones, the Cascadia megathrust is characterized by infrequent but large rup-tures with an average recurrence time of 500-600 years (e.g., Goldfinger et al., 2013). Only afew rare earthquakes of small or intermediate magnitude (Tréhu et al., 2012) has been detectedin the locked seismogenic zone of the megathrust and the latest significant seismic release ofaccumulated plate convergence stress occurred during the great earthquake of 1700 (e.g., At-water et al., 2005), indirectly inferred to have been of magnitude Mw 8.7-9.2 (e.g., Hyndmanand Wang, 1995).11.1. Seismicity along the Cascadia subduction zoneFigure 1.1: Simplified cartoon of the Cascadia Subduction Zone. Source: US GeologicalSurvey.(Retrieved on May 31 2018 from (or “slab”) earthquakes are located within the subducting oceanic crustand/or mantle, caused by the sudden seismic rupture of pre-existing normal faults and otherstructural weaknesses. Although their position within the slab mutes their effect on surfaceground motion, they present a significant hazard due to their inland epicentral location undermajor cities. The 2001 Mw 6.8 Nisqually earthquake is the latest example of a large destructiveintraslab earthquake in the region (Kao et al., 2008). The distribution of intraslab earthquakesin combination with thermo-petrologic models of slab structure provide insight into the meta-morphic reactions occurring within the slab as it encounters higher temperature and pressureregimes. However, the physical mechanisms leading to deep intraslab earthquakes are still de-bated (e.g., Preston et al., 2003; Hasegawa and Nakajima, 2017). Intraslab earthquakes can beclassified as intermediate or deep focus earthquakes, depending on if they occur shallower ordeeper than 300 km respectively. Cascadia has no recorded deep focus earthquakes.Crustal earthquakes are located within the North America plate along existing, tectonicallyactive faults such as the Seattle Fault. In northern Cascadia, the large-scale stress regime,obtained from moment tensor solutions of regional crustal earthquakes (Balfour et al., 2011),is dominantly margin parallel, indicating N-S oriented compression and caused by the north-ward migration of the forearc Oregon Block (e.g., Balfour et al., 2011). Slip along crustalfaults accommodate this large scale stress accumulation. The seismic hazard caused by theseearthquakes is significant due to the proximity of large crustal faults near population centres.A fourth type of seismicity named tectonic tremor was first discovered in 2002 (Obara, 2002)in southwest Japan and accompanies the periodic phenomenon of slow slip along the transitionzone between megathrust rupture and stable sliding (Rogers and Dragert, 2003). Tremor is anoisy, semi-continuous signal recorded in the 1-10 Hz frequency band that lacks a clear phaseonset, but is inferred to occur in the same area as slow slip as it migrates slowly over a period ofdays to weeks along-strike. Figure 1.2 presents a graphical representation of tremor epicenters21.1. Seismicity along the Cascadia subduction zoneoccurring below northern Washington during an Episodic and Slow Slip (ETS) episode. Shellyet al. (2007) have shown that tremor is produced by swarm of low frequency earthquakes (orLFEs), that share the same bandwith as tremor but diplay discrete P and S arrivals. Becauseof the low signal-to-noise ratio of individual LFE waveforms, stacking techniques are needed tobuild an “LFE template” that reveals clear P, S and converted phase onsets (e.g., Bostock etal., 2012).Figure 1.2: Simplified cartoon of Episodic Tremor and Slow Slip. Source: IRIS, with cour-tesy of K. Creager.(Retrieved on May 31 2018 from focal mechanisms of LFE templates are consistent with shear slip on the plate interface(e.g., Royer and Bostock, 2013). The physical mechanisms responsible for the generation ofLFEs, their limited bandwith and distribution are still debated (e.g., Beroza and Ide, 2011).One reason is that this relatively newly discovered phenomenon is challenging to detect in-strumentally and requires advanced signal processing techniques and dense seismic networks tocharacterize. Hence, a complete and consistent portrayal of the locations, spatio-temporal be-haviour and seismic characteristics of tremor and LFEs is still in development. Theoretical andconceptual modelling of ETS depends upon a comprehensive set of observations. Consequently,there is a need for comprehensive catalogues of tremor epicenter and LFE hypocenters, whichis one of the issues addressed in this thesis.31.2. Fundamentals of traveltime inversion1.2 Fundamentals of traveltime inversion1.2.1 Seismic body wavesAn earthquake can be represented as a point dislocation source. The solution to the elasticdisplacement field generated by such a source in an homogeneous, isotropic and unboundedmedium generates two types of body waves: compressional (P or Primary) waves with onelongitudinal polarization direction and shear (S or Secondary) waves, with two transverse po-larization directions (SH in the horizontal plane and SV in the vertical plane). In the far-field,high-frequency approximation, body wave propagation in heterogeneous media can be stud-ied using ray theory: the propagation of a P or S body wave is represented by a ray pathwith a local speed Vp or Vs, an amplitude and direction determined by geometrical spreadingand Snell’s law, respectively. This thesis exploits traveltime inversion techniques where thetraveltime T is defined as the time a seismic wavefront takes to travel along a ray path fromearthquake source to seismograph. Traveltime is computed through solution of the eikonal equa-tion (∇T)2 = 1v2 = u2 where v is the local wavespeed (Vp or Vs) and u is the local slowness,the inverse of the local velocity.The velocities Vp and Vs depend upon the elastic properties and density of the materialthey travel through and are given by:Vp =√l+ 2mr=√k + 43mr(1.1)Vs =√mr(1.2)where k is the bulk (or compression) modulus, m the shear modulus, l the first Lamé parameterand r the density. Poisson’s ratio s is a related elastic parameter that is particularly importantas a diagnostic for different rock types. For an isotropic material subjected to uniaxial stress,it is defined as the negative ratio of transverse strain to axial strain. Poisson’s ratio can bedetermined from measurements of Vp and Vs, assuming an isotropic medium, as:s =12(1− 1(Vp/Vs)2 − 1)(1.3)Rocks that compose the Earth’s crust and mantle are made of mineral aggregates. The averageVp, Vs and s of a rock depend upon the volumetric proportion of each constituent mineral.Elastic velocities also vary as a function of porosity, fluid pressure, temperature and pressure.Poisson’s ratio has the advantage that it is largely independent of pressure and temperaturethereby enhancing its value for interpretation purposes. Poisson’s ratio for most common crustallithologies varies between 0.25 and Fundamentals of traveltime inversion1.2.2 The earthquake location problemAn earthquake hypocenter yields essential information on its source. For small rupture areas,the earthquake source can be approximated as a point source and defined by four parameters:three spatial coordinates (latitude, longitude and depth contained within a location vectorx⃗) representing the location from which seismic energy was radiated, and an origin time trepresenting the onset of seismic radiation from this source. These four parameters need to beinferred indirectly from the measurement of m seismic phase arrivals at a network of seismographstations. Given a velocity model of the Earth and a ray tracing algorithm with which we cancompute predicted arrival times, a minimum of 4 arrival times is needed from a given phaseto obtain an hypocenter solution. However in practice, the earthquake location problem iscomplicated by:• accurate identification of the onset of seismic phases in a station’s record. S-waves andother reflected and refracted phases always appear within the P-wave coda, and theironset is more challenging to identify.• the choice of an appropriate velocity model accurately representing Earth’s structurethrough which the seismic waves travel and the accuracy of the ray tracing algorithmused to compute theoretical ray travel times.• the geometry of the stations at which arrivals were picked, which can greatly affect theaccuracy and reliability of the hypocenter estimate. The better the azimuthal distributionaround the earthquake source, the better will be the estimate.The goal is to recover the hypocenter m = (⃗x, t) that best fits the m arrival time measurementstk given a chosen reference Earth model and a ray tracing algorithm that allow forward mod-elling of predicted arrival times tPk = F(m). Because the operator F that computes predictedarrivals depends on m itself and because we have typically m > 4 , earthquake location is anoverdetermined nonlinear inverse problem. Most routine earthquake location algorithms use theiterative least squares method of Geiger (1912) that linearizes the inversion problem by lookingat a first order perturbation of the hypocentral parameters dm around an initial guess m0 closeto the true solution. The vector of residuals rk = tk − tPk is written as r = Gdm, when G isthe m× 4 Jacobian matrix of partial derivatives ¶tPk¶mj for the kth observation and jth hypocentralparameter . We invert for dm that minimizes the L2 norm of r −Gdm. This least squaresinversion is iterated using a new initial guess m0 + dm at each iteration until a residual cutoffcriterion is reached. The system of normal equations is commonly solved using the singularvalue decomposition method (SVD), that allows calculation of uncertainties via the covariancematrix of the generalized inverse of G. The least squares solution is statistically the solutionof maximum likelihood in the presence of normally distributed noise. However, one drawbackof this method is that it is dependent upon a good initial guess close to the true value and isnot guaranteed to reach a global minimum in the misfit function.51.2. Fundamentals of traveltime inversionSince the effect of a perturbation in the origin time and of a perturbation in depth aresimilar, there is a strong trade-off between the depth estimate and the origin time that is bestresolved by using both P-wave and S-wave arrivals. In the absence of sufficient P and S arrivaltimes, the earthquake’s depth is commonly fixed to a designated value and one can invert forthe epicenter alone (latitude, longitude and origin time).1.2.3 The double-difference techniqueThe double-difference location algorithm was developed by Waldhauser and Ellsworth (2000)and is nowadays widely used to obtain highly precise relative earthquake hypocenters, usefulfor delineating the presence of faults and other seismogenic structures. The double-differencealgorithm minimizes the effect of unmodelled 3D velocity structure without the use of stationcorrections and has been shown to improve uncertainties in relative hypocenter locations bymore than an order of magnitude.It builds upon the result of Fréchet (1985) that relates the relative hypocentral parametersbetween two events i and j to their differential traveltime residual for an observation k: ¶tijk¶mdmij =drijk , where drijk is the so-called double difference, the residual between observed and calculateddifferential travel for two events drijk = (tik − tjk)obs − (tik − tjk)calc. The observed differentialtravel time may be obtained via waveform cross-correlation or by subtracting absolute catalogarrival times. Generalizing this equation for the case where the slowness vector is not assumedconstant between the two events gives:¶tik¶midmi − ¶tjk¶mjdmj = drijk (1.4)Using this equation for a set of M double-difference observations for N events yields the sparsesystem of equations:W[GlI]m = W[d0](1.5)whereW represents the a priori quality weights for each observation andG is the M× 4N matrixof partial derivatives, d is the vector of M double difference observations, m is the 4×N vectorof hypocentral parameters and l is the damping factor. For a large number of earthquakes,this problem is solved as a regularized nonlinear least squares optimization problem using theconjugate gradient algorithm LSQR of Paige and Saunders (1982).The inversion is iterated until a stopping criterion is reached and the data are re-weighted ateach-iteration using biweight and biexponential functions to control the effect of large residualsand large inter-event distances, respectively.61.2. Fundamentals of traveltime inversion1.2.4 Basics of local earthquake traveltime tomographyPioneered in 1976 by Aki and Lee (Aki and Lee, 1976), the local earthquake traveltime to-mography problem is a generalization of the earthquake location problem to the simultaneousinversion of earthquake hypocenters and 3D velocity structure in the volume containing theearthquakes and stations. A standard way to parametrize 3D velocity structure is by dividingthe 3D volume in L blocks of constant velocity. The model parameter vector to solve for nowincludes L velocity parameters in addition to 4× n hypocentral parameters of n earthquakes.Given m observations and n earthquakes, there are m · n equations to solve. For an observationk from event j:rjk = dtj +¶Tjk¶xdxj +¶Tjk¶ydyj +¶Tjk¶zdzj +Lål=1¶Tjk¶vldvl (1.6)Because of the last term, there is coupling between the velocity structure and hypocenter param-eters that introduces significant nonlinearity. Hence one commonly completes many iterationsthat alternate between joint hypocenter-velocity inversions and hypocenter relocation inver-sions to minimize this effect. Due to the linearization of the tomographic nonlinear problem,the choice of an appropriate 1D, 2D or 3D starting model is often critical to achieve convergenceto a global minimum in the misfit function.The ability to resolve sharp, spatially localized features in velocity structure depends stronglyon the model parametrization (such as velocity block size), the quality of arrival time picks, thesmoothing constraints and on the sampling of anomalies by rays with a wide range of incidentangles. Because of the latter requirement, the distribution of seismicity and stations is criticalto achieve a desired tomographic resolution. Moreover, because the shallow structure of theEarth immediately beneath the stations has poor ray directionality coverage and is usuallycharacterized by strong, vertical velocity gradients, station correction terms commonly needto be introduced. The forward problem of ray tracing within complex 3D structure, which isnecessary to compute predicted traveltimes, is a complex computational problem in itself thatrepresents the biggest part of the computational burden of most tomographic methods. As aresult, a wide variety of algorithms for 3D ray tracing have been proposed, all with individuallimitations and pitfalls.As for the location problem, the system can be written in matrix form as r = Gdm.Regularization is again necessary due to ill-conditioning. Commonly it penalizes the roughnessof the model, and can be solved via the damped least squares optimization problem:min‖Gm− d‖22 + l2‖Lm‖22 (1.7)where L is an operator measuring the roughness of the model such as the Laplacian ∇2, and lis the smoothing parameter that balances data misfit and model roughness.In this thesis, we employ the double-difference tomography technique of Zhang and Thurber(2003). It extends the double-difference location algorithm described above to the seismic71.3. Thesis outlinetomography problem. Since the double-difference measurements also contain information on thenear source velocity between pairs of events, it has proven to significantly improve resolution inthe source region of clustered earthquakes. The use of differential traveltime measurements alsoobviates the need for station corrections. Forward modelling of the ray paths is accomplishedusing the approximate pseudo-pending algorithm of Um and Thurber (1987).1.3 Thesis outlineThe following chapters address three main research objectives.In Chapter 2, I address the difficulties in identifying the initial LFE waveforms necessaryto obtain high signal-to-noise ratio (SNR) LFE templates such as those obtained by Bostocket al. (2012) using network cross-correlation and stacking techniques. I develop an efficienton-the-fly detection algorithm for LFEs that exploits cross-station waveform similarity in orderto find new source locations and expand the spatial extent of the LFE template catalogue inthis region.Since high SNR LFE templates present clear P and S arrival onsets, they can be usedfor traveltime earthquake tomography. Because LFEs are situated in a region of very lowbackground seismicity in Cascadia, they present a unique research opportunity to extend raycoverage in the vicinity of the ETS region. Their inclusion within tomographic inversions of theregion holds the potential to provide additional insight into the location of the plate interface,the position of the mantle wedge corner and the metamorphic reactions that produce slabfluids and the generation of ETS. In Chapter 3, I use the improved LFE catalogue to invertfor 3D velocity structure under southern Vancouver Island via double-difference tomography.The choice of double-difference tomography is motivated by the fact that it provides improvedresolution near the source region, the subducting oceanic crust in the case of LFEs. In additionto the LFE catalogue, I include regular crustal and intraslab earthquakes and active sourcedata. My results show that LFEs provide a unique insight into the velocity structure near andabove the ETS source region, which in turn allows us to make inferences on the role of fluidsand metamorphic reactions taking place near the plate boundary and into the overlying forearc.Encouraged by the success of our study for southern Vancouver Island, I develop a double-difference tomographic model for northern Vancouver Island, where new LFE templates werealso identified using our cross-station detection method. Northern Vancouver Island is of par-ticular interest due to its unique tectonic setting, an end-member of warm subduction environ-ments, and because no previous 3D tomographic study has been performed in the region. I alsoidentify the effect of major tectonic structures such as the Nootka Fault Zone and the positionof the plate edge on the forearc and near plate boundary velocity structure. Our results arepresented in Chapter 4.Finally, in Chapter 5, I conclude this thesis by identifying my main research contributionsand some future research directions that naturally follow from my work.8Chapter 2Detection and Location of Low FrequencyEarthquakes Using Cross-station Correlation2.1 IntroductionDuring the last decade, it has been established that the transition between locked and smoothlysliding portions of some subduction plate boundaries is the site of a complex, periodic stressrelease that manifests seismically as tremor in the 1-10 Hz frequency band and geodetically asseaward motion of the overriding plate suggesting slow shear slip at or near the plate interface.Every 8 to 16 months, tectonic tremor and slow slip occur between 25 and 40 km depth indifferent along-strike segments of Cascadia (Brudzinski and Allen, 2007) and last for periodsof several days to weeks. This phenomenon, designated episodic tremor and slow slip (ETS),may harbour the potential to improve our understanding of megathrust seismogenesis. Someauthors have suggested that the slow slip episodes in the upper part of the transition zonecould transfer stress to the seismogenic zone or shift its lower limit landward (Dragert et al.,2004; Wech and Creager, 2011). Chapman et al. (2009) have derived a seismogenic couplingprofile from ETS data for the Washington segment of Cascadia and concluded that the lowerlimit for future megathrust rupture could extend to a depth of 25 km, or 10 km deeper thanpreviously thought (Hyndman and Wang, 1993). During ETS, it has been reported that tremorepicenters evolve as primary and secondary fronts showing complex migration patterns (Houstonet al., 2011; Ghosh et al., 2010; Rubin and Armbruster, 2013), notably sudden changes in themain front along-strike propagation direction, first identified and termed rapid tremor reversals(RTR) by Houston et al. (2011).Low frequency earthquakes (LFEs) are small (M ∼ 1-2) events sharing the frequency spec-trum of tectonic tremor. The seismic signal is noisy, but often shows discernable S-wave, andsometimes P-wave, arrivals. Shelly and coworkers have demonstrated that tremor consists ofswarms of LFEs and confirmed that LFE focal mechanisms correspond to shear slip parallel onor near the plate interface and parallel to the plate motion (Shelly et al., 2007; Ide et al., 2007).92.2. MethodShelly et al. (2006) used LFE detections, manually picked by network analysts of the JapanMeteorological Agency, as templates for waveform matched filtering on seismograms from thewestern Shikoku region, to reveal highly repetitive sources distributed along the top of thesubducting Philippine Sea plate. The temporal record of repetition of an LFE template hasunveiled tremor modulation by tides (Shelly et al., 2006; Thomas et al., 2012) and triggeringby transient stresses caused by surface waves generated by large earthquakes (Miyazawa andBrodsky, 2008; Rubinstein et al., 2007; Peng and Chao, 2008). It also facilitates the assem-bly of LFE Green’s functions for structural studies using high SNR templates formed throughstacking of multiple LFE detections (Bostock et al., 2012; Nowack and Bostock, 2013; Matharuet al., 2014).The observation that LFEs occur repeatedly at localized asperities with the same sourcemechanism has led to multi-channel network detection approaches to search for waveform rep-etitions and automatically identify LFEs for use as new templates (Aguiar and Beroza, 2014;Brown et al., 2008). A second class of approach, more akin to the manual/visual approachtaken by analysts and termed “cross-station” approaches by Rubin and Armbruster (2013),relies on the similarity of waveforms from an individual LFE recorded across multiple stations.Although such waveform similarity might not be anticipated a priori, Royer and Bostock (2013)found that it clearly holds for S-waves generated by LFEs in Cascadia for near-vertical sourcegeometries owing to simple Green’s functions and the dominantly shallow thrust mechanism ofthe LFE source. Armbruster, Rubin and co-workers (Armbruster et al., 2014; Rubin and Arm-bruster, 2013) have exploited this similarity to examine details of slow slip process in southernVancouver Island. The method of Frank and Shapiro (2014) also falls into this class of approachbut the location method is cast in terms of a grid search rather than linearized inversion.Motivated by these findings, we employ cross-station waveform similarity in the developmentof a technique to detect individual, impulsive LFEs that can be used to assemble high SNRLFE templates using iterative network cross-correlation. An advantage of this approach is itsefficiency and its utility in the selection of prospective templates from targetted epicentral sites,in contrast to e.g. network autocorrelation where location is conducted only after expensivepattern recognition approaches are applied.In the following sections, we will first explain our method, then we present sample resultsfor southern Vancouver Island and northern Washington and compare with existing tremorcatalogues obtained using other methods (Wech and Creager, 2008; Wech, 2010; Kao et al.,2006).2.2 MethodOur cross-station algorithm shares some elements with that described by Armbruster, Rubinand co-workers (Armbruster et al., 2014; Rubin and Armbruster, 2013) but is tailored specif-ically to the detection and location of discrete, impulsive LFEs using a full complement ofstations available for a given region and period. The primary focus is the detection and epicen-102.2. Methodtral location of LFEs using S-waves with a subsequent analysis devoted to hypocentral locationusing P-waves where these can be confidently identified. The S-wave detection strategy com-prises the following steps: i) selection of optimal search window parameters based on networkconfiguration, ii) window admission/rejection based on spectral SNR, iii) LFE detection basedon multichannel cross-correlation (MCCC), iv) S-wave arrival time estimation and location.These steps are organized in the flowchart in figure 2.1 and are detailed in the following sub-sections using examples from the POLARIS array on southern Vancouver Island (Nicholson etal., 2005).2.2.1 Search Window Parameters and Station NetworkThe cross-station method we present here places specific demands on both station density andnetwork geometry with respect to a given source region. Because post-critical, free surfaceinteractions can distort waveforms (Booth and Crampin, 1985), stations should ideally be po-sitioned near the LFE epicentral region in order to maximize likelihood of waveform similarity.The pre-critical regime ends where the ray parameter exceeds the P-wave slowness at the sur-face, that is ∼0.2 s/km (or equivalently, an S-wave incidence angle exceeding ∼35.3◦ for the1D velocity model used by the Geological Survey Canada for routine earthquake location innorthern Cascadia.Figure 2.2a shows a map of stations available for 2004 on southern Vancouver Island (Nichol-son et al., 2005) along with contours of the number of stations with pre-critical arrivals availablefor a grid of hypocenters located on the plate interface model of Audet et al. (2010). This plotreveals how the network geometry affects the detection of LFEs as a function of location since aminimum of three S-wave arrivals is needed to determine an accurate epicenter and origin time.It should be noted that this map is a general indicator of network coverage and does not takeinto account the specific period of availability of each station. Nevertheless, it highlights thoseregions where we should expect enhanced detection of LFEs. For 2003 and 2005, the stationdistribution is slightly different, as shown in figures A.1 and A.2 of appendix A.In our initial detection procedure we employ S waves only and search across horizontalcomponent time series at different stations to search for similar waveforms associated with anindividual LFE. In order to distinguish waveforms from distinct LFEs occurring closely in time,we need to scan the time series with a window length that is as small as possible. However,there is a trade-off between the duration of LFE S-waveforms (typically 4 s or less at stations onVancouver Island) and their moveout as determined by the network station configuration. Asa guide to selection of an appropriate window length for a region of interest, one can estimatethe maximum travel time difference among all pairs of stations in the network. A contourmap of this quantity is displayed in figure 2.2b for the same grid of hypocenters as figure 2.2a.For southern Vancouver Island, we choose the scanning window length to be 18 s with 17.5 soverlap, as a good compromise between these 2 considerations. A 24 hour record thus comprisesa total of 172,764 18-s windows.112.2. MethodIt should be noted that prior to 2003 and after 2005, the network geometry is less suited toour approach (see figure A.3 for 2013 in Appendix A). Stations are farther apart and the delaytime range between stations is large. Hence the search window must be longer, increasing thelikelihood of including multiple LFE events and complicating the analysis of cross-correlationfunctions. In this work, we restrict consideration to southern Vancouver Island data from 2003to 2005 only.See figure 2.3aSee figure 2.3bSee figure 2.3cDe!ne absolute S arrival times as t0 +δt1+δt2 and use a bogus P time to locate using HYPOINVERSE if less than 4 stations, go to next windowSet the reference absolute time t0 to the maximum of the stacked en velopes and take a window of ± 4 s around t0 for the wave formsCompute optimal dela y times δt1 for the en velopes using  multichannel cross-correlation (MCCC) and shift theenvelopes and wave forms by δt1Find the largest subset of traces with delay time uncertainties σi under the tolerance σmax. Shift according to corresponding delay times δt2.    Compute cross-correlation coe#cients Φij of the envelopes and selectstations corresponding with  Φij > C envChoice of window parameters Choice of station listExtract waveforms for all stations in search windowKeep traces with su#cient spectral SNRusing thresholds T 1 and T 2Use the epicenters and a starting depth of 32 km to  compute the predicted S and P arrival times at all stations.Extract the S waveforms with a window of ± 1.5 s around the Sarrival times for the stations used in epicenter location and both  horizontal channels. For all stations in the network, e xtract long windows on the vertical  channel starting from -4 s before the predicted P arrival time andending to the S arrival time.A P-wave arrival is registered if the normalized cross- correlation coe#cients with each S waveform have a mean  superior to 0.5 and the delay time of the P window relative to   the aligned S windows has an uncertainty smaller than 0.1 s.  LFE epicentersKeep w ave forms with timing uncertainties smaller than 0.1 s to  use as templates for P detection and align using MCCC.Use previously selected set of S waveforms as matched !lters andcross-correlate with the long P windows. Look for both a positive ornegative maximum to account for polarity uncertainty.Add newly found P times to previous S times to locate.LFE epicenters LFE hypocenterstraces S 1traces S 3traces S 2Figure 2.1: Cross-station detection work flow for LFE epicenter and hypocenter detection.2.2.2 Search Window Admission/RejectionParticle velocity waveforms for the 2 horizontal components (N,E) are extracted day by dayand the records are scanned according to the window parameters and station selection specifiedin the previous section. The N and E components are treated as independent traces since LFEwaveform similarity varies on a channel-by-channel basis, as noted for templates on in southernVancouver Island (Royer and Bostock, 2013). Armbruster, Rubin and co-workers (Armbrusteret al., 2014; Rubin and Armbruster, 2013) advocate a grid search through rotation angle todetermine a single, optimum horizontal component waveform as an additional consistency mea-sure. However, this approach becomes computationally expensive with a full complement (20+)of stations and so we forego this step. The raw seismograms are band-pass filtered between 1and 8 Hz, resampled at 40 samples-per-second and corrected for instrument response. For eachstation, we subtract the mean and normalize to the maximum absolute amplitude of the twohorizontal channels since we are interested in similarity of the waveform shape across stations,independent of amplitude.Tremor and LFEs on Vancouver Island are especially rich in frequencies in the 2-5 Hz band122.2. Method 20’  124 oW  40’  20’  123oW   48oN  12’  24’  36’  48’   49 oN 0 20 40 kmPFB TWGB TWBBTSJBLCBCLZB TWKBMGCBKLNBSSIBSNB GOWBSHDBPGC KHVBSILBSHVBVGZ   Number of stations2 4 6 8 10 12 14 16 124 oW  123 oW   49 oN 0 20 40 km  Max S−wave delay time [s]0 2 4 6 8 10 12 14 16 18 2020 km 20 km30 km 30 km40 km 40 km 12’  24’  36’  48’  20’  40’  20’ Figure 2.2: a) Map of the number of stations with pre-critical arrivals for hypocenterslocated on the plate interface model of Audet et al. (2010) for southern VancouverIsland. The 20, 30 and 40 km depth contours correspond to the dashed lines fromleft to right respectively. The inverted triangles indicate the stations used in thisstudy. LFE detections are more likely in regions of high station density statisfyingthe pre-critical arrival criterion. b) Same as a) but for the maximum delay betweenS-wave arrival times among all pairs of stations.and depleted in higher frequencies (Beroza and Ide, 2011). This property can be exploited toeliminate traces with poor spectral SNR using two thresholds T1 and T2. We compute the powerspectral density of each trace and eliminate the traces with more than T1 % of their total energyat the frequencies lower than 2 Hz or higher than 5 Hz. In addition, since the power spectraldensity of noisy traces is often dominated by a single frequency peak, we further exclude tracesusing a second threshold T2 on the percentage of the total energy contained within ± 0.5 Hz ofthe dominant frequency. The thresholds T1 and T2 may depend on the spectral characteristicsof the tremor data set under consideration. To determine appropriate values, we comparedthe spectra of traces with high noise levels (including instrumentally generated transients) withthose of traces with discernable LFEs and found that values in the range 50-90% for both ofthe two thresholds characterize low SNR traces. For southern Vancouver Island, we chose thevalue of 70% for both T1 and T2 as an acceptable trade-off between excluding excessively noisytraces but retaining traces that show potential for LFE detection despite a lower SNR. Thesesimple spectral constraints may not be sufficient to distinguish tremor from some intraplateearthquakes. However, in southern Vancouver Island the GSC catalogue recorded less than2 events per day during the ETS episodes studied, whereas LFE template catalogues showhundreds to thousands of detections per day. Hence, contamination of our catalogue by regularearthquakes is negligible.132.2. MethodAt the end of this stage, we have for each search window a set of LFE candidate traces S1 tobe further analysed. An example extracted from southern Vancouver Island data is illustratedin figure 2.3a.Figure 2.3: a) Horizontal particle-velocity waveforms showing impulsive S arrivals from anLFE on July 18 2004 (N component on the left panel, and E component on theright). Waveforms are band-passed filtered between 1 and 8 Hz, and sampled at 40sps. Asterisks indicate the traces selected for the next step. b) Envelope functions ofwaveforms in a) corresponding to envelope cross-correlation coefficients above Cenvand aligned using MCCC. The vertical dotted line indicates the absolute referencetime t0 corresponding to the maximum of the stack and the 2 continuous lines definethe time window used to mask the waveforms for subsequent fine-scale alignment. c)Subselected traces from figure b) aligned with delay times dt2 and timing uncertaintiessi < smax computed by MCCC. This is the final selection of traces used for epicenterlocation.2.2.3 LFE Detection via Multi-channel Cross-correlationFor selection S1, we seek an optimal alignment of the LFE waveforms across a minimum of4 different stations. To achieve this objective, we exploit the multi-channel cross-correlationalgorithm of Vandecar and Crosson (1990) hereafter referred to as MCCC. MCCC was originallydesigned for semi-automatic phase picking of teleseismic body waveforms and was shown toyield relative times of quality comparable or better than manual phase picking. The optimaltime delays between a set of waveforms are found through least-squares inversion of preliminarycross-correlation functions. The inversion minimizes the errors due to distorted cross-correlationfunctions caused by multi-pathing, provided that this effect can be considered as an uncorrelatednoise. We employ this algorithm within our LFE detection/location procedure because we142.2. Methodare searching for similarity across waveforms and because it provides objective, quantitativeestimates of timing uncertainties. Following Vandecar and Crosson (1990), for n traces, wecan write n (n-1)/2 overdetermined equations for the relative delay times and add an arbitraryzero-mean constraint.ti − tj = ∆tij, (2.1)nåi=1ti = 0, (2.2)where ∆tij is the lag of the maximum of the cross-correlation function Fi,j for each pair of traces,i.e. Fi,j(∆tij) = maxFi,j, and ti and tj are the arrival times at stations i and j, respectively. Sincestation separation distances are known and only pre-critical angles of incidence are allowed, weconstrain the search for cross-correlation maximum (∆tij) to a specific lag window that honoursthe criterion of pre-critical incidence for each station pair. If the channel pair corresponds tothe two horizontal components from the same station, we use a lag value threshold of 0.3 sthereby accommodating the potential for S-splitting (Bostock and Christensen, 2012; Matharuet al., 2014).The timing uncertainty si for station i is given by equation (2.3) with res denoting the residualas defined in equation (2.4).si =vuut 1n− 2(i−1åj=1res2ij +nåj=i+1res2ij), (2.3)resij = ∆tij − (ti − tj). (2.4)si can be used to discriminate consistent LFE waveforms (low si) from noisy amplitude bursts(high si). We monitor both Fi,j and si to use as selection/detection criteria.Optimal Lag Time EstimationWe employ the MCCC algorithm outlined above in 2 steps: first, to compute coarse delay timesto approximately align waveform envelopes, and second, to perform a fine-scale alignment usingwaveforms.Cross-correlation coefficients are calculated for each pair of trace envelopes in S1. To targetenergetic, impulsive wavelet envelopes, we retain those traces with coefficients above a thresholdCenv, that we set typically within the range 0.4-0.5 depending on the noise level. Alignedenvelopes are stacked and the absolute time corresponding to the stack maximum, upon whichS arrivals of the candidate LFE should be aligned, is assigned as a reference time t0. We narrowour time window to ± 4 seconds around t0 as shown in figure 2.3b by applying a temporal maskto the waveforms of the selected set S2.In order to keep only traces with internally consistent times, i.e. with a small timing152.2. Methoduncertainty si, we find the largest subset S3 of waveforms that yields a prescribed maximumtiming uncertainty, smax by computing cross-correlation coefficients for all pairs of waveformsin S2. We force the search for the cross-correlation maximum to lags within ± 2 s only,since waveforms should be roughly aligned through the envelope correlation stage. Due to thecommon shallow thrust mechanism, we expect that S waves recorded on horizontal componentchannels should exhibit consistent polarity at small epicentral distances. We set a thresholdCmin to a token maximum value (say Cmin = 1.0) and enter a loop in which Cmin is iterativelydecreased by a constant decrement. At each iteration, we admit waveforms correspondingto cross-correlation coefficients larger than Cmin, thereby increasing the number of waveformswithin S3, and compute the optimal delay times and their timing uncertainty using MCCC. Ifthe maximum uncertainty is smaller than the tolerance value smax (again set to 0.3 s to accountfor possible splitting), we continue to decrease Cmin to test additional waveforms, otherwise webreak the loop with a final set S3, which meet both the Cmin and smax criteria. This procedureensures that our final selection of stations and channels S3 and their computed time delaysdisplay similar waveforms and internal consistency to within ±smax.2.2.4 S-wave Arrival Time Estimation and Epicenter LocationFigure 2.3c shows the final selection S3 of channels for the example, aligned using the envelopedelay times dt1 and waveform delay times dt2. We evaluate the absolute S wave arrival timefor each station i of the final selection of admitted waveforms as tSi = t0 + dt1i + dt2i. Ifboth horizontal channels are represented for a given station, we adopt the mean value as thestation time. Epicentral location is performed using HYPOINVERSE, a standard linearizedinversion program (Klein, 2002), and the 1D GSC velocity model used for routine locations.We fix hypocentral depth at 32 km which is approximately the mean depth determined for LFEtemplates on southern Vancouver Island (Bostock et al., 2012; Royer and Bostock, 2013). Eachstation arrival time is weighted proportionally to its timing uncertainty si for dt2 as computedby MCCC. The epicenters with a final nominal horizontal error of less than 3 km, as computedin HYPOINVERSE, are logged as detections.2.2.5 Hypocenter DeterminationFor each LFE epicenter, we attempt to find at least one P-wave arrival on vertical componentchannels of any of the available stations in the network. We employ the most coherent S-waveforms as templates to search the vertical channel for matching waveform prior to theexpected or observed S arrival, as detailed below.First, for each station that passed the epicenter detection criterion, we extract with short(-1 s to +3.5 s) windows the S waveforms from both horizontal components. In order to retainonly S waveforms with high SNR as templates, timing uncertainties are again estimated withMCCC, and we apply a more stringent threshold of 0.1 s on them. If less than 4 waveformssatisfy this criterion, we skip the event.162.2. MethodNext, predicted P and S arrival times are computed using the GSC velocity model for afixed hypocentral depth of 32 km. For each station, we cross-correlate the short S windows withthe vertical channel inside a conservative window from -4 s before the predicted P time to theS time on the vertical component. We find the optimal time shift between the P long windowand the selected S short windows using MCCC. In contrast to S-waves, the shallow thrustmechanism produces a null in the radiation pattern of P waves that results in low amplitudesand variable polarity in the immediate epicentral region. Hence, we must search for either astrong positive or negative correlation coefficient in the normalized cross-correlation functions.A P-wave is registered if three conditions are met: 1) The cross-correlation coefficients betweenthe P window and each of the selected S windows have a mean larger than 0.4, 2) the P windowhas a correlation coefficient larger than 0.5 with at least one S window, and 3) the timinguncertainties registered by MCCC are less than 0.1 s for all traces.Delay times for the S and P windows are inverted for hypocentral location using HY-POINVERSE, and we retain hypocenters with final nominal horizontal and vertical errors ofless than 3 km. Figure 2.4 shows an example of a positive P detection. The 3 channels arealigned according to their S arrival times. The solid lines indicate the S and P times input tothe hypocenter inversion whereas the dashed lines demark the predicted arrival times for thehypocentral solution.172.2.MethodFigure 2.4: Example of a hypocenter detection. Waveforms are aligned according to their predicted S arrival times. Left and middlepanels are the N and E horizontal channels and right panel is the vertical channel. Solid lines indicate the times used forhypocenter inversion, and dotted lines indicate predicted arrival times.182.3. Results2.3 ResultsWe applied our cross-station approach to the ETS episodes of March 2003, July 2004 andSeptember 2005 on southern Vancouver Island since these events were sampled by stations of thePOLARIS-BC array (Nicholson et al., 2005) in addition to permanent stations of the CanadianNational Seismograph Network(CNSN). These three events had durations of approximately 15,16 and 20 days (Kao et al., 2009) for which we obtained about 3162, 8600 and 8413 epicentersrespectively. We also analysed data from the Washington state CAFE and Array of Arrays(Ghosh et al., 2012) experiments from ETS episodes in 2010 and 2011. The results are availablein Appendix A.2.3.1 Epicenters for Southern Vancouver Island and Comparison with OtherCataloguesEpicenters for LFEs below southern Vancouver Island are compared with the tremor catalogueobtained by the source-scanning algorithm of Kao et al. (2006) in figure 2.5. The migration ofthe slow slip main front from the southeast to the northwest is evident in both catalogues. Ourepicenters do not extend as far to the northwest as the tremor locations because the numberof stations registering arrivals at pre-critical incidence falls off in this direction, as reflected infigure 2.2. However, our epicenters extend farther downdip, along the east coast of VancouverIsland and south of Victoria (about longitude -123◦20’). We also applied the cross-stationapproach to inter-ETS periods between January 1 2003 and June 30 2006 and obtained 3300epicenters which are compared with Kao’s tremor catalogue in figure 2.6. Except for the periodof March to May 2005, we find increased LFE activity for each period of enhanced tremorbetween ETS episodes for approximately the same along-strike and along-dip distances, but wealso resolve more inter-ETS LFE “bursts”.Our LFE detection approach is not expected to replicate the template catalogue, becausenetwork cross-correlation approaches can detect LFEs at much lower SNR levels than ourmethod allows. Nevertheless, high SNR LFEs should result in some overlap between the 2catalogues. We found 3% of our epicenters have at least 3 arrivals within 0.5 s of templatedetection arrivals. For these common detections, the difference in epicenter location shows adistribution centered at 4 km falling steeply to 15 km. This is because the HYPOINVERSEerror calculation erroneously assumes that arrivals were picked with a precision comparable tomanually picked phases for regular earthquakes. Since our absolute arrivals are determined byleast squares minimization, the absolute location error can be larger due to outliers caused bycycle-skipping and low SNR.2.3.2 Spatio-temporal Clustering of LFE EpicentersThe LFE detections display strong spatio-temporal clustering. The spatial densities of epicen-ters were estimated using the kernel density estimator of Botev et al. (2010). For southernVancouver Island, the density for the three ETS episodes and the inter-ETS periods are shown192.3. Resultsin figure 2.7. For the 3 ETS episodes the spatial concentrations are similar with the exceptionof the large western “asperity” A centered at (-124◦, 48◦30’) that is largely absent for 2003.There are fewer tremor epicenters from Kao’s catalogue in this area in 2003 than in 2004 and2005. However, we have verified that the lack of detections could also be due to the absenceof stations LCBC, GLBC JRBC along the south coast and stations TSJB and LZB near thecenter of the POLARIS line during 2003. When these stations are removed from considerationin the 2004 and 2005 episodes, the western asperity A is no longer as evident, thereby un-derlining the importance of network configuration and sampling in defining spatial resolution.The inter-ETS activity is largely distinct from that determined during the major ETS episodes,with two asperities, D and E, located close to the 40 km depth contour of the plate interface,and approximately aligned in the plate motion direction with the asperities rupturing updipduring regular ETS episodes (A, B and C). These two downdip asperities are temporally andspatially distributed in a manner consistent with the reported “minor tremor” by Armbrusteret al. (2014). Figure 2.6 displays the along strike location of detections as a function of timefrom January 2003 to June 2006 and reveals a temporal clustering of inter-ETS activity. Thetwo downdip asperities rupture in alternation several times between ETS episodes.202.3. Results 20’  124oW  40’  20’  123oW  10’  20’  30’  40’  50’   49oN 20 km 20’  124oW  40’  20’  123oW 20 km  Feb−26 Feb−27 Mar−01 Mar−02 Mar−04 Mar−05 Mar−07 Mar−08 20’  124oW  40’  20’  123oW  10’  20’  30’  40’  50’   49oN 20 km 20’  124oW  40’  20’  123oW 20 km  Jul−11 Jul−13 Jul−14 Jul−16 Jul−17 Jul−19 Jul−21 Jul−22 Jul−24 20’  124oW  40’  20’  123oW  10’  20’  30’  40’  50’   49oN 20 km 20’  124oW  40’  20’  123oW 20 km  Sep−06 Sep−08 Sep−10 Sep−12 Sep−14 Sep−16 Sep−18 Sep−20 Sep−22Figure 2.5: LFE epicenters for the ETS episode of March 2003, July 2004 and Septem-ber 2005 with an horizontal error of less than 3 km. Left panels show the LFEdetections of the present study, and right panels show tremor detections from thesource-scanning algorithm of Kao et al. (2006). Triangles indicate stations used forthe detections in left panels. The shade of each detection dot is coded accordingto the time elapsed since the beginning of the episode. The along-strike migrationof tremor and LFE activity from the southeast to the northwest is evident in bothcatalogues. 212.3.ResultsDistance along strike [km]Distance along strike [km]Distance along dip [km]Figure 2.6: a) Distance along-strike of the LFE epicenters as a function of time from January 2003 to June 2006. The colorbar isproportional to along-dip distance: dark shades are updip, clear shades are downdip. The reference lines used to computealong strike and along-dip distances are found in figure 2.10. The continuous line is the cumulative number of detections. b)Same as a) but for the tremor catalogue of Kao et al. (2006).222.3. ResultsWithin an ETS episode, there is also strong temporal clustering. Our detections for the2004 episode are compared with template detections (Bostock et al., 2012; Royer and Bostock,2013) in figure 2.8. Comparisons for the 2003 and 2005 ETS episodes are shown in figures A.6and A.7 available in Appendix A. We observe that tremor migrations patterns are consistentbetween the 2 catalogues, but the cross-station approach extends the range in both the along-dip and along-strike directions. In addition to the main front migration, smaller migrations onshorter time scales can be observed. Four migrations are identified in black outlines in figure2.8 and mapped in space in figure 2.9. Events C and D (∼10 km/h) correspond to the RTRsthat were identified from the template detections by Royer et al. (2015). Events A (∼10-15km/h) and B (∼10-20 km/h) were not individually identified within the template detections. Aappears to be an RTR, with a small component of velocity in the downdip direction whereas Bsuggests a migration path proceeding from updip to downdip. Since the cross-station detectionsare not spatially limited by template availability, they provide a better indication of true RTRdimensions which appear to be of the order of 50 km in length and somewhat narrower in width.232.3.Results 20’  124 oW  40’  20’  123 oW  20’  30’  40’  50’   49 oN 0 20 40 km  0246810121416 20’  124 oW  40’  20’  123 oW  20’  30’  40’  50’   49 oN 0 20 40 km  012345678910 20’  124 oW  40’  20’  123 oW  20’  30’  40’  50’   49 oN 0 20 40 km  24681012 20’  124 oW  40’  20’  123 oW  20’  30’  40’  50’   49 oN 0 20 40 km  05101520253030 km20 km40 kma) b)c) d)ACEDBRelative densityRelative densitySalt Spring Isl.San Juan Isl.Figure 2.7: Epicenter density maps of southern Vancouver Island computed using the kernel density estimator of Botev et al. (2010)for a) 2003 ETS episode, b) 2004 ETS episode, c) 2005 ETS episode and d) inter-ETS detections from January 2003 to June2006. Triangles are the stations used to locate the events for each respective period. The dashed closed contour lines indicateapproximate outlines of the 5 main ”asperities” we qualitatively identify. Dotted lines are the 20, 30 and 40 km depth contoursof the plate interface model of Audet et al. (2010).242.3.Results020406080100120140160 a)Event count−60−40−200204060b)Distance along strike [km]Jul−12 Jul−14 Jul−16 Jul−18 Jul−20 Jul−22 Jul−24−60−40−200204060c)Distance along strike [km]  Distance along dip [km]−50 −40 −30 −20 −10 0 10 20 30 40RTRRTRRTRRTRRTR RTRRTRFigure 2.8: Cross-station detections compared with template detections for the 2004 ETS episode. a) Number of cross-stationdetections with respect to time. b) Distance along-strike is plotted against time for our cross-station detections, and datapoint shades indicate distance along dip: dark shades are updip, clear shades are downdip. The profiles used to computealong-strike and along-dip distances are shown in figure 2.10. c) Same as b) but for template detections. RTRs identified byRoyer and co-workers (Royer and Bostock, 2013; Royer et al., 2015) are shown in black outlines.252.3. ResultsA20’ 124oW 40’ 20’ 123oW 20’ 30’ 40’ 50’ 49oN 0 20 40 km05−Sep−1409:3910:2511:1111:5612:42B20’ 124oW 40’ 20’ 123oW 20’ 30’ 40’ 50’ 49oN 0 20 40 km05−Sep−1623:2400:4602:0903:3204:55C20’ 124oW 40’ 20’ 123oW 20’ 30’ 40’ 50’ 49oN 0 20 40 km04−Jul−1810:2211:4313:0414:2615:47D20’ 124oW 40’ 20’ 123oW 20’ 30’ 40’ 50’ 49oN 0 20 40 km04−Jul−2100:1001:4203:1404:4706:20ABDECevent A event Bevent C event DFigure 2.9: Examples of migration patterns in southern Vancouver Island. The shade ofeach epicenter dot is coded according to time. Dashed lines are, from left to right,the 20, 30 and 40 km contours of Audet’s plate interface model. For event C andD, templates used by and co-workers (Royer and Bostock, 2013; Royer et al., 2015)to identify these RTRs are shown (diamonds). Closed dashed contour lines with aletter label indicate asperities defined in figure Hypocenters and Depth DistributionA total of 818 hypocenters were obtained using P detections over the 3 ETS episodes from south-ern Vancouver Island. However, there is a high percentage of false detections due to significantlylower SNR on the vertical channel, occasional instrumental transients, and ambiguity causedby multiple sources separated by delays comparable to S-P delay times during intense tremoractivity. In order to extract meaningful depth information, hypocenter detections were visu-ally inspected and obvious false detections (e.g. from instrumental transients) were removed.318 hypocenters remain and are displayed in figure 2.10, projected onto a vertical dip plane.Hypocenters agree well with their associated epicenters. The horizontal location difference issmaller than 5 km for 90% of the detections and less than 3 km for 70%, which falls within the262.3. Resultslocation accuracy estimated by HYPOINVERSE. For comparison, template hypocenters fromRoyer and Bostock (2013) are also shown in figure 2.10. Only 12 detections are common tothe template catalogue and the depth difference with the corresponding templates is less than3 km.We see a strong clustering of detections along a plane consistent with the template detectionsand that parallels but sits between the plate interface models of Audet et al. (2010) and McCroryet al. (2012). At the downdip end, the hypocenters shallow slightly and approach the plateinterface of Audet et al. (2010).272.3.ResultsFigure 2.10: LFE hypocenters for southern Vancouver Island. a) Circles are the hypocenters obtained in this study using P detections,and their shade indicate the vertical error. The two perpendicular dashed lines and solid arrows indicate the referencecoordinates used to compute along-strike and along-dip distances, with the direction of arrows corresponding to positivevalues. Squares are templates (Bostock et al., 2012; Royer and Bostock, 2013). b) Depth with respect to along-dip distance.Circles are the hypocenters detected with the cross-station approach and squares are the template hypocenters. Dashed lineis the plate interface model of McCrory et al. (2012) and solid line is the model of Audet et al. (2010). Background shade isproportional to the density of hypocenters, with darker shades indicating high concentration.282.4. Discussion2.4 DiscussionA primary objective of this research was to automate the tedious manual selection of initialtemplates for network matched filtering in order to build a comprehensive template cataloguefor southern Vancouver Island and potentially other regions of Cascadia. Our method can targetspecific locations at the outset and is computationally efficient (∼1-3 hours per day of data)as compared to autocorrelation approaches (Brown et al., 2008) that require expensive patternrecognition with location as a final step. Using our cross-station detection time windows asinitial templates in the iterative network matched detection/stacking method of Bostock andRoyer (Bostock et al., 2012; Royer and Bostock, 2013), we have compiled additional high SNRtemplates for southern Vancouver Island. An example of initial and final templates are shownin figure 2.11 and 2.12, respectively. The final template reveals clearer P and S arrivals at manystations and represents detections from a location at (-124◦4’, 48◦29’, 30.5 km) at the seawardmargin of the large asperity A in figure 2.7 below southern Vancouver Island. The resultingtemplate has sufficient SNR to allow for determination of waveform polarity and so is suitablefor moment tensor inversion (Bostock et al., 2012; Royer and Bostock, 2013) and structuralstudies (Nowack and Bostock, 2013).When compared to catalogues that do not discriminate LFEs from tremor (Kao et al., 2006;Wech, 2010), our LFE catalogue displays similar spatio-temporal behaviour and in particularthe characteristic migration of the slow slip main front. Since LFE waveforms are of shorterduration, we can achieve a higher spatial resolution than tremor hypocenters, provided thatfalse detections can be effectively culled, for example by setting a threshold on permissiblelocation errors. The enhanced resolution is useful because it enables a more detailed study ofrupture area, small scale migration patterns (Rubin and Armbruster, 2013) and spatio-temporalclustering of seismic radiation at the plate boundary transition zone. However, the requirementof waveform similarity also limits the spatial range outside the network confines over whichdetections can be made relative to tremor detection via envelope similarity for longer durationbursts as in Obara (2002) and Wech and Creager (2008). We also note that whereas thesimplicity of LFE Green’s functions for southern Vancouver Island promotes detections basedon cross-station waveform similarity, this condition need not exist everywhere. For example,there is less consistency in waveforms across stations in northern Washington than in southernVancouver Island. However, even in regions where Green’s functions are more complex, ourapproach can still yield a useful suite of automatic detections based on the similarity of strong,abbreviated bursts of energy on LFE envelopes, although false detection rates may be higherand locations estimates poorer as a result of waveform distortions or cycle skipping. Theseautomatic detections can be subsequently evaluated visually and should be useful for templateconstruction.Examination of spatio-temporal density patterns in LFE detections reveals a highly het-erogeneous distribution during ETS, as previously observed for tremor in Washington stateusing seismic array analysis (Ghosh et al., 2009). The largest asperity (or near contiguous292.4. DiscussionFigure 2.11: Example of detection from our catalogue that is used as an initial templateto produce a high SNR template using the network cross-correlation and stackingprocedure of Bostock et al. (2012). The stations with an asterisk were used by thecross-station algorithm for detection and location.assembly of multiple asperities) A in southern Vancouver Island is located to the southwest,updip of the 30 km depth contour of Audet et al. (2010) plate boundary model . East of thisregion is the area of enhanced activity B for which Rubin and Armbruster (2013) detailed thebehaviour of small-scale slip fronts . A third asperity C is located below and immediately tothe north of station TWKB, less than 10 km updip from one of the 2 downdip asperities Dand E active between ETS episodes. The concentration of inter-ETS activity at the downdipend of the tremor zone in Cascadia has been previously reported by Wech and Creager (2011)and Armbruster et al. (2014), and by Obara et al. (2010) in Japan . In particular, the twointer-ETS asperities below Salt Spring Island and west of San Juan Island are also documentedby Armbruster et al. (2014). The fact that these asperities appear to be aligned with the ETSasperities in the along-dip, plate convergence direction prompts us to speculate that they arestructurally related to the shallower asperities through variable topography on the downgoing302.5. Conclusionsplate, as previously proposed for the Nankai subduction zone in Japan (Ide, 2010).We have also attempted to automatically detect P waveforms to gain insight into the LFEdepth distribution. Because the false detection rate is high, we have visually filtered out falsepositives and confirmed 318 valid detections. As shown in figure 2.10, our hypocenters dis-play a broader distribution in depth than templates hypocenters, but are generally consistentacross Vancouver Island (Bostock et al., 2012; Royer and Bostock, 2013). We explain the depthvariability in figure 2.10 by a) projection effects from 3-D structure, in particular due to thestrong curvature in the subducting slab and b) errors in our LFE P arrival times caused bythe generally high noise levels on the vertical component. Nevertheless, a pronounced con-centration of hypocenters corresponds with LFE template hypocenters and falls between theplate boundary models of Audet et al. (2010) and McCrory et al. (2012). The latter model wasassembled from a combination of constraints from intraplate seismicity and structural studieswhereas the former interprets the plate boundary to be the top of a well defined low velocityzone (LVZ) extending across the Cascadia margin (Hansen et al., 2012). Depending on whichplate boundary model is preferred, the depth distribution can be interpreted in different ways:if one considers the McCrory model to reflect the true plate position, then LFE sources aredistributed within the overriding plate, as argued by Kao et al. (2005). In contrast, if theAudet et al. model is adopted as correct, LFEs can be identified as a manifestation of currentplate motions at the plate boundary. In this context, we appeal to the interpretations Hansenet al. (2012) and Bostock (2012) wherein the LVZ is inferred to be weak upper oceanic crustalmaterial that acts as a distributed shear zone within which LFEs occur.RTRs on southern Vancouver Island were previously identified from templates by Royeret al. (2014). However, since spatial resolution was limited to the template distribution, thecorresponding time-distance plots possess a discretized, and possibly aliased appearance. Ourcross-station approach provides a more spatially continuous illumination of RTRs, and we areable to isolate the same events, but also new ones, as shown in figures 2.8 and 2.9. As evidentin Figure 2.9, the RTR migration direction can be quite complex. For example, RTR C seemsto migrate first updip and then downdip, as it proceeds along strike. Typical dimensions are oforder 50 km in the propagation direction and somewhat narrower in width.2.5 ConclusionsWe have developed an LFE detection strategy exploiting cross-station similarity and traveltimeconsistency of S-waveforms that differs from other algorithms (Armbruster et al., 2014; Rubinand Armbruster, 2013) in its use of a full complement of stations to survey extended regions. Todemonstrate its utility, we applied the method to datasets from southern Vancouver Island andWashington state. The LFE catalogues compare favourably with tremor catalogues for bothregions providing improved spatio-temporal resolution with somewhat reduced spatial cover-age. Spatio-temporal coverage is also consistent with information derived from a more discretespatial representation afforded by LFE templates (Bostock et al., 2012; Royer and Bostock,312.5. Conclusions2013). The more continuous spatial sampling of the cross-station approach allows us to bettercharacterize RTRs which can display complex propagation patterns. The lower SNR levels dis-played by P-waveforms exacerbate identification and phase measurement, however hypocenterdeterminations using P arrival time estimates cluster about those for LFE template waveformsand lie between the plate interfaces specified in two recent models. Cross-station detection ofimpulsive LFE waveforms has proven useful for template construction in northern Cascadia; itsutility in other regions will depend upon the degree of waveform similarity exhibited by stationsat smaller epicentral distances that, in turn, is influenced by the strength and scale of crustalheterogeneity.322.5. ConclusionsFigure 2.12: Final template obtained from the initial template of figure 2.11. The stationswith an asterisk were used by the cross-station algorithm for initial detection andlocation. The high SNR final template reveals P-wave and S-wave LFE waveformsat many more stations.33Chapter 3Seismicity, Metamorphism, and FluidEvolution Across the Northern CascadiaForearc3.1 IntroductionThe Cascadia subduction zone is a warm end-member among subduction zones globally. Theoceanic crust is 300-600◦C warmer in the forearc than cold subduction zones such as northeast-ern Japan, and is characterized by feeble volcanism, localized adakitic arc lava composition,and shallow seismicity. Because of warmer slab temperatures, eclogitization of the oceaniccrust proceeds at shallow depths, between 40-90 km (Kirby et al., 1996; Peacock and Wang,1999). Shallow dewatering reduces the water available for partial melting of the asthenosphere,and hence diminishes volcanic vigor. Warmer slab temperatures also influence seismogenesisthrough the prominent expression of slow-slip phenomena and restriction of intraslab earth-quakes to shallow depths. Slow-slip phenomena include slow slip, tectonic tremor and lowfrequency earthquakes (LFEs) (Dragert et al., 2001; Rogers and Dragert, 2003; Peacock et al.,2011; Bostock et al., 2012) along the plate boundary at depths of ∼ 30-40 km in northernCascadia. Tremor occurrence parallels the coast from northern California to northern Vancou-ver Island, with highest density beneath southern Vancouver Island/northern Washington, andnorthern California (Wells et al., 2017). LFEs are primary constituents of tremor (Shelly et al.,2007) and can be isolated using array processing techniques. LFE templates reveal thrust focalmechanisms (Royer and Bostock, 2014) consistent with plate motions, and moment magnitudesof 1.0-2.6 Mw in southern Vancouver Island (Bostock et al., 2015).Tremor epicenters (Kao et al, 2009; Boyarko and Brudzinski, 2010) and LFE hypocenters(Royer and Bostock, 2013; Plourde et al., 2015) in Cascadia are anticorrelated with those ofregular seismicity including both intraslab and crustal earthquakes. The subsurface expressionof the subducting Juan de Fuca plate provides a structural context for the different classes of343.1. Introductionseismicity. It is dominated by a thin low-velocity zone (LVZ) of extreme Vp/Vs ratio of ∼2.35 (Poisson’s ratio ∼ 0.39) extending under southern Vancouver Island between depths of 20and 45 km (Audet et al., 2009, 2010; Hansen et al., 2012). LFEs and tremor are inferred tooccur near the top of the LVZ, which is interpreted here as overpressured upper oceanic crustfollowing Hansen et al. (2012) and Bostock (2013). Where intraslab earthquakes occur (GeorgiaStrait/Puget Sound and northern California), they reside ∼ 7 km deeper near the oceanic Mohoand into the uppermost mantle (Preston et al., 2003; Bostock et al., 2012; Plourde et al., 2015;Kao et al., 2008) as in southwest Japan (Shelly et al., 2006). Moreover, their density is reducedin the tremor epicentral region. The epicentral distribution of crustal earthquakes also tendsto skirt that of the tremor generating region (Kao et al., 2009; Bostock, 2013).In this study we focus on structure and seismicity in southern Vancouver Island and theadjacent Washington-British Columbia mainland (see Figure 3.1). Most of Vancouver Island iscomprised of intrusive and volcanic assemblages, and sedimentary sequences that make up theWrangellia terrane, remnants of an oceanic plateau accreted during the Late Jurassic or EarlyCretaceous (e.g., Monger and Brown, 2016). At the southwestern tip of Vancouver Island, twosmaller terranes were later accreted during the Eocene. The Pacific Rim terrane, composedof metamorphosed sedimentary and volcanic rocks, is sandwiched between the San Juan andSurvey Mountain faults to the north and the Leech River Fault to the south. The Crescent(aka Siletz) terrane, composed of sea-floor pillow basalts and gabbros, is delimited to the northby the Leech River Fault, extends at depth to the south beneath a thick column of sedimentsmarking the Strait of Juan de Fuca and outcrops again in the Olympic Mountains in northernWashington. Sediments of the Nanaimo group overlie the southeastern coast of the Island,north of Victoria, and the Strait of Georgia. The eastern reaches of our study region includethe thick sediments of the Fraser Valley, bordered to the north by the mountains of the CoastPlutonic complex. Among previous studies of seismic velocity structure in the region, we citethe work of Clowes et al., (1995), Brocher et al. (2001), Zelt et al. (1996), Ramachandran etal. (2005), Ramachandran and Hyndman (2012) and McCrory et al. (2014). The depth of theforearc Moho is mapped near 36 km depth beneath the eastern reaches of the Fraser valley byZelt et al. (1996) consistent with receiver function estimates by Nicholson et al. (2005).The purpose of this work is to build upon these previous studies and investigate in moredetail the relationship between the major tectono-structural elements, different styles of seis-micity and fluid pathways within the subduction complex. We employ double difference tomog-raphy (Zhang and Thurber, 2003) that affords improved resolution of velocity structure nearearthquake source areas and improved seismicity relocation over conventional tomographic ap-proaches where seismicity is highly clustered such as in southwestern B.C. (Balfour et al., 2012).We image P-wave and S-wave velocity from multiple data sets to deduce Poisson’s ratio for astudy area similar to Ramachandran and Hyndman (2012). In the following sections we discussthe data and methods employed, and follow with a presentation of earthquake relocations andkey cross-sections of our model as represented by Poisson’s ratio and P-velocity (Vp). In the353.1. IntroductionFigure 3.1: Large-scale geologic map of study region. Semi-transparent colored patchesindicate geological terranes from BC Geological Survey data base: Crescent (CR),Pacific Rim (PR), Wrangellia (WR) and Coast Plutonic Complex (CPC). Orange ter-rane that includes the San Juan Islands and labelled SJ is classified as “Unassigned”by BCGS, but from Monger and Brown (2016) the San Juan Islands comprise mostlymid-Mesozoic and younger oceanic rocks from the Constitution formation (affiliatedto the Pacific Rim terrane) and the Fidalgo ophiolite complex. Gold solid lines areknown faults from the BC Geological Survey database (Leech River Fault (LRF), SanJuan Fault (SJF), Tofino Fault (TF) and Cowichan Fold and Thrust Belt (CFTB)and pink solid lines are from the USGS Fault and Fold database (Darrington-DevilsMountain Fault (DDMF), Hurricane Ridge Fault (HRF) and South Whidbey IslandFault (SWIF). Inverted triangles with red edges are permanent stations or tempo-rary arrays except POLARIS stations which are identified with green edges. Squareswith yellow outlines are stations deployed during the March 1998 SHIPS experiment.Solid Dark grey lines are plate depth contours (20, 30, 40 km) from Audet et al.(2010) and the dashed dark grey lines, from McCrory et al. (2012).discussion, we combine our observations with constraints from recent studies to present a newmodel for fluid evolution and crustal composition in this region that delivers new insights intothe nature of forearc crustal seismicity.363.2. Data and Methods3.2 Data and MethodsOur dataset comprises P- and S-wave arrival times for three types of events: regular intraplateearthquakes (within the forearc crust and subducting Juan de Fuca plate), low frequency earth-quakes (LFE) located near or at the plate boundary at depths of 30-40 km and air gun shotsfrom the 1998 Seismic Hazards Investigation in Puget Sound (SHIPS) experiment (Fisher etal., 1999; Brocher et al., 2001), described in more detail below. The study region encompassessouthern Vancouver Island along with the Juan de Fuca and Georgia Straits (latitude between47.75° and 50° N and longitude between 125.25° and 122.25° W).We employ the P- and S-wave first arrival picks of 4725 regular earthquakes cataloguedby Geological Survey Canada (GSC) between 1992 and 2012, with reported magnitudes (ML)between 0.01 and 3.46. Horizontal and vertical location errors and RMS residuals reported inthe catalogue are extracted and assigned as uncertainties for these data. We also employed theautomated cross-station algorithm of Savard and Bostock (2015) and modified its parametersettings to detect 333 new regular earthquakes between 2002 and 2006 that were not reportedby GSC. P- and S-wave first arrivals were manually picked at all available nearby stations froma collection of permanent and temporary networks that include the POLARIS array (Nicholsonet al., 2005), the Sooke, Sequim and Lopez small aperture arrays from the 2004 Deep TremorProject (La Rocca et al., 2005), the Plate Boundary Observatory Borehole Seismic Network, theCanadian National Seismograph Network, and some stations from the Pacific Northwest SeismicNetwork. For uncatalogued events, we assign uncertainties computed by HYPOINVERSE(Klein, 2002) for an initial hypocenter location based on manual picks using the 1D velocitymodel used by GSC for routine earthquake location. The final contribution to the earthquakedata set includes P- and S-wave picks from 274 events between 1980-2000 obtained from TomBrocher from a previous tomography study in N Washington (Preston et al., 2003) which werelocated in the southeastern corner of our study area. Initial hypocenters and uncertainties wereagain determined using HYPOINVERSE. For the entire earthquake dataset, the horizontalepicentral error has a mean of 0.7 km (standard deviation of 0.8 km) and the depth error hasa mean of 1.3 km (standard deviation of 1.8 km).LFE templates are stacked waveforms comprising 100’s-1000’s of independent LFE eventsdetected using network cross-correlation and yielding high signal-to-noise ratio waveforms withdiscernable P- and S-wave first arrivals at many stations (Bostock et al., 2012). Our LFEcatalogue for southern Vancouver Island comprises a total of 276 templates, 160 of these areconsidered independent, i.e. with less than 10% shared cross-correlation detections, obtainedusing the methods described in Bostock et al. (2012) and Savard and Bostock (2015). Thiscatalogue includes 2814 P-wave arrivals (average of 10 picks per template) and 3205 S-wavearrivals (average of 12 picks per template). LFE template hypocenters are valuable for tomog-raphy because their epicenters coincide with a region of the forearc with particularly low levelsof regular seismicity.Finally, we used P-wave first arrival picks for selected shot lines of the active-source 1998373.2. Data and MethodsSHIPS experiment (Ramachandran et al., 2005) recorded at temporary land stations. Theseshot lines run along the Strait of Juan de Fuca south of Vancouver Island into the Straitof Georgia to the west, extending as far north as Texada Island, and are assigned a timinguncertainty of 40 ms. We uniformly decimated this dataset by retaining only every thirdconsecutive shot for a given line, for a total of 2352 shots and 59329 P-wave picks, and foundthat using 33% of the SHIPS travel time dataset was a reasonable compromise between theinversion computation time and recovered model resolution. The SHIPS dataset is used to helpstabilize the tomographic inversion by constraining upper crustal velocity structure where fewearthquakes occur.We employ the double-difference tomography code tomoDD of Zhang and Thurber (2003),an adaptation of the double-difference hypocenter relocation algorithm of Waldhauser andEllsworth (2000) to the local earthquake tomography problem. Double-difference tomogra-phy exploits the proximity of event pairs and their traveltime difference at a given station toresolve nearby velocity nodes without bias from velocity anomalies near the stations (Zhangand Thurber, 2003). We simultaneously invert for P- and S- wave velocity and exploit bothcatalog and waveform cross-correlation time differences.Waveform cross-correlation delays are computed for all earthquakes between 2002-2006 andall LFE templates, for which we can exploit POLARIS array data. POLARIS and LFE datawere not employed by Ramachandran et al. (2005) or Ramachandran and Hyndman (2012),and provide significantly improved resolution of the forearc region. To reduce free surfacedistortion of waveforms at larger epicentral distances and to improve signal-to-noise ratios forcross- correlation analysis, we decompose 3-component seismograms into upgoing P-SV-SHwaveforms using the transfer matrix method of Kennett (1991). We use P waveforms to cross-correlate P-wave arrivals and SH for S-wave arrivals. Waveforms are filtered between 3-12 Hz,an effective pass band for phase arrival picking in southern Vancouver Island. We select 3-swindows within the predicted P and S wave arrival times (-1 s to +2 s relative to the predictedarrival), with a cosine taper of 0.5 s. Normalized cross-correlation functions are quadraticallyinterpolated near the maximum to provide delays at subsample precision. We use a lowerthreshold of 0.8 on normalized cross-correlation coefficients. This approach resulted in 30938and 58428 cross-correlation delays for P- and S-waves, respectively.The 3D P-wave velocity model of Ramachandran et al. (2003) is linearly interpolated on toour inversion grid with a horizontal spacing of 12 km and vertical spacing of 3 km. The choiceof a coarser grid was motivated by the heterogeneous distribution of seismicity in this regionand the large size of the dataset. This interpolated coarse model is used as our initial model forthe inversion, along with a uniform Vp/Vs ratio of 1.73 (Poisson’s ratio of 0.25). Several testinversions using a 1D model and a coarse grid as input produced final solutions converging tothose recovered with the 3-D model, albeit with larger anomalies and less detailed structure.The inversion employed 28 iterations, alternating hypocenter relocation with joint inversionof hypocenters and velocity. Our weighting scheme, based on the hierarchical approach of Zhang383.3. Resultsand Thurber (2003), employed absolute arrivals in the first iterations, then placed progressivelygreater weight on differential measurements, using first catalog delay times and then cross-correlation delay times, to refine small-scale structure. The weighted root mean square (RMS)traveltime residual for the catalogue data was reduced from 1.4747 s to 0.1195 s (-91%) whilefor the cross-correlation data it decreased from 0.2862 s to 0.0442 s (-84%) compared to theinitial interpolated 3D model with HYPOINVERSE locations.To image Poisson’s ratio, we culled our dataset by retaining arrival time data only forstations at which both P- and S-wave first arrivals were recorded for a given event (therebyexcluding, e.g., all SHIPS arrival time data). As a consequence the number of absolute P-wave traveltimes used is reduced by 70%, but the S-wave dataset is nearly the same (-15%).This reduction in data was required to prevent bias in measured Vp/Vs (Poisson’s ratio) valuescaused by uneven ray distribution and density (Evans et al., 1994; Tryggvason and Linde, 2006).For this culled dataset, the weighted catalog RMS residual is reduced from 1.1411 s to 0.0924s (-92%).We assessed the spatial resolution of our models using synthetic reconstruction tests since thescale of our inverse problem renders the construction and examination of the resolution matrixor the use of statistical methods such as bootstrap resampling (Rawlinson and Spakman, 2016)computationally impractical. Relevant sections for the recovery of a checkerboard pattern withvelocity anomaly amplitudes of 10% and the recovery of a synthetic slab with an imposed Vp/Vsratio of 2.35 are described and presented in Appendix B.3.3 ResultsFigures 3.2 and 3.3 show three depth profiles of Poisson’s ratio and Vp within the best- resolvedarea of the model, on which we will focus our attention. These profiles roughly coincide withstations of the linear POLARIS array across Vancouver Island (see Figures 3.1 and 3.4) andintersect a high density of low frequency earthquakes that improve our resolution of the overlyingforearc crust relative to previous models. One striking feature consistent in all 3 profiles is thehigh Poisson’s ratio region at the inferred mantle wedge corner formed by the plate interfacemodel of Audet et al. (2010) and the continental Moho at 36 km depth extrapolated fromregions to the east (Zelt et al, 1996; Nicholson et al., 2005). Poisson’s ratio hovers near 0.27-0.28 in this region, which is consistent with the range of 0.26-0.28 reported by Ramachandranand Hyndman (2012). However, a depth-localized low Poisson’s ratio region draped atop theforearc mantle corner in their model inferred to be silica enriched, is not evident in our results.Instead, Poisson’s ratio is notably low (∼ 0.22) through much of the forearc crust above themantle wedge corner and is associated with high average Vp, near 6.65 km/s (Figure 3.3).This low Poisson’s ratio region correlates spatially with increased levels of seismicity withinthe forearc crust, concentrated between 15-25 km depth. Previous teleseismic receiver functionstudies using POLARIS array data have shown that the average Vp/Vs of the forearc crustdisplays a statistically significant decrease with depth above the plate or, equivalently, landward393.3. Resultsdistance into the forearc (Audet and Bürgmann, 2014; Audet et al., 2010), consistent with theprofiles in Figure 3.2 between 60-200 km distance.Figure 3.2: Model Poisson’s ratio versus depth along the sections AA’, BB’ and CC’ shownin Figure 3.4. These profiles are taken parallel to the linear POLARIS array. Blackdiamonds are LFE hypocenters, black and pink dots are relocated earthquakes within10 km of the cross-section, with pink dots identifying previously uncatalogued earth-quakes detected by the cross-station algorithm. The thick gray line is the plateinterface model of Audet et al. (2010), extrapolated to depths >40 km as dashed,and dotted gray line is the plate interface model of McCrory et al. (2012). The hor-izontal dash-dot line is the nominal Moho depth of 36 km. Inverted black trianglemarks the SW coast of Vancouver Island.The extremely high Vp/Vs (∼ 2.35 or Poisson’s ratio ∼ 0.39) that characterizes the LVZas reported in Audet et al., (2009) and Hansen et al. (2012) and its spatial correspondencewith LFEs in Cascadia, is evident as a band of more modest (but elevated with respect tosurroundings) Poisson ratio of 0.26-0.27 (Vp/Vs of 1.76-1.78) that closely aligns with the slabmodel of Audet et al. (2010) from 20 km depth to its merger with high Poisson’s ratios inthe mantle wedge corner at ∼ 40 km depth. Due to regularization constraints, simultaneoushypocenter relocation, a paucity of underlying intraplate events, and the limited capacity oftravel-time tomography to resolve low-velocity zones compared to receiver-function methods,403.3. ResultsFigure 3.3: Model Vp along the sections AA’, BB’, and CC’ shown in Figure 3.4. Seecaption to Figure 3.2 for details.our tomographic model is expected to significantly smooth out sharp velocity contrasts (anddrive Vp/Vs to more normal values), but the feature is nevertheless recognizable across all threeprofiles. Almost all LFE hypocenters occur within this band of high Poisson’s ratio, betweenthe plate interface model of McCrory et al. (2012) and Audet et al. (2010).Relocated seismicity is shown in Figure 3.4 along with locations of major crustal faults.Crustal seismicity (in blue) is highly clustered east of the 30 km slab depth contour of theAudet model that sits slightly seaward of the inferred forearc mantle corner at 36 km depth.Seismicity landward into northern Washington and the BC lower mainland appears as well-defined clusters. Double-difference relocations indicate that these clusters are “streaks” ofseismicity with a more pronounced elongation in depth, consistent with an earlier study byBalfour et al. (2012). Seaward of the LFE source region, forearc seismicity is more evenlydistributed, with only 2 shallow clusters identified near [125.5°W, 49°N] and [125°W, 48.75°N]close to the west coast of Vancouver Island. Above the 20 and 30 km slab depth contoursthere is little to no forearc seismicity despite the seaward extension of Leech River Fault. Eastof approximately the 30 km depth contour, the Leech River Fault, recently shown to have413.3. Resultsbeen active in the Quaternary (Morell et al., 2017), is characterized by increased levels ofmicroseismicity.Figure 3.4: Relocated epicenters after inversion of the full travel time dataset. Red dots areintraslab earthquakes located below the plate interface model of Audet et al. (2010),and blue dots are crustal events located above. Positions of profiles AA’, BB’ andCC’ in Figures 3.2 and 3.3 are indicated. Black solid lines indicate known forearcfaults and thick dashed gray lines are the plate contours at 20, 30 and 40 km fromleft to right. Green diamonds are LFE epicenters. Yellow line is SHIPS ship track(e.g., Preston et al., 2003).Intraslab seismicity (in red) also appears to be divided into two distinct zones: diffuseseismicity seaward of the LFE source region and sparse intraslab seismicity that reappearslandward of the LFE source region, where the LVZ is expected to disappear due to increasedseismic velocities accompanying the onset of eclogitization of the subducting crust (Nicholsonet al., 2005, Bostock, 2013). Of interest in Figure 3.4, is a deep cluster of low magnitude(ML 0.06-3) intraslab events located southeast of Texada Island near 65 km depth below theGeorgia Strait. A profile through the cluster (not shown) indicates that this feature intersectsthe McCrory plate boundary model but lies ∼ 10 km below a smooth landward extrapolationof the Audet plate boundary which would place it near the oceanic Moho of the Juan de Fucaplate.Based on linkages between waveforms with high cross-correlation coefficients, 29 distinctclusters of events, including LFEs, can be identified using the GrowClust hierarchical clusteringapproach of Trugman and Shearer (2017). These clusters are plotted in Figure 3.5 along withtheir tomoDD relocations, with all events in a given cluster shown in a single color. Note that423.4. Discussionwaveform cross-correlation was only performed for events between 2002 and 2006, due to thelimited duration POLARIS deployment. With the exception of the shallow crustal cluster #1 and the aforementioned mantle cluster # 4, earthquake clusters are concentrated above themantle wedge, reaching depths of ∼ 25 km. The average local Poisson’s ratio in the vicinity ofthese forearc continental crustal clusters is 0.23 ± 0.02, for the deep slab mantle cluster (# 4)it is 0.22 ± 0.02, and for the shallow western cluster (# 1) it is 0.22 ± 0.01.Several forearc clusters in Figure 3.5 lie in close proximity to the Leech River Fault and/orDarrington-Devils Mountain Fault zone; namely clusters # 5, 6, 7, 10, 14, 15, 17, 18 and 24.Figure 3.6 shows 8 cross-sections of the Vp model, produced using the entire traveltime dataset, along with relocated seismicity. Morell et al. (2017) used field observations and lidar toargue that the Leech River fault is an active fault zone up to 1 km wide, extending 30 to 60 kmin length onshore. In eastern sections of the Leech River fault, where its trace changes strikefrom east-west to southeast-northwest and profiles aa’ and bb’ cross, these authors identifiedmany steep, reverse fault strands dipping northward between 60° and 90°, suggesting tectonicactivity during the late Pleistocene accommodating both strike and dip slip motion. Crosssection cc’ falls between the mapped traces of the Leech River and Darrington-Devils Mountainfault zones. Continuity of seismicity across bb’, cc’ and dd’ provides strong evidence that thetwo systems are the same. Many new events detected by our cross-station algorithm (bluestars in Figure 3.6) are 20-30 km deep at locations relative to the surface trace that would beconsistent with a steeply dipping fault in the range suggested by Morell et al. (2017). Lateralgradients in Vp in the upper 15 km also suggest the presence of a steeply north dipping faultbounding lower velocity material to the north. In profiles aa’ to ee’, the high velocity anomalyof ∼ 6.7 - 7 km/s south of the fault line at depths of ∼ 10-20 km likely corresponds to theCrescent (Siletz) terrane (Ramachandran et al., 2005). Many new events are located within theCrescent terrane at depths of 20-25 km beneath the Strait Juan de Fuca as shown in profile dd’.3.4 DiscussionIn combination with relevant constraints from previous studies, we proceed to frame our to-mographic Vp/Vs (Poisson’s ratio) anomalies and earthquake relocations in the context of amodel for fluid evolution as one proceeds landward off the Vancouver Island coast into the northCascade forearc. The conceptual model is illustrated schematically in Figure Slab Fluid Evolution From Trench to TremorFollowing Hansen et al. (2012), we assume that extensive hydrothermal circulation at/nearthe Juan de Fuca ridge results in elevated fluid concentrations within the upper oceanic crust.This inference is consistent with seismic refraction profiling of the Juan de Fuca plate offthe coast Washington and central Oregon (e.g., Horning et al., 2016; Canales et al., 2017).Canales et al. (2017) examined P-velocities within the Juan de Fuca plate along a margin-parallel refraction profile just seaward of the trench off Washington and northern Oregon. Their433.4. DiscussionFigure 3.5: Earthquake clusters linked by high waveform cross-correlation coefficients.Subfigures a) and b) show map locations and subfigure c) a depth profile alongthe red line in a) and b). Relocations reveal a tendency for crustal forearc seismicityto occur along vertical lineations consistent with the subset of clusters relocated byBalfour et al. (2012). Green diamonds are LFE hypocenters. The solid gray lineis the slab model of Audet et al. (2010) which is extrapolated as dashed to depthsgreater than 40 km. The dotted gray line is the model of McCrory et al. (2012).modelling indicates the Juan de Fuca lower crust and upper mantle are dryer than those for anyother subducting plate with most of the incoming water stored in sediments and upper oceaniccrust. In particular, Vp estimates for the Juan de Fuca plate off northern Washington areconsistent with a nominally anhydrous lower crust and upper mantle. Nedimović et al. (2009)documented ridge parallel normal faulting extending to Moho depths some 60 km landward ofthe Juan de Fuca ridge, indicating that localized hydration of the Juan de Fuca lower crust andmantle is possible. These authors further noted the potential of propagator wakes, anomalouslyhighly fractured areas of the plate, as localized conduits for mantle hydration.As the Juan de Fuca plate subducts and encounters higher pressures and temperatures, freewater is expelled from sediments, and prograde metamorphic reactions gradually liberate fluidsfrom hydrous minerals within the upper oceanic crust. These fluids are trapped within the plate443.4. DiscussionFigure 3.6: Vp model profiles and earthquake relocations along the eastern section of theLeech River Fault and western end of the Darrington-Devils Mountain Fault Zone.The middle of each section centers approximately on the fault trace, or, in the caseof profile cc’ which passes through Victoria, on an interpolation between the twofault systems. The thin vertical dashed line in each profile is the reference to thesurface fault trace. The oblique dashed lines represent northward dips of 63◦ and45◦ starting from the surface fault trace as labelled in the map plot. Thick graylines and horizontal black dash-dot lines delimit the Audet and McCrory slab modelsand continental Moho level (36 km), respectively. Events within 5 km of each profileare plotted: black dots are relocated earthquakes from the GSC catalogue, greendiamonds are relocated LFEs, blue stars are earthquakes newly detected in thisstudy; and larger colored circles are clustered events from Figure 3.5.soon after subduction through the development of a low permeability seal at the plate interface(Audet et al. 2009). The seal develops through some combination of deformation-inducedgrain-size reduction, the formation of phyllonites, and/or mineral precipitation thereby reducingpermeability to < 10−20m (Peacock et al., 2011). Fluid pressures rise to near-lithostatic values,leading to the seismic expression of upper oceanic crust as a high Poisson’s ratio and low VsLVZ (Christensen, 1984; Hansen et al., 2012), exhibiting strong reflectivity at seismic reflection453.4. DiscussionMt BakerDeformation frontMohoSerpentinizedforearc mantleHigh Poisson’s ratio / pore pressure reflective LVZ, ETS and LFEsCrustAccretionary prism100 km100 kmVancouver Island Northern Washington SGPacific Ocean   Accreted forearc terranesEarthquake nest: locally hydrated oceanic mantle (propagator wake?)Quiet zone:reduced forearc seismicityLow Poisson’s ratio / high Vp anomaly: quartz veins in greenschist/amphibolite facies mafic crust; enhanced & clustered forearc seismicityEclogitization of oceanic crusta)Upper crustLower crustHMBDMG ECLb)Figure 3.7: Schematic model of fluid evolution A) across northern Cascadia forearc, andB) within dehydrating oceanic crust where HMB is hydrated metabasalt, DMG isdry metagabbro and ECL is eclogite. SG is Strait of Georgia. Figure modified fromPeacock et al. (2011) with permission.frequencies (Hyndman, 1988; Bostock, 2013). Strong reflectivity is observed to commenceseveral 10’s of km offshore southern Vancouver Island coast at depths of ∼ 20 km (Nedimovićet al., 2003) and continues across the island. It and the LVZ persist to depths of ∼ 45 kmbeneath western Strait of Georgia (e.g., Nicholson et al., 2005). Recent teleseismic receiver-function work offshore southern Washington by Janiszewski and Abers (2015) suggests theLVZ there may persist farther seaward into the locked zone to depths of ∼ 16 km or possiblyshallower.Slow slip and its seismic expression in tremor and LFEs originate near or at the top ofthe LVZ and are observed at depths of 27-40 km (Royer and Bostock, 2014) in this region.It has been shown that tremor and LFEs are sensitive to stress changes of a few kPa at theplate interface caused by tides (Houston, 2015; Royer et al., 2015) and passing surface waves(Rubinstein et al., 2007; 2008), providing further evidence for the maintenance of high porepressures and low effective stresses at the top of the LVZ at these depths.Previous studies (e.g., Kao et al., 2009; Bostock et al., 2013) have noted a tendency towardanti- correlation in epicentral distribution of tremor/LFEs with regular forearc seismicity. Thisanticorrelation is evident in the relative paucity of regular crustal earthquakes above the tremorzone as displayed in Figures 3.2-3.4 is a consequence of the absence of slab-derived fluids enteringthe overlying North American crust due to the low permeability seal atop the slab. Hence weposit that fluids play a major role in promoting crustal seismicity within the Cascadia forearc.We note that Wells et al. (2017) have highlighted a statistically significant, along-strike anti-correlation between tremor density and the presence of major crustal faults, suggesting that463.4. Discussionthese deep faults are able to disrupt the plate boundary seal and, locally, tap fluids from theLVZ, thereby reducing tremor density. The Leech River Fault that extends to at least 25 kmdepth below southeast Vancouver Island, was one of the structures considered in their analysis.The deep clusters of microseismicity below this structure in the vicinity of Victoria, displayedin Figure 3.6, may therefore manifest localized fluid delivery from the slab. Hasegawa andNakajima (2017) have also suggested that the spatial anti-correlation of LFEs and high levelsof crustal seismicity observed at the Kii and Ise Gaps in southwest Japan, can be explained bya permeable forearc that drains fluids at the megathrust, with these fluids migrating throughcrustal faults, reducing their shear strength and promoting metamorphic reactions.3.4.2 Fluid Expulsion Into Forearc CrustAt depths near 45 km, the top of the LVZ begins to lose its velocity contrast with the overridingplate (Nicholson et al., 2005). As suggested by Audet et al. (2009), this process is thought tooccur as a result of eclogitization and a concomitant breach of the plate boundary seal thathas trapped fluids over ∼ 100 km of along-dip distance. Eclogitization results in a total solidphase volume reduction of ∼ 10% and is accompanied by production of fluids (e.g., Ahrensand Schubert, 1975; Hacker, 1996). In warm subduction zones, the reaction also involves areduction in the total product volume (fluid plus solid phases) resulting in void space (Aberset al., 2013). As a consequence, it will lower pore pressure, thereby discouraging generationof ETS. This line of reasoning implies that the onset of eclogitization controls the landwardlimit of tremor generation through the combination of reduced pore pressure and breached seal.Alternative explanations such as a change in frictional properties or compositional controls onpermeability originating at the mantle wedge corner, as has been proposed for southwest Japanby Katayama et al. (2012) and Okazaki and Katayama (2015) and for Cascadia by Hyndmanet al. (2015), are inconsistent with the observations that significant tremor generation occurs atcrustal depths (27 km and deeper) in Cascadia (Royer and Bostock , 2014) whereas it is confinedto dominantly mantle depths (40-60 km depth) in eastern Alaska (Chuang et al., 2017).Evidence for the breach of plate boundary seal comes from three sources. First, the absenceof continental Moho near the mantle wedge corner implies that fluids have penetrated themantle wedge and caused extensive serpentinization, effectively erasing the velocity contrastwith overlying crustal material in the North American plate (Nicholson et al., 2005; Bostock etal., 2002; Brocher et al., 2003). Second, a pronounced landward increase in 3He/4He ratiosmeasured in mineral springs above the mantle wedge corner in northern Washington Statesuggests a punctuated influx of fluids from the Juan de Fuca plate (McCrory et al., 2016).This observation also requires that slab-derived fluids flux through the mantle wedge into thecontinental crust. Thus, although it has been demonstrated experimentally that antigorite candevelop strong permeability anisotropy under shear (Kawano et al., 2011), this effect does notappear to prevent fluid progress through the Cascadia mantle wedge. Finally, the appearanceof highly clustered, crustal seismicity with swarm-like characteristics to depths of 25 km above473.4. Discussionand landward of the downdip termination of the LVZ and tremor/LFEs (Figure 3.5) impliesthat these earthquakes are facilitated by fluid ingress at depth (e.g., Vidale and Shearer, 2006).Before proceeding, we remark that increased levels of forearc crustal seismicity in the southernVancouver Island – Puget Sound region (as well as northern California) are approximatelymirrored updip by the regions of highest tremor density (Wells et al., 2017). The inferencethat crustal seismicity is enabled by fluids suggests that tremor density too is dependent uponabsolute fluid volumes transported within the subducting plate.3.4.3 Forearc Crustal CompositionWe now consider the effect of fluid ingress from the Juan de Fuca plate on the compositionof the overriding continental crust in the vicinity of the wedge corner. As evident in Figure3.2, this region of continental crust is characterized by Poisson’s ratios (nominal value 0.225,see Appendix B, Figure B.6) that are considerably lower than adjacent regions and the globalaverage (Christensen, 1996). Moreover, this signature occurs in combination with high veloci-ties, that is, high crustal Vp (∼ 6.65 km/s; see Figure 3.3 and Figure B.6 in Appendix B) andhence proportionately higher Vs (∼ 3.9 km/s; see Figure B.1 in Appendix B), an associationthat favors a compositional origin over one dominated by a simple fluid presence since velocitiesshould decrease in the latter case. If the two degrees of freedom are represented by Vp andPoisson’s ratio, we may examine candidate compositional mixtures by plotting these nominalvalues with laboratory measurements by Christensen (1996) at 400 MPa (results at 200 and600 MPa yield similar results). Assuming linear, two-phase mixing using average end-membervalues, it becomes evident that there are limited scenarios as to mixtures that can explain ourmeasurements. All scenarios require quartz to elevate Vs since quartz has, by a wide margin, thelowest Poisson’s ratio of commonly occurring minerals at 0.1 (Christensen, 1996). It is notablethat the two best matching scenarios involve greenschist facies metabasalt (BGR) and/or am-phibolite (AMP) as the dominant endmember. The tomographic anomaly can thus be modelledas “average” greenschist facies metabasalt and/or amphibolite with the addition of substantialquartz (∼ 15% when averaged through traveltimes). Compositional deviations from these av-erages will modify the required quartz content. Surface geology above the anomaly is markedby a transition from Wrangellia (island arc assemblages of intermediate to mafic compositions)in the Gulf Islands in the northwest to oceanic rocks and ophiolitic complexes (dominantlymafic compositions) in the San Juan Islands to the southeast (Monger and Brown, 2016), andso is broadly consistent with a model interpretation involving greenschist facies metabasalt andamphibolite.3.4.4 Quartz Concentration MechanismsThe suggestion that low Poisson’s ratios in crustal forearcs could be the result of silica precipi-tation from slab-derived fluids is due to Ramachandran and Hyndman (2012) and Hyndman etal. (2015). They employed this mechanism to explain a low Poisson’s ratio anomaly just above483.4. DiscussionFigure 3.8: Plot of forearc crustal anomaly Poisson’s ratio and Vp (red ellipse) on cor-responding laboratory measurements for various crustal lithologies as measured byChristensen (1996). Labels AMP, BGR, DUN, SER, QTZ correspond to amphibo-lite, greenschist basalt, dunite, serpentinite and quartz. See Christensen (1996) forother label definitions. The red ellipse represents average values for the low Poisson’sratio anomaly in the forearc crust above the mantle wedge corner that is discussedin the text. The blue dashed line indicated the approximate expected values fortwo-phase mixing of quartz with either greenshist basalt or amphibolite.the mantle wedge obtained using the Vp model of Ramachandran et al. (2005) and a subse-quently derived Vs model. Hyndman et al. (2015) suggested that rising fluids, containing asmuch as 1-10% silica from the subducting oceanic crust landward of the wedge corner are chan-neled toward the wedge corner by an impermeable, serpentinized, forearc mantle. Since silicasolubility increases rapidly with temperature, these authors suggest that silica-saturated fluidsaccess the lower continental crust updip of the mantle wedge corner and precipitate quartz asthey rise. We note several aspects of their silica precipitation model that need to be reconciledwith observations made here. First, it does not acknowledge the presence of an impermeableseal at the top of the plate up- (and down-) dip of the mantle wedge corner to maintain thefluid overpressures that accompany LFEs and tremor (Peacock et al., 2011). Second, it relieson the mantle wedge corner as a hydrological control that, as noted above, is inconsistent with493.4. DiscussionLFE/tremor observations across different subduction zones (Chuang et al., 2017; Royer andBostock, 2014; Shelly et al., 2006). Third, it employs a slab geometry across southern Vancou-ver Island that locates the majority of intraslab seismicity within oceanic crust implying thatLFEs occur within the mantle wedge and that the wedge corner is located ∼ 40 km seawardof its inferred position in this study (see Ramanchandran and Hyndman, 2012, their Figure 4).Fourth, their low Poisson’s ratio anomaly lies deep within their crustal model, extends laterallyover ∼ 100 km and, although continuous, is due to high Vs at its seaward edge and low Vp at itslandward edge. This complex signature may result from a limited resolution that is improvedupon in the profiles of Figures 3.2,3.3 through our inclusion of i) LFE data corresponding toa region largely devoid of regular earthquakes and ii) data from diverse classes of seismicityrecorded on the densely sampled POLARIS array (Nicholson et al., 2005). Note that the profilesin Figure 3.2 suggest that low Vp/Vs persists vertically through much of the North Americancrust, and over a more limited horizontal extent, near 40 km (see supplementary information).This spatial distribution is not consistent with the slab silica precipitation model where a strongpositive gradient in silica content with crustal depth is predicted. Finally, as discussed by Hyn-dman et al. (2015), there exists a significant mass-balance discrepancy between the amountof quartz required to explain the volume over which low Vp/Vs is observed and that availablefrom the downgoing plate (see below). This discrepancy is further exacerbated by the presenceof a silica- impoverished mantle wedge that likely filtered significant quantities of silica fromslab fluids en route to the crust, to become serpentinized. The imbalance need not affect theargument by Audet and Bürgmann (2014) that quartz deposition controls the periodicity ofETS recurrence since only minor volumes (but expansive areas) of quartz are necessary to affectfrictional properties along the lower edge of the plate boundary.Let us consider for a moment, the slab silica deposition model in the context of the fluiddelivery constraints discussed in section 4.2. More specifically, we rely on i) expulsion of fluidsfrom the slab due to eclogitization at P/T conditions predicted to initiate below seaward por-tions of the mantle wedge in Cascadia (Peacock and Wang, 1999; Rondenay et al., 2010), andii) the focused, upward migration of these fluids into the overriding crust where they depositquartz and facilitate seismicity. The mass-balance problem discussed by Hyndman et al. (2015)arises as follows. If we assume mean values for representative endmember lithologies, traveltimeconstraints imply the addition of ∼ 15% quartz by volume (or a 6 km thick layer) within acrustal column of 40 km width corresponding to the landward extent of the low Poisson’s ratioanomaly. This number must be matched by the silica budget of the subducting crust over 10to 50 Ma of subduction at an assumed rate of 4 cm/a. The solubility of silicon at depths belowthe mantle wedge is 0.356mol kg−1 of H2O (Manning, 2004), the molar volume of quartz is2.268× 10−5m3mol−1 and the analysis of Peacock et al. (2011) suggests 3% H2O by volumewithin a ∼ 3.0 km thick upper oceanic crust. These numbers imply 4× 105m3 of SiO2 per me-tre of strike for the 10 Ma figure. When distributed over a 40 km along-dip interval, we arriveat an equivalent layer thickness of SiO2 of 10 to 50 m depending on duration of subduction.503.4. DiscussionThese estimates are between 2 and 3 orders of magnitude too small to explain the 6 km seismicprediction, and support the argument that silica addition to the forearc continental crust fromslab fluids is not the major contributing factor to the Poisson’s ratio anomaly in figure 3.2.The work of Breeding and Ague (2002) was cited by Hyndman et al. (2015) to support theirthesis of quartz addition to the forearc crust through precipitation from slab fluids. The formerstudy employed SiO2−Zr systematics from greenschist-facies metamorphosed sediments of theOtago Schist, containing 10-30 vol% quartz veins, to argue that silica was externally derived. Acontrasting study of quartz vein networks in a similar accretionary environment at greenschistconditions (Kodiak Island, Fisher et al., 1995) found no evidence for long distance-transport ofsilica, although the fluids responsible for quartz mobilization may have travelled long distances.Rather, geochemical measurements indicated local diffusion of silica from the rock matrix toquartz veins with the amount of wall rock silica depletion comparable to the amount of quartzprecipitated. This latter study raises the possibility of significant in situ quartz generationthrough metasomatic reactions facilitated by ingress of slab fluids. An illustrative reactionrelevant to metabasalt at greenschist facies conditions that may characterize the Gulf / SanJuan Islands region at depth, is the sericitic alteration of albite modeled by:NaAlSi3O8albite+ 2H+(aq) +K+(aq) = KAl3Si3O10(OH)2sericite+ 6SiO2quartz+ 3Na+(aq) (3.1)that produces significant silica along with a ∼ 7% volume reduction. Given the massbalance issues associated with large-scale addition of slab-derived silica to the forearc, we suggestthat the tomographic Vp and Vp/Vs anomaly results from in-situ quartz generation throughmetasomatism aided by fluids fluxed from the slab downdip of tremor (Figure 3.7).3.4.5 Crustal Forearcs and Orogenic Gold DepositsOur observations and interpretations are reminiscent of the association between quartz, seismic-ity and greenschist facies metamorphism in orogenic gold quartz-vein systems that Sibson andcoworkers (e.g., Sibson et al., 1988; Sibson and Scott, 1998) have explained using a fault-valvemodel. Quartz vein networks develop syntectonically in horizontal compressional or transpres-sive regimes. They involve vein emplacement on high-angle reverse faults and subhorizontalextensional fractures. High-angle reverse faults are unfavorably disposed to failure in this stressenvironment. Their activation and the generation of extensional fractures are possible only inthe presence of supralithostatic, pore-fluid pressures exceeding the minimum (vertical) com-pressive stress. Accordingly, the fault-valve model involves a cycle of pore- pressure build-up,fault-strength exceedance accompanied by seismic rupture, fluid ascent enabled by increasedpermeability, and pore pressure reduction accompanied by mineral precipitation that in turndecreases permeability. Estimates of earthquake magnitudes for the Val d’Or gold field of theAbitibi Greenstone Belt from quartz vein observations by Robert et al.(1995) correspond to ML<3-4. Our observations can be considered in this context since the maximum principal stress513.4. Discussionin northern Cascadia is horizontal and oriented along strike (e.g., Wang, 2000; Balfour et al.,2011). The elongated vertical extent of most forearc crustal earthquake clusters in Figure 3.5is consistent with slip on high-angle reverse faults, requiring locally high pore-fluid pressures.Furthermore, the range of magnitudes for events in the crustal forearc region used in this studyis comparable to those quoted above at ML 0.01-3.46.Quartz vein networks within orogenic gold deposits range from Archean to Tertiary inage. The processes involved in concentrating gold to form these deposits must act over re-gional scales, that is, dimensions of ￿100’s of km, and involve greenschist to amphibolite faciesmetamorphism, close to the brittle-ductile transition (Phillips and Powell, 2010; Goldfarb andGroves, 2015). Goldfarb and Groves (2015) have evaluated various models proposed over thepast few decades for the genesis of orogenic gold deposits. Models invoking gold-bearing fluidsreleased in the middle crust through greenschist-to-amphibolite facies metamorphism (Phillipsand Powell, 2010) and focused upward into regional fault systems, are generally consistent withavailable data, whereas magmatic-hydrothermal origins are dismissed on the basis of geologicand geochronologic constraints. Given a significant representation of orogenic gold depositsback into the Precambrian and a likely greater prevalence of warm subduction conditions inearlier times, it is interesting to speculate as to whether these deposits may have arisen withincrustal forearcs through metamorphism and metasomatism enhanced and enabled by shallow,warm slab dehydration.3.4.6 Leech River Fault and Texada Island SeismicitySibson (2009) has applied the fault-valve model to larger (M￿6) earthquakes in the backarcof northeast Japan. These events occur with unusual frequency on high-angle thrust faultsthat are reactivated normal faults formed in earlier extensional episodes. Local geophysicalstudies often indicate an association with high-pore pressures. It is interesting to note thatthese seismogenic structures are observed within the backarc of a cold subduction zone settingwhere the dominant release of water occurs at considerably greater depths (Peacock and Wang,1999; Kirby et al., 1996), contrasting distinctly with the focused, forearc fluid flux envisagedfor the warm subduction conditions of northern Cascadia. We note that comparably largecrustal earthquakes may also be of regular occurrence within the Cascadia forearc. The LeechRiver fault forms the boundary between the Pacific Rim terrane (Leech River Schist) andthe Crescent (or Siletz terrane). Offshore Vancouver Island, its trace adopts a strike-parallelgeometry that bends abruptly as it enters the study area to assume an east-west and thensoutheast-northwest trajectory on eastern Vancouver Island (Figure 3.1; see also Hyndmanet al., 1990). Consequently, it is likely that the Leech River fault has undergone significantdeformation and reorientation along its onshore segments (e.g., Fairchild and Cowan, 1982).The seismicity profiles in Figure 3.6 imply that the Leech River fault is steeply dipping at ≥60 ◦possibly with multiple strands, thus requiring high pore-fluid pressures to allow microseismicityin the presence of an arc-parallel, compressional stress field (Sibson, 2009; Wang, 2000). Recent523.5. Conclusionswork by Morell et al. (2017) using lidar and field observations has revealed that this sectionof the Leech River fault has experienced at least two surface-rupturing M≥6 earthquakes inthe past 15 ka. This region lies along the eastern margin (km 120-140) of the low Poisson’sratio anomaly in profile AA’ of Figure 3.2. We suggest the deeper (∼25 km) microseismicity inprofiles aa’ to dd’ within Figure 3.6 occurs where medium-grade metamorphic reactions definethe base of the brittle-ductile transition and where future large crustal earthquakes will likelynucleate again (Ito, 1999; Sibson, 2009).To complete our discussion of the North Cascadia water cycle we note that the Juan de Fucaintraplate earthquake cluster # 4 between (58-69 km) below Texada Island is well represented(≥ 88 events between 1992-2012) and isolated from other intraslab seismicity that tends to beconcentrated to the south and east below Puget Sound. It lies along the subsurface extrapolationof a prominent propagator wake (Nedimović et al., 2009). These structures are likely to produceincreased permeability and localized fluid access to mantle depths within a Juan de Fuca platethat is otherwise anyhdrous at these levels (Canales et al., 2017). Accordingly, we suggest thisearthquake cluster originates through dehydration embrittlement as locally hydrated mantletransforms from antigorite back to peridotite (e.g., Peacock, 2001).3.5 ConclusionsWe have employed double difference tomography to elucidate the relations between microseis-micity, fluid pathways and metamorphism in the northern Cascadia forearc. Much of ourinterpretation draws from results of earlier studies; for example, the identification of LVZ withoverpressured upper oceanic crust (Audet et al. 2009; Peacock et al., 2011; Hansen et al., 2012;Bostock, 2013), fluid infiltration of the mantle wedge (Bostock et al., 2002; Brocher et al., 2003),the identification of low Poisson’s ratio crust with high concentrations of quartz (Ramachan-dran and Hyndman, 2012; Hyndman et al., 2015) and the temporal modulation of slow slipphenomena through silica precipitation (Audet and Bürgmann, 2014). Here we have focused onunresolved issues concerning i) the long-recognized spatial anticorrelation between epicenters oftremor/low frequency earthquakes and crustal earthquakes in northern Cascadia (e.g., Kao etal., 2009), and ii) the mass balance discrepancy that arises in explaining low forearc Poisson’sratios through silica precipitation of slab-derived fluids (Hyndman et al., 2015; Ramachandranand Hyndman, 2012; Audet and Bürgmann, 2014). Specifically, we suggest eclogitization ofoverpressured upper oceanic crust leads to a breach in the plate boundary seal and focusedexpulsion of slab fluids into the overlying mantle wedge and North American crust. These slab-derived fluids may deposit minor silica but their primary effect is to promote metasomatismat greenschist and lower amphibolite facies within deforming continental crust. Their passagethrough a ductile and low permeability wedge and lower continental crust promotes high-porepressures that enable microseismic reactivation of steeply dipping thrust faults (e.g., Sibson etal., 1998; Sibson, 2009), ensuing fluid ascent and concentration of quartz through pore pressurecycling. Accordingly, the forearc crust in warm subduction zone settings is posited to represent533.5. Conclusionsa global extreme as regards sustained exposure to incoming fluids. The anticorrelation of micro-seismicity and tremor observed in other locations (southwest Japan, Nakajima and Hasegawa(2016); northern California, Boyarko and Brudzinski (2010), Plourde et al. (2015); and easternAlaska, Wech (2016), Chuang et al. (2017)) and genetic associations between seismicity, fluidsand quartz vein networks suggest that the model detailed here may have broad relevance towarm convergent margin processes and the development of orogenic gold deposits.54Chapter 4Double-difference Tomography of NorthernVancouver Island4.1 IntroductionThe area of northern Vancouver Island is of particular interest as it is situated at the northernterminus of Cascadia subduction and its tectonic setting differs considerably from southernVancouver Island, as shown in Figure 4.1. The northern edge of the Juan de Fuca platewas broken into a micro-plate, called the Explorer plate, through the creation of the Nootkatransform fault zone around 4 Ma (e.g., Rorh and Furlong, 1995). To the north, the QueenCharlotte Triple Junction (QCTJ) is formed by the Queen Charlotte Transform, the Dellwood-Tuzo-Wilson-Knolls spreading system and the subducting Explorer micro-plate. In contrastto southern Vancouver Island, there is little seismicity in the onshore forearc, but extensiveoffshore seismicity, principally along the Nootka fault zone and the western boundary of theplate, but also throughout the Explorer plate, indicating a high degree of intraplate stress. Itis suggested that major internal shear deformation within the western and northern part of theExplorer plate is effectively fragmenting the plate as it accommodates an ephemeral adjustmentbetween the Pacific and North America plates (Rorh and Furlong, 1995; Kreemer et al., 1998;Dziak, 2006). Since its formation at approximately 4 Ma, the convergence rate of the Explorerplate has significantly decreased to a nearly steady value of about 19 mm/yr for the past 1 Maaccording to sea-floor magnetic anomalies, which contrasts with the convergence of 43 mm/yroff central and southern Vancouver Island (Riddihough, 1977; Wilson, 1993). The difference inconvergence rate is accommodated by deformation along the Nootka Fault zone, a 20 km wideshear zone with left-lateral transform motion (Hyndman et al., 1979). The age of the Explorerplate is estimated to be about 5 Ma at the deformation front and underneath Vancouver Island(Wilson, 1993, 2002; Davis and Riddihough, 1982; Davis and Hyndman, 1989; Gao et al., 2017)making it an end-member “warm” subduction zone.554.1. IntroductionFigure 4.1: Tectonic context of the study region. Dark green dots are earthquakes ofminimum magnitude 2 from 1960 to 2018 retrieved from the IRIS DMC FDSN WebService (U.S. Geological Survey, Earthquake Hazards Program (2017); ISC On-lineBulletin (2015); Dziewonski et al. (1981); Ekström et al. (2012), and light green dots are earthquakes recorded from1991 to 2002 by the Sound Surveillance System (SOSUS) (Fox et al., 1994). Principalstructures are labelled: JdF—Juan de Fuca plate, GA—Gorda micro-plate, EX—Explorer micro-plate, PA— Pacific plate, NA—North America plate, BTF—BlancoTransform Fault zone, STF—Sovanco Transform Fault zone, DTWK—Dellwood-Tuzo Wilson Knolls spreading system, QCF—Queen Charlotte Fault, NFZ—NootkaFault Zone, MTJ—Mendocino Triple Junction and QCTJ—Queen Charlotte TripleJunction. Pink triangles are Cascade volcanoes and purple triangles are exposedAlert Bay volcanics of Tertiary age (Armstrong et al., 1985; Lewis et al., 1997; Cuiet al., 2017). Thick gray lines are the 20, 30 and 40 km depth contours from theJuan de Fuca slab from left to right respectively Audet et al. (2010). The yellowcontours indicate the extent of onshore tectonic tremor (Wech, 2010). Ridges arelabelled and shown as the thickest lines. The blue rectangle is the study area ofChapter 3, whereas the green rectangle highlights the study region for this Chapter.Plate convergence rates from Wilson (1993)564.1. IntroductionFigure 4.2: Regional seismo-tectonics of northern Vancouver Island. Same as figure 4.1but with moment tensor solutions from Braunmiller and Nábělek (2002); Splinder etal. (1997); Cassidy et al. (1988); Rogers and Hasegawa (1978).The basement of Vancouver Island is composed mainly of the Wrangellia terrane, an oceanicplateau formed during the Paleozoic-early Mesozoic and accreted onto North America during themid-Cretaceous, ca. 100 Ma (e.g. Monger and Journeay (1994) and references therein). In theEarly Eocene, two smaller terranes collided and partly underthrust Wrangellia: the Mesozoicsedimentary Pacific Rim terrane and the Eocene mafic volcanic Crescent terrane (references inHyndman et al. (1990)). The West Coast fault, Tofino fault and inferred Crescent Thrust faultdefine the boundaries between the three terranes and the accreted sedimentary wedge (Hyndmanet al. (1990), also see Figure 4.3). The Pacific Rim terrane is mapped along the entire margin ofVancouver Island, but the northern end of the underthrusting Crescent terrane is situated justoff Hesquiat peninsula (see figure 4.3 and 4.4). Recently, Johns et al. (2012) found evidence ofa fragment of the Crescent terrane that detached and moved north and landward, and is nowlocated just under the coast, between Brooks Peninsula and Nootka Island. On the continentalshelf, the Pacific Rim and Crescent terranes and the accreted sedimentary wedge are overlainby an up to 4000 m thick blanket of Eocene to Holocene sediments that form the Tofino basin,extending from southern Vancouver to Brooks Peninsula (Hyndman et al. (1990); Johns et al.574.1. Introduction(2012) and references therein). These sediments outcrop on the western edge of Nootka Islandand Hesquiat Peninsula as the Carmanah Group (see Figure 4.4). North of Brooks Peninsula,the Winona Basin is characterized by young (Pliocene-Pleistocene) sediments overlying a young,warm and buoyant fragment of oceanic crust with age varying from 1 Ma at its southern end to4 Ma in the north where it meets the Dellwood Knolls, forming a structure called the WinonaBlock. The Winona Block likely formed and potentially became independent from Explorerwhen the Pacific-America-Explorer triple junction jumped from its previous position prior to1 Ma just off Brooks Peninsula to its current position just south of Haida Gwaii (Davis andRiddihough, 1982; Riddihough et al., 1980; Murray and Tiffin, 1974)) (see Figures 4.1 and 4.3).These authors also argued that northwesterly striking folds and ridges in the sediments areindicative of continued underthrusting of the Winona Block beneath the continental margin ata rate likely equal to that of the Explorer plate at its southeast end transitioning to null at itsnorthern end.A more detailed geological map of Vancouver Island is shown in figure 4.4. The Wrangelliaterrane is mainly characterized by the flood basalts of the Karmutsen Formation (KV), un-derlying and overlying limestone from the Quatsino Formation. It overlies the older Devonianvolcanic rocks of the Sicker and Buttle Lake groups (Sk). Volcanic activity during the EarlyJurassic formed the calc-alkaline volcanic rocks of the Bonanza formation (Bo). Finally, thelast component of the terrane was formed during the Jurassic, when volcanic intrusions createdthe granodioritic batholiths of the Island Plutonic Suite (IP) and the West Coast CrystallineComplex (WC) (Jones et al., 1977; Greene et al., 2010; Earle, 2016). Along the east coast nearPort Hardy and Campbell River, the Nanaimo Group (Na) is composed of a 5000 m thick layerof sediments deposited during the Late Cretaceous.584.1. IntroductionFigure 4.3: Regional faults and major terranes of the study region. Terranes are indi-cated by colored shades: WR, Wrangellia terrane (orange); PR—Pacific Rim terrane(blue); CR—Crescent terrane (purple) and CPC—Coast plutonic Complex (green,at the northeast edge of the map). Faults are colored according to their classificationwithin the BC Geological Survey Fault database (Cui et al., 2017): red lines arethrust faults, purple are steeply-dipping faults (type unspecified), green are normalfaults, yellow are extension faults and black are faults of unspecified nature in thedatabase. Bold fault lines are named and their abbreviation is indicated by blackarrows: WCF—West Coast Fault, BRTF—Beaufort Range Thrust Fault, BCF—Beaver Cove Fault, BLF—Bonanza Lake Fault, BIF—Browning Inlet Fault, GRF—Goospeed River Fault, HF—Holberg Fault, HLF—Hustan Lake Fault, KCF—KwaisCreek Fault, LMLF—Le Mare Lake Fault, MCF—Mahatta Creek Fault, KRF—Klaskish River Fault, NRF—Nahiwitti River Fault, NIF—Neroutsos Inlet Fault,NF—Nimpkish Fault, SCF—Spencer Cove Fault, WLF—William Lake Fault. Thedotted bold black line defines the Brooks Peninsula Fault Zone (BPFZ) comprised ofa corridor of steeply-dipping, northeasterly-trending, sub-parallel faults such as KRF(Nixon et al., 2011). Red triangles are exposed Alert Bay volcanics.594.1. IntroductionFigure 4.4: Major geological groups of northern Vancouver Island. As in figure 4.3, CR—Crescent terrane, PR—Pacific Rim terrane. The major units of the Wrangelliaterrane are shown: Bo—Bonanza volcanic rocks, IP—Island Plutonic Suite, KV—Karmutsen volcanics and Quatsino limestone, WC—West Coast crystalline complex,Na—Nanaimo sediments, Ca—Carmanah sedimentary rocks, Sk—Sicker and ButtleLake volcanic rocks, QC—Queen Charlotte group sedimentary rocks. Layers ex-tracted from BC Geological Survey database (Cui et al., 2017). Pink triangles areexposed Alert Bay volcanics.The complex tectonic history of the region is reflected by a contrast in physical properties ofthe continental forearc as defined by the boundary formed by the Brooks Peninsula Fault Zone(BPFZ), (see Figure 4.3) a corridor of northeasterly-trending, steeply-dipping faults landwardof the Brooks Peninsula that cuts through several northwesterly trending faults (Muller et al.,1973). Lewis et al. (1997) note an abrupt change in heat flow on northern Vancouver Island604.1. Introductioncoinciding with this boundary, ranging from 46mWm−2 just south of BPFZ similar to centraland southern Vancouver Island, to a value of 67mWm−2 north of the BPFZ, similar to theQueen Charlotte Basin measurements. The increase in heat flow coincides with a lower meanelevation and more subdued topography north of the BPFZ, evident in Figure 4.5. Slightlyhigher Bouguer gravity anomalies and lower seismicity are also noted for the northern tip ofthe Island. The low heat flow for southern Vancouver Island is due to the heat sink createdby the subducting Juan de Fuca plate at 30-40 km depth beneath the Island (Lewis et al.,1988; Hyndman and Wang, 1993). By considering the time constant for this thermal effectto reach the surface, Lewis et al. (1997) argued that the change in heat flow represents thenorthern edge of the subducting Explorer plate, previously inferred by Riddihough (1977) fromseafloor magnetic anomalies. The heat flow observations also suggest that the triple junctionformerly off Brooks Peninsula was stationary for 40 Ma. This proposed location for the plateedge also approximately coincides with the occurrence of the Alert Bay Volcanic Belt, a seriesof isolated volcanic flows, pyroclastics, plugs and dikes of Tertiary age. It has been suggestedthat faulting activity in the BPFZ of similar age is related to this volcanic activity (Muller etal., 1973). Armstrong et al. (1985) proposed that this late Tertiary volcanic episode, startingca. 8 Ma and peaking at 4 Ma, is a plate edge phenomenon, caused by a disruption in steady-state subduction patterns at the plate edge. Lewis et al. (1997) argued against this hypothesisbased on a geochemical affinity of the Alert Bay volcanics with volcanic rocks from the MassetFormation in the Queen Charlotte Basin. They argued instead that the northern tip of theIsland shared a Tertiary episode of extension with the Queen Charlotte Basin, implying athinning crust north of the BPFZ.In 1980, the two 2D controlled-source refraction seismic profiles with both onshore andoffshore components, were shot on Vancouver Island as part of the Vancouver Island SeismicProject (VISP). The margin-parallel profile stretched for 160 km along the spine of VancouverIsland. A second profile was located perpendicular to the margin across southern VancouverIsland, extending from the volcanic front to the offshore, through the deployment of severalocean-bottom seismographs (OBS) placed on the continental shelf. Velocity structure in theupper 20 km of the crust, was well constrained. It displayed lateral homogeneity along strikethrough the Wrangellia terrane, apart from a localized discontinuity (Vp∼7.8 km/s) at 23 kmdepth, located approximately at the latitude of Campbell River within the lower crust. It wasspeculated that this feature represents a remnant of a subducted slab (McMechan and Spence,1983; Spence et al., 1985).Cassidy et al. (1998) deployed a 5 element broadband array to map crustal structure innorthern Vancouver Island using receiver functions. Their recovered S-wave velocity profilescorroborate the position for the northern limit of the Explorer plate, but contradict the hy-pothesis of a thinned crust north of Brooks Peninsula. They imaged a well-defined Moho at anormal continental depth of ∼37-39 km north of Brooks Peninsula along with an upper crustalS velocity discontinuity interpreted as the top of the underlying high velocity Wrangellia ter-614.1. IntroductionFigure 4.5: Relief and station map for northern Vancouver Island. Bathymetry and to-pograpy from ETOPO1 model (Amante and Eakins, 2009). Depth and elevationcontours are shown at 500 m intervals. Inverted triangles are stations used in thisstudy: black— CNSN stations, yellow—POLARIS arrays and red— SEAJADE IIland and OBS stations. For clarity, only CNSN stations are labelled on this map.Pink triangles are exposed Alert Bay volcanics. WB—Winona Block, EX—Explorerplate, NFZ—Nootka Fault Zone, JdF—Juan de Fuca plate.rane rocks. South of Brooks Peninsula, they imaged two dipping low velocity zones identicalto the southern Vancouver Island profiles, indicating a northward continuation in subductionprocesses.Braunmiller and Nábělek (2002) examined the source parameters of 84 large earthquakes(M>4) from the Explorer region (see figure 4.2). Their plate motion estimates suggest thatthe Winona block is part of the Explorer plate and that the latter is moving independently,624.2. Data and Methodsslowly converging with North America at a rate of 5mm/yr near the QCTJ to 20 mm/yr nearNootka Island. They observed that the easternmost seismicity along the Nootka Fault Zone,ending near Nootka Island, is deeper and with source mechanisms more similar to seismicityoff Brooks Peninsula suggestive of compression induced by the Explorer plate. Earthquakesoff Brooks Peninsula are inferred to be within the overriding North America plate and showsource mechanisms with strike-slip motion on shallow dipping, east-west striking fault planes(Splinder et al., 1997) (see figure 4.2).Audet et al. (2008) exploited the dense POLARIS-BC Northern Vancouver Island arraydeployed between June 2005 and early 2008 to image the LVZ and slab morphology using re-ceiver functions and to map upper mantle structure to 300 km depth via P-wave teleseismictomography. The sharp edge of the Explorer plate is evidenced by an abrupt disruption ofthe subducted oceanic crust (imaged as a low velocity zone, or hereafter LVZ) beneath theBPFZ and a marked contrast in upper mantle P-wave velocity at 300 km depth. They noted aflattening of the slab at crustal levels beneath the middle of array, near the town of Woss, ac-companied by a pronounced, localized anisotropy. Their proposed model explains this anomalyas a result of thermo-mechanical erosion by a toroidal mantle flow around the slab edge leadingto plate thinning at shallow levels. A resultant stretching of the Explorer slab at this locationwould explain the localized strong anisotropy signature within the slab as shear-induced. Slabcontours deepen slightly where the Nootka Fault is subducted, suggestive of a possible tear inthe slab.Royer and Bostock (2013) exploited the same POLARIS-BC dataset to detect and locateLFEs. LFE hypocenters provide a local and precise estimate of plate boundary depth, andsuggest slab flattening as proposed by Audet et al. (2008). These LFE epicenters are morelocalized along dip and along strike compared to the southern Vancouver Island LFE catalogue,and are located near the northernmost extent of automated tremor epicenters detected by Wech(2010) (see tremor contours in Figure 4.1 and 4.6).Gao et al. (2017) used these LFEs as constraints on slab geometry to calculate a 2D thermalprofile of the Explorer plate. Their model suggests that the seismogenic zone of the Explorerplate is located entirely offshore, along a potentially ∼60 km wide zone along dip with significanttsunamigenic potential, similar to the Juan de Fuca megathrust models of Hyndman and Wang(1995); Wang et al. (2003).4.2 Data and MethodsIn order to better understand how the northwest termination in subduction under VancouverIsland affects the forearc crust, we apply double-difference local earthquake tomography toderive P-wave and S-wave velocity models for the region.634.2. Data and Methods4.2.1 Data collection and selectionWe obtained the arrival time catalogue of 4607 earthquakes recorded by the CNSN stationnetwork between July 2005 and November 2014 based on picks by analysts at the GeologicalSurvey of Canada (GSC). These events are located at latitudes between 49◦N and 51◦N andlongitudes between 130◦W and 125◦W (see Figure 4.6). Approximately 35% of theses eventsare offshore earthquakes with a designated fixed (10 km) depth due to their large minimumepicentral distance and large azimuthal gap. Since these events do not afford any additionalspatial ray coverage for our area of interest compared to offshore events closer to the coast withcalculated depths, we systematically exclude them. Another 10% of the GSC catalogue eventsare blasts with an additional 3% labelled as potential blasts, that together constitute the vastmajority of catalogued events located landward of the coastline.644.2.DataandMethodsFigure 4.6: LFE hypocenters (cyan diamonds) obtained from Royer and Bostock (2013) and L.Y. Chuang (personnal communication)and event catalogue obtained from GSC. Blue dots are events with a designated fixed depth, red dots are events with a computeddepth estimate, green dots are blasts, light pink dots are designated as possible blasts. a) Map view. The CNSN, POLARIS-BCNVI and SEAJADE-II station networks are shown with the same symbols as in Figure 4.5. The perpendicular thin black linespassing through Woss are the reference axes we used to calculate along strike and along dip distances. b) Along strike depthprofile. c) Along dip depth profile. Both b) and c) have vertical exaggeration for clarity.654.2. Data and MethodsWe also acquired waveform data for stations of the POLARIS-BC NVI array that recordeddata between June 2005 and early 2008 and again briefly in late 2009-early 2010. We also ex-ploited the OBS array and land stations deployed during the SEAJADE-II experiment betweenJanuary and September 2014, and October 2013 and November 2014, respectively. Using theseadditional station networks, we attempted to detect additional small, uncatalogued earthquakesthat could potentially enhance the spatial ray sampling of the study region onshore, by modify-ing parameters for our cross-station detection algorithm as we did for our southern VancouverIsland study (Chapter 3). We detected a total of 128 new events with POLARIS data and 1110events with SEAJADE-II. In the latter instance, most of new detections are aftershocks of theMW 6.5 earthquake of April 23 2014 that occurred 20 km seaward of Nootka Island. We alsoran a standard “short-time-average/long-time-average trigger” (STA/LTA) picking algorithm(Allen, 1978; Trnkoczy, 2012) that detected 482 new events with the CNSN, POLARIS-BC andSEAJADE-II networks. Since the cross-station algorithm does not detect P-waves on the fly,the detected events require an input fixed depth to obtain epicenters and these are not expectedto be accurate since background seismicity depths vary widely. Hence, detected events requiredvisual inspection and manual picking. STA/LTA located events also required visual inspectionand manual re-picking.A total of 700 events occurring during the POLARIS-BC NVI and SEAJADE-II deploymentwere manually picked for these station networks. This set includes: CNSN catalogue eventslocated landward of the coast, the largest offshore events with unfixed depth estimates, and newuncatalogued events detected with the cross-station and STA/LTA algorithms. We attemptedto pick as many events as possible landward and far from clusters, since redundant ray pathsare not expected to significantly improve the tomographic resolution.Blasts provide an important constraint on the upper crustal velocity structure onshore.While many of these blasts are clustered around the location of known quarries, a significantnumber of blasts are likely generated in the construction of logging roads throughout the Islandand cannot be confidently associated with a known source location. We selected all events within5 km of known active quarries identified within the MINFILE database of BC Geological Surveyand visual inspection on google Earth. We then relocated these events with their identifiedsource locations and inverted for the origin time with HYPOINVERSE using the routine NVIonshore model. Events with a traveltime RMS residuals smaller than 0.25 s were retained, fora total of 66 blasts.Since onshore seismicity is sparsely distributed and Wadati-Benioff seismicity is almost non-existent landward of the coast, LFE hypocenters represent an important dataset for imagingforearc structure above the LVZ. We used the 55 LFE hypocenters of the Royer and Bostock(2013) study in addition to 68 new LFE templates obtained by Lindsay Yuling Chuang viathe application of our cross-station algorithm explained in Chapter 2 (unpublished data). Wemanually re-picked P- and S-arrival times for all LFE templates. In general, arrival time pickson LFE templates from this region are poorer quality than those for southern Vancouver Island664.2. Data and Methodsdue to the smaller number of detections used in stacking.4.2.2 Minimum 1-D Vp model inversionThe choices of initial model and initial hypocenter locations are critical to obtaining a reasonablesolution and avoiding artifacts in local earthquake tomography (e.g., Kissling et al., 1994).This sensitivity is due to the mathematical formulation of the tomography problem: a linearapproximation to the nonlinear, coupled hypocenter-velocity problem, solved via a dampedleast squares system (see Chapter 1). An iterative inversion scheme can get trapped inside alocal minimum if the input parameters are not sufficiently close to their true values. (Kisslinget al., 1994). In Chapter 3, we had access to a 3D Vp velocity model from a previous study(Ramachandran et al., 2005) that we used as initial model. However, in northern VancouverIsland, we only have access to the 2D controlled-source refraction profiles from the VISP (Spenceet al., 1985; Drew and Clowes, 1990). Offshore, an accretionary sedimentary wedge lies atop anoceanic crust with shallow Moho overlain by a thick blanket of low-velocity sediments, whereasonshore we have the high-velocity volcanic rocks of the Wrangellia terrane (Figure 4.4). As aconsequence of this contrast in velocity structure and the station event coverage (most eventsare offshore, most stations are onshore), the selection and development of an appropriate initialmodel is complicated. GSC uses two different 1D models for routine earthquake location foroffshore and onshore events, respectively (Figure 4.7a).There is a significant risk of biassing the tomographic solution by introducing a prioristructure in the initial model (e.g., Eberhart-Phillips, 1990; Kissling et al., 1994; Evans et al.,1994). In the absence of previous 3D velocity information, there are typically two approachesfor selecting an initial model. One approach is to construct an a priori model based on availablegeological data and previous 2-D refraction profiles, evaluate the fit to the traveltime data set,and then perturb the model to improve fit. A second approach is to start with a simple 1-Dmodel developed from the same traveltime dataset and allow the tomographic inversion to addperturbations only as required by the input data. Eberhart-Phillips (1990) compared bothapproaches for the Coalinga region in California, using coincident refraction 2-D profiles toconstruct an a priori model. She noted that the first approach produced an overly complexsolution for which it was difficult to assess what features were required by the data. In contrast,the second approach produced better results for which recovered structures were clearly requiredby the dataset, within the context of model regularization/parametrization, and not on anypre-conceived notions. This “minimum” 1D approach is the recommended standard, especiallyif one has little or uncertain a priori knowledge of 3D structure. (Ellsworth, 1977; Roecker,1981, 1993; Ellsworth, 1991; Kissling et al., 1994; Evans et al., 1994). However, in some cases, ifthe problem is undetermined and the ray distribution particularly inhomogeneous, an a priorimodel may produce a more credible solution (e.g., Chiarabba et al. , 1995).In the context of northern Vancouver Island, the danger of biasing the tomographic solutionis high due to the presence of major discontinuous structures in the study region (e.g., the674.2. Data and MethodsNootka Fault zone, the edge of the Explorer slab at BPFZ, underthrust terranes) and the rangeof different slab models (Audet et al., 2010; McCrory et al., 2012; Gao et al., 2017). We exploredthe option of constructing a 2D model by stitching together various 1D models, appropriatefor the OBS stations offshore (Spence et al., 1985; Obana et al., 2015) and for the onshoreWrangellia terrane (e.g. routine GSC 1D model, initial 1D model used by Ramachandran etal. (2005)). However, the solutions for 2D models constructed this way did not converge to acommon structure after multiple iterations. The chosen transition interval between the onshoreand offshore 1-D models, (i.e. the range of nodes through which we stitch the models) also has amajor effect on the solution. It is difficult to render a smooth transition due to the tomographicconstraints on grid spacing forced by the inhomogeneous data distribution. Moreover, initialhypocenters cannot be accurately calculated for such 2D models using earthquake locationalgorithms that rely on 1-D velocity models, such as HYPOINVERSE. Consequently, we electedto perform a preliminary 1D inversion using the program VELEST (Kissling, 1988; Kissling etal., 1994) to assemble an initial “minimum 1D model” that best fits our tomographic dataset.This program also computes station corrections and accounts for station elevations (as opposedto HYPOINVERSE) leading to more accurate initial hypocenters, especially since we includeboth land and OBS stations (Husen et al., 1999; Arroyo et al., 2009). It also allows us toidentify systematic errors in the dataset prior to tomographic inversion (Maurer et al., 2010).Hypocenters are more likely to be accurate when determined for a minimum 1D model thanfor a 3D model that does not represent lateral variations accurately (Theunissen et al., 2018).VELEST solves the coupled hypocenter-velocity problem by simultaneous inversion of hypocen-ter locations, 1D velocities and station corrections. It iteratively performs the full inversion viadamped least squares by penalizing the squared traveltime error and the Euclidian norm of themodel, with different damping values for hypocenter parameters, velocities and station correc-tions. In the resulting “minimum 1D model” solution, the velocity in each layer is the averagevelocity within this depth range as weighted by the ray distribution. We employ the VELESTalgorithm with a representative subset of well-constrained events (see Figure 4.9). The recov-ered station corrections reflect both near-surface geology but also large-scale deviations fromthe 1D model due to 3D structure (e.g., the subducting slab).Application of VELEST involves a trial-and-error approach as we invert for many candidate1D input models, informed in part by previous studies and available geological information butalso allowing for large model variations so as to probe the model space efficiently and identifythe global minimum. We follow the “recipe” described by (Kissling et al., 1994). We select 510well-constrained earthquakes from our dataset, including LFEs but excluding blasts, with anazimuthal gap <180◦, a minimum of 8 recorded arrivals and a minimum depth of 1 km (to reducethe risk of including blasts mislabelled as earthquakes). Kissling et al. (1994) recommends usingblasts only for the assessment of uncertainties and accuracy. Station elevations are particularlyimportant due to the use of OBS stations at depths up to about 2 km below sea level (seeFigure 4.5). However, VELEST requires all stations to be located within the first layer of the684.2. Data and Methodsmodel. Due to the large elevation difference of 4.2 km between the deepest (J03) and highest(BTB) stations, this approach would require an unreasonably thick surface layer that wouldnot accurately represent the high velocity gradients known to exist at shallow depth (Spenceet al., 1985). Hence we proceed as Husen et al. (1999) and fix all stations to zero elevation.Thus, station corrections will incorporate the effect of elevation relative to reference stationEDB (chosen due to the high number of picks at this station and its central location in theoverall station distribution). The 1D model is initially parametrized as 2 km thick layers in theupper crust and 4 km thick layers in the lower crust. The initial models we tested are shownin Figure 4.7a. For each input model, we use a simple iteration scheme, in which the solution(1D model and hypocenters) from the previous iteration is used as input for the next. Figure4.7b displays solution convergence after 100 iterations for each initial model. As expected fromthe lack of earthquakes in the shallow crust (Figure 4.7c), the convergence is not as good forthe upper 6 km as it is for the deeper layers. Converged solutions differ primarily within theshallow layers. We select for tomography the minimum 1D model that best fits the shallowvelocities from the 2D refraction profiles of Spence et al. (1985) near Woss, since most of ourpicks are at CNSN and POLARIS-BC stations, located within the Wrangellia terrane.Station corrections and relocated hypocenters determined by VELEST are shown in Figure4.8. Negative station corrections signify that real velocities are faster than in the model. OBSstations all have strongly negative corrections due to their depth relative to EDB, situatedrelatively close to sea level (163 m). Positive corrections at VI30, VI31 and VI32 could be acombined effect of low velocity Nanaimo Group sediments and 3D structure effects across theBPFZ. Overall, western land stations tend to have negative corrections, especially along theextension of the Nootka Fault and eastern land station display positive corrections, and mayreflect the effect of the subducting slab.We relocated our earthquake, LFE and blast catalogue with HYPOINVERSE, using theVELEST minimum 1D Vp model and associated P arrival station corrections, which compensatefor station elevation and depth, and a fixed Vp/Vs ratio of 1.73. This best fitting Vp/Vs ratiowas computed from a Wadati plot, i.e. by plotting the S-P arrival time difference as a functionof the P-wave traveltime, and computing Vp/Vs from the slope of the least-squares regressionfit (e.g., Chatterjee et al., 1985). HYPOINVERSE scales the S arrival corrections according tothis constant Vp/Vs ratio.After relocation, we culled the dataset to exclude poorly constrained events from tomogra-phy by removing fixed depth events (except 66 blasts identified earlier), events with convergenceproblems, less than 5 readings, an azimuthal gap larger than 300◦ or a RMS residual largerthan 0.4 s. Hypocenters determined for the resulting dataset using the minimum 1D model areshown in Figure 4.9 and compared with the locations computed with the routine GSC onshoremodel without consideration of station correction or elevations. Relocated hypocenters appearslightly more tightly clustered than those determined from the GSC model with a tendencyfor deeper events to be relocated to shallower depths. The difference in the RMS traveltime694.2. Data and MethodsVelocity (km/s)2 4 6 8-60-50-40-30-20-100# events0 10 20 30-60-50-40-30-20-100Velocity (km/s)2 4 6 8Depth (km)-60-50-40-30-20-100GSC offshore modelGSC onshore modelb) c)a)Figure 4.7: Minimum 1D Vp model. a) Different starting models used. The offshore andonshore models used for routine location by GSC are identified. Other models areconstructed from 2D refraction profiles in Spence et al. (1985) and by arbitrarilyvarying velocities to maximize sampling of the model space. b) Convergent 1D modelsobtained from a). Gray lines indicate different iterations and thick colored linesindicate final, stable solutions obtained for the different starting models. The thickblack line indicates the chosen minimum 1D model. c) Depth profile of relocatedhypocenters in the minimum 1D model.residual distribution is shown in Figure 4.10. The heavy tail of events with high residuals issignificantly reduced and the locations show an overall decrease in the average RMS residual ofabout 0.03 s.704.2. Data and MethodsFigure 4.8: Station corrections computed by VELEST for the minimum 1D model and re-located hypocenters (black dots). The color of station symbols is proportional to theP-wave station correction in seconds as shown by the colorbar (blue are negative cor-rections, red positive). The reference station is EDB. Station corrections contain theeffect of near surface geology, unmodelled 3D structure effects and station elevationrelative to EDB.714.2.DataandMethodsFigure 4.9: Tomographic dataset and relocations within the minimum 1D model (green dots) compared to locations from the routineGSC model (red dots). a) Map view, with thin lines indicates axes of along dip/strike projection. b) Along strike depth profile.c) Along dip depth profile. d) depth distribution.724.2. Data and MethodsFigure 4.10: RMS traveltime residual distribution after relocation within minimum 1Dmodel. a) Residuals for routine locations. b) Residuals for relocations with theVELEST solution.4.2.3 3D graded inversion using TomoDDUsing the minimum 1D Vp model from VELEST, we follow a graded inversion approach(Eberhart-Phillips, 1993; Foulger et al., 1995; Evans et al., 1994) to invert for 3D velocitystructure. We start the inversion for Vp and Vs on a coarse grid with a fixed Vp/Vs, andgradually increase the number of velocity nodes in order to obtain average values for poorlysampled volumes. This approach will improve convergence and simplify the interpretation ofthe final model. From the final, detailed 3-D Vp model at the finest grid spacing allowed bythe ray distribution, we invert again for Vp and Vs with a fixed Vp/Vs ratio, but with a culleddataset as we did in Chapter 4, retaining only pairs of P and S arrival times recorded at thesame station. We then divide the resulting Vp model by the Vs model to avoid issues dueto differences in ray density and sampling for P- vs S- waves (Eberhart-Phillips, 1990). SinceVp/Vs variations are small, this final inversion should recover the correct anomalies despite aconstant Vp/Vs as input (Foulger et al., 1995; Evans et al., 1994).The selection of the smoothing (regularization) parameter depends upon the grid spacing.For each inversion, we choose the smoothing parameter based on a trade-off curve (Eberhart-Phillips, 1986; Evans et al., 1994) between the weighted data variance and the model norm.The first inversion with coarsest grid has a horizontal spacing of 40 km and depth nodessituated at z = 0,6,12,18,28,38,45,55 km, selected by merging layers of similar wavespeed inthe minimum 1D model. The second inversion has the horizontal spacing reduced to 30 km734.3. Resultsand unchanged vertical spacing of 6 km. Horizontal spacing is decreased to 20 km with verticalspacing maintained at 6 km in the third inversion. The grid for the final, fourth inversion hasa horizontal spacing of 15 km and vertical spacing of 4 km.Our final dataset comprises a total of 16573 P arrival picks, 13196 S arrival picks and acatalog of 399070 and 306346 differential delay times for P- and S- waves, respectively. We donot employ waveform cross-correlation delays in this study due to time constraints.4.3 ResultsThe initial weighted RMS for the differential catalogue data for the coarser grid was 468 msand after the fourth inversion was reduced to 163 ms, a reduction of 65%. Figure 4.11 showsrelocated earthquakes and LFEs for the final model. Before interpreting these relocations andthe recovered 3-D velocity structure we present a series of resolution tests that will aid in theirinterpretation.4.3.1 Checkerboard resolution testsWe evaluate the resolution of the model by performing a synthetic checkerboard recoverytest (CRT) for two different scales of anomalies of ±10% in velocity with random Gaussian-distributed noise of 0.025 s for synthetic P traveltimes and 0.05 s for S traveltimes. Figures4.11, 4.12, 4.18 and 4.20 show the recovery of coarse anomalies of dimensions 40x40 km2 hori-zontally and 20 km vertically. Figures 4.13, 4.14, 4.15, 4.16, 4.18 and 4.19 show the recoveryof finer-scale anomalies of dimensions 20x20 km2 horizontally and 10 km vertically. Seismicityrelocated within the recovered checkerboard model is shown in the left column for reference.For depth slices, only seismicity located at or below the depth of the slice is shown. For depthprofiles, seismicity is shown within 10 km of the slice for coarse anomalies, and 5 km of the slicefor finer anomalies.The recovery of Vp and Vs is practically identical for coarse anomalies and the pattern is wellrecovered in shape and amplitude between 10 and 40 km in depth. Areally, recovery is goodbetween Brooks Peninusla and Hesquiat Peninsula, landward from the deformation front tojust seaward of the easternmost limit of the tremor contours. Recovery of finer-scale anomaliesreveals that detailed resolution is more restricted spatially. Resolution drops abruptly landwardand just below the LFE hypocenters. Depth profiles indicate that the shape of the recoveredpatterns tends to be smeared in upward and northeast direction due to most events beinglocated offshore.We compared these checkerboard patterns with the Derivative Weight Sum (DWS) values,a measured of the density of rays sampling each grid node (Thurber and Eberhart-Phillips,1999). We determine a threshold criterion of DWS/max(DWS)>0.01 that we employ in ourrecovered tomographic velocity anomaly plots to define the well resolved volume.744.3. ResultsFigure 4.11: Map view of earthquake and LFE relocated with tomoDD. A grid is shownand labelled for reference for the depth profile figures discussed in section 4.3. Theorigin is centered on the town of Woss, which is also the centre of the POLARIS-BC array. Note that the Nookta Fault Zone falls between y = −60 and −40 km.The BPFZ is located in the vicinity of x = 40 and 60 km. Red stars indicate thelocations of blasts used. Green squares are stations. Pink arrows identify lineartrends in the epicenter locations just north of the NFZ. Yellow dotted line indicatethe outline of tremor epicenters.754.3.ResultsFigure 4.12: Depth slices at 10 (top row) and 20 (bottom row) km depth for the checkerboard test with anomalies sized 40x40x20 km3.The middle column shows the recovery for Vp, and right column shows the recovery for Vs.764.3.ResultsFigure 4.13: Depth slices at 30 (top row) and 40 (bottom row) km depth for the checkerboard test with anomalies sized 40x40x20 km3.The middle column shows the recovery for Vp, and right column shows the recovery for Vs.774.3.ResultsFigure 4.14: Depth slices at 5 (top row) and 10 (bottom row) km depth for the checkerboard test with anomalies sized 20x20x10 km3.The middle column shows the recovery for Vp, and right column shows the recovery for Vs.784.3.ResultsFigure 4.15: Depth slices at 15 (top row) and 20 (bottom row) km depth for the checkerboard test with anomalies sized 20x20x10 km3.The middle column shows the recovery for Vp, and right column shows the recovery for Vs.794.3.ResultsFigure 4.16: Depth slices at 25 (top row) and 30 (bottom row) km depth for the checkerboard test with anomalies sized 20x20x10 km3.The middle column shows the recovery for Vp, and right column shows the recovery for Vs.804.3.ResultsFigure 4.17: Depth slices at 35 and 40 km depth for the checkerboard test with anomalies sized 20x20x10 km3. The middle columnshows the recovery for Vp, and right column shows the recovery for Vs.814.3.ResultsFigure 4.18: Depth profiles at y=-40,-20,0,20,40 km (from top to bottom) for the checkerboard test with anomalies sized40x40x20 km3. The middle column shows the recovery for Vp, and right column shows the recovery for Vs.824.3.ResultsFigure 4.19: Depth profiles at y=-40,-20,0,20,40 km (from top to bottom) for the checkerboard test with anomalies sized20x20x10 km3.834.3.ResultsFigure 4.20: Depth profiles at x=-80,-60,-40,-20,0 km (from top to bottom) for the checkerboard test with anomalies sized40x40x20 km3. The middle column shows the recovery for Vp, and right column shows the recovery for Vs.844.3.ResultsFigure 4.21: Depth profiles at x=-80,-60,-40,-20,0 km (from top to bottom) for the checkerboard test with anomalies sized20x20x10 km3. The middle column shows the recovery for Vp, and right column shows the recovery for Vs.854.3. Results4.3.2 Relocated SeismicityA striking feature that emerges from the tomoDD hypocenter relocations and one that is par-ticularly evident in map view is the alignment of epicenters near Nootka Island starting justnorth of y = −60. There exist at least 3 sets of parallel lineaments (labelled L1-L3 in Figure4.11) within the rectangular box defined by minimum and maximum x,y coordinates (-100,-60)and (-40,-20). It is also conceivable that similar structures exist further seaward of this region,outside of the model domain, as suggested by the fainter feature labelled L4. These lineamentstrend in a slightly more northeasterly direction than the nominal (y) strike coordinate and wespeculate that they may be the expression of strike-parallel faulting generated within/near theNFZ along the outer rise. We will refer below to the two most prominent lineaments L1 and L3in reference to velocity structure. The focal mechanism solutions of a few large (M>5) regionalearthquakes coincide spatially with the position of these lineaments, and these events exhibitstrike slip moment tensor solution with one of the focal planes having a strike parallel to thestrike of the seismicity distribution of these lineaments (see figure 4.2). The other major con-centration of seismicity, below Brooks Peninsula is more limited in its along-strike content withlittle or faint indication of internal lineation, albeit locations are less accurate due to a moresparse distribution of stations in that region. With the exception of minor possible clusteringat the northern end of Vancouver Island, most remaining onshore seismicity is very sparselydistributed.4.3.3 Vp, Vs and Poisson’s ratio modelsNootka Fault Zone VicinityFigure 4.22 displays several along-dip depth profiles in the vicinity of the NFZ. The Vp and Vsvelocity contours appear to flatten to the northwest in agreement with the slab model of Audetet al. (2010). Just south of the NFZ, at y = − 65 km, a thick, high Poisson’s ratio anomalyappears near the top of the Juan de Fuca slab that may manifest the presence of overpressuredfluids in an oceanic crustal LVZ. This signature is disturbed to the northwest within the NFZ.Lineaments L1 and L3 are more or less evident within the seismicity distribution. The strongestconcentration of seismicity in L1 appears to favour regions of high Vp and Poisson’s ratio, andto shift to shallower levels to the northwest, although this tendency is not in evidence for L3.The complex co-variation in velocity represented by the Poisson’s ratio image may represent in-terplay between multiple mechanisms including strong plate bending and dehydration of mantleserpentinite that was preferrably formed along the NFZ seward of the trench. More generally,seismicity near x = −80 km, extends between 10 and 40 km depth, and most would appear toreside within the oceanic crust and mantle with more minor or absent levels in the deep forearccrust, dependent upon which slab model (Audet et al. (2010) versus McCrory et al. (2012) is amore accurate portrayal of true structure. The validity of either model in the offshore regionmay brought into question given the presence of the NFZ, a recently formed plate boundary.864.3. ResultsLFE vicinityFigure 4.23 displays structure along dip profiles that transect the LFEs, near the middle ofthe Explorer plate. A flattening of the slab contours is again reflected in the Vp and Vscontours, in reasonable agreement with the Audet et al. (2010) slab model. LFE relocationsat y =20 km display more scatter than those at y =−5 km, perhaps due to large uncertaintiesfor some template picks. Those LFEs along the y =−5 km profile are well clustered andconsistent with flattening of the Audet et al. (2010) slab interface although they fall a fewkm deeper. In both depth profiles containing LFEs, Poisson’s ratio images reveal associatedanomalies of up to 0.29 that, as in southern Vancouver Island, likely manifest the ppresenceof high pore fluid pressures within the oceanic crust. This anomaly is confined to the flatslab region disappearing updip where regular earthquakes occur. The appearance of a lowPoisson’s ratio (0.20-0.22) anomaly within the crust immediately landward of LFE terminationin profiles y = −5,20 km is interesting in the context of the observations and interpretationof a similar feature in southern Vancouver Island described in Chapter 3. It is reasonable tospeculate that slab fluids and metasomatism are also at play here; however, this feature is notassociated with crustal microseismicity and the density mapped surface faults is relatively low.Ray sampling of this feature is likely to be directionally biased and limited to LFE sources, sosome caution in its interpretation is warranted. The profile at y =10 km lies just beyond thenorthwestern edge of LFE hypocenters, where the LVZ is no longer in evidence although thismay reflect a decrease in ray coverage since recovery of checkerboard patterns near the Audetet al. (2010) slab model is marginal away from the LFE source region. Tremor epicenters fromthe automated detection algorithm of Wech (2010) do extend farther to the northwest, but theyare not manually reviewed and, due to their fixed depth and absence of clear onsets, generallyhave larger horizontal location errors. The high Vp anomaly of 6.8-7.0 km/s near 10 km depthbetween -60 and -80 km in x present on depth profiles at y=-5 and 10 km may represent thefragment of Crescent terange identified by Johns et al. (2012).Brooks Peninsula VicinityVelocity structure near the BPFZ is shown in Figure 4.24 for values of y=25,45,55 km. Thisregion likely encompasses the edge of the subducting Explorer plate and the transition to a moretypical crustal structure (Cassidy et al., 1998; Audet et al., 2008). At y=25 km, just south ofthe BPFZ, velocity contours and display some dip between x=-60 and x=0 km and remain inrough agreement with both plate models suggesting northward persistance of the slab edge.Clustered seismicity between 20-30 km offshore lies approximately 10-15 km south of the BPFZand much of it appears to reside within the forearc crust, 1-6 km seaward of the Cape Cookfault. Most northwesterly trending faults (e.g., the West Coast fault) like the Cape Cook faultare expected to dip northward so the origin of this seismicity is unclear. In contrast, seismicityat y=40 km is deeper seated and appears to straddle the crust mantle boundary in its seawardportion. Localized depressions in Vp and Vs contours at the landward edge of this seismicity874.3. Resultsmay signal partial underthrusting of the southernmost Winona block or a deep expression of itsboundary with the Explorer plate. Concentrated seismicity ceases at the landward initiationof Brooks Peninsula and coincides with high Poisson’s ratios of 0.27-0.29. North of the BPFZ,for example along dip profile y=55 km, both Vp and Vs values are rather lower than they areto the south, consistent with Lewis et al. (1988) observations of higher heat flow in this region.North and the BPFZ and landward (x >−40 km), where our model is well resolved, weobserved lower Vp and Vs, consistent with higher heat flow measurements.Along strike trendsFigure 4.25 displays five along-strike profiles at x=-60,-45,-30,-15, 0 km. One feature of interestis the presence of 2 clusters of crustal seismicity approximately 20 km deep in northermostVancouver Island, north of the plate edge in profiles at x=-60,-45,-30 km at an alongstrikeposition of y=90 km. In map view (see Figures 4.11 and 4.26a), they are evident just west ofstation HOLB and located within 6 km northwest of the Holdberg Fault (Figure 4.26a). Thesetwo clusters contain 6 and 10 events. A dip profile oriented across the fault (Figures 4.26b and4.27) reveal a linear trend for the bigger cluster, indicative of a plane dipping at ∼70 degrees thatwould intersect the surface near the trace of the Holberg Fault. The depths of these events issurprising given that the maximum depth of crustal seismicity as defined by the brittle-ductiletransition is inversely proportional to heat flow (e.g., Hauksson, et al., 2011). Along-strikeprofiles in Figure 4.25 and along dip profile in Figure 4.27 reveal lower Poisson’s ratios in theupper crust in the vicinity of the Holberg fault. It may also be of interest to note that theHolberg Fault hosts significant Cu, Au, Mo, Ag and perlite deposits, all classified as being ofporphyry-copper/hydrothermal origin (Figure 4.26a). There also appear to be some seismicitythat could be associated to the Cape Cook Fault trace if this fault is dipping northeastward at30 degrees, the expected dip angle for a splay fault (Figure 4.27).Most of the velocity profiles in Figure 4.25 indicate a decrease in crustal velocity occurringnear y=40-50 km approximately in line with the BPFZ, indicative of the change in heat flowand tectonic environment as mentioned above. The slab LVZ is evident only where LFEs weredetected, although the localized distribution of LFE templates in the vicinity of the POLARIS-BC array probably limit its detection to the south and east.884.3.ResultsFigure 4.22: Along dip depth profiles at y=-65,-50,-35 km, from top to bottom, near the Nootka Fault Zone. The first column showVp, the second Vs and the last Poisson’s ratio. The white dashed contour line is the extent of well-resolved areas, defined asDWS/max(DWS) > 0.01. The solid gray line shows the plate interface model of Audet et al. (2010) and the dotted grayline, the McCrory et al. (2012) plate interface model. Inverted yellow triangles are stations within ± 10 km. Seismicity isshown within ± 7 km. The y-axis is depth, and the x-axis represents the along dip distance, with zero corresponding to thetown of Woss. See 4.11 for a map view of the reference axes.894.3.ResultsFigure 4.23: Same as 4.22, but for along dip depth profiles at y=-20,-5,10 km, from top to bottom, near LFE seismicity and themiddle of the Explorer plate. The y-axis is depth, and the x-axis represents the along dip distance, with zero correspondingto the town of Woss. See 4.11 for a map view of the reference axes.904.3.ResultsFigure 4.24: Same as 4.22, but for along dip depth profiles at y=25,40,55 km, from top to bottom, near the Brooks Peninsula FaultZone. The first column show Vp, the second Vs and the last Poisson’s ratio. The y-axis is depth, and the x-axis representsthe along dip distance, with zero corresponding to the town of Woss. See 4.11 for a map view of the reference axes.914.3.ResultsFigure 4.25: Same as 4.22, but for along strike depth profiles at x=-60,-45,-30,-15,0 km, from top to bottom. The first column showVp, the second Vs and the last Poisson’s ratio. The y-axis is depth, and the x-axis represents the along strike distance, withzero corresponding to the town of Woss. See 4.11 for a map view of the reference axes.924.3. ResultsFigure 4.26: Close-up of the area around the Holberg Fault. a) Tectonic context: thickblue lines are named faults in the area and thinner blue lines are secondary faults.Red circles are known and developed mineral deposits. Black circles is relocatedseismicity. Thick black line is the depth section across the Holberg Fault and shownin b. Red star is the location of quarry blast events used in this study. b) Depthsection perpendicular to the Holberg Fault. BIF—Browning Inlet Fault, GRF—Goodspeed River Fault, NRF—Nahwitti River Fault, NIF—Neroutsos Inlet Fault,SCF—Spencer Cove Fault, SRF—Stranby River Fault, WLF—William Lake Fault,RD—Red Dog Cu deposit, H—Hushamu Cu-Au deposit, AB—Apple Bay silica(perlite) quarry, ICM—Island Copper Mine.934.4. Concluding RemarksFigure 4.27: Velocity structure near the Holberg Fault at y=70 km. Relocated seismicitywithin 30 km north or 10 km south of the depth section is included to representall events in the Holberg Fault region. Deformation front and intersecting faulttraces labelled as inverted black triangles. Inverted yellow triangles indicate nearbystations. Black dashed lines indicate linear trends in the seismicity and possiblefault planes. White dashed contours indicate the well-resolved area of the model.CCF—Cape Cook fault, WCF—Westcoast fault, HF—Holberg fault4.4 Concluding RemarksWe have used microseismicity, quarry blasts and LFEs to develop the first 3-D crustal scalemodel of Vp, Vs and Poisson’s ratio for northern Vancouver Island. Although resolution is944.4. Concluding Remarksrather more limited than it is for southern Vancouver Island, a number of interesting featureshave been documented that warrant further investigation. These include:• Distinct, parallel lineaments that trend slightly oblique to strike are resolved within seis-micity concentrations of the NFZ immediately offshore Nootka Island. We speculate thatthese features may represent faulting developed as result of high strains associated withslab bending, amplified in the vicinity of the newly formed plate boundary.• Velocity contours near the top of the Explorer plate display a flattening near the town ofWoss, corroborating the documentation of a flat slab made by Audet et al. (2010) usingreceiver functions.• LFEs below northern Vancouver Island coincide with an LVZ characterized by Poisson’sratios near 0.27-0.29, consistent with strong fluid overpressures within the Explorer plateoceanic crust.• A low Poisson’s ratio anomaly of 0.20-0.22 within the forearc crust is located just landwardof the LVZ and LFE hypocenters, a configuration which closely resembles that of southernVancouver Island. However, this feature is not accompanied by significant microseismicity.• Deep (∼20 km) crustal seismicity along the Holberg Fault north of the Explorer Plateedge that suggests that this structure is deep seated and seismogenic.95Chapter 5ConclusionsIn this section, we list the main contributions and limitations of the research described inChapters 2, 3 and 4. We also list revelant future research directions for each project.5.1 Detection and Location of Low Frequency Earthquakes UsingCross-station CorrelationThe automated cross-station algorithm we developed provided a fast and efficient way of ob-taining individual LFE waveforms suitable for the construction of LFE templates. It has beenused to extend the existing LFE catalogue in southern and northern Vancouver Island, andalso allowed the identification of the first reported LFEs in Alaska (Chuang et al., 2017) andnorthern California (Plourde et al., 2015). The more extensive LFE catalogue for southernVancouver Island has proven useful for the study of LFE magnitudes and their spatiotemporaldistribution during Rapid Tremor Reversals (Bostock et al., 2015; Royer et al., 2015). TheLFEs have also provided insights into the spatiotemporal slip migration patterns during andin-between ETS episodes and have provided useful metrics that can be exploited to developtheoretical and conceptual source models of slow slip (e.g., Hawthorne et al., 2016). We havealso demonstrated the flexibility of the method by successfully employing it in the detection ofuncatalogued microseismicity.One limitation is that the method currently detects only one wavetype, assuming that theother is of negligible amplitude in comparison. If this assumption is violated and if waveformsimilarity is not highly consistent due to simple source-time functions such as in the case ofLFEs in southern Vancouver Island, then the on-the-fly epicentral estimates have a large error.Extending the algorithm with an automated detection of the second type of wavetype wouldhelp allow hypocentral (versus epicentral) location and improved application to the detectionof micro-seismicity.At present, detected microseismic events require visual inspection due to the high rate offalse positives. This is due to the difficulties in treating glitches and other high amplitude noisesignals via an automated workflow. In order to make our algorithm suitable for stable, long965.2. Seismicity, Metamorphism, and Fluid Evolution Across the Northern Cascadia Forearcterm monitoring, we would need to develop a robust filtering scheme that can classify detectionsas LFE, micro-seismicity, blast or noise. Supervised machine learning methods have recentlyproven to be useful and efficient to solve this type of seismic classification problem (Reynenand Audet, 2017; Paitz et al., 2018).5.2 Seismicity, Metamorphism, and Fluid Evolution Across theNorthern Cascadia ForearcOur Vp, Vs and Poisson’s ratio models of southern Vancouver Island have significantly improvedthe resolution near and above the LFE hypocenters as compared to the previous 3D model ofRamachandran and Hyndman (2012). Our results allowed us to build a conceptual model of therole of slab-derived fluids in LFE generation and forearc microseismicity. It will be interestingto see if observations of Vp and Poisson’s ratio anomalies at other warm subductions zonecorroborate our findings.One important question raised in this study concerns the source processes at work for theanomalous Texada cluster of deep Wadati-Benioff earthquakes. An analysis of focal mechanismsof these events would provide additional insight as to whether a propagator wake is respon-sible for this cluster. In addition, it will be important to ascertain whether these events aredominantly located within the oceanic crust or mantle.The inferred role of fluids in triggering deep forearc seismicity also raises a question aboutwhether the source mechanisms of these earthquakes are accompanied by a non-double compo-nent indicative of collapse or expansion of the fault plane that would corroborate the presenceof high pore fluid pressure within the fault zone. Because of the low magnitude of the eventscorresponding to these clusters, the identification of first motion polarities on seismograms ischallenging. This challenge may be addressed through development of an algorithm that ap-plies waveform cross-correlation to the inversion of moment tensor solutions. My colleagueAlexandre Plourde is currently developing such an algorithm that shows promising preliminaryresults.5.3 Double-difference Tomography of Northern Vancouver IslandOur research on northern Vancouver Island is a significant contribution since it representsthe first 3D Vp, Vs and Poisson’s ratio models to date for the region. While the resolution weachieved in this region is limited by the lack of crustal seismicity within the North America plate,we have nonetheless identified several interesting features that warrant further investigation,such as large oblique faults that appear within the subducting Explorer slab close to the coastand which correlate with both high and low Poisson’s ratio anomalies. As we prepare theresults of this study for publication, we will conduct further investigation in the coming monthsto find an explanation to these anomaly patterns. We have also identified a potentially activeseismogenic fault north of the Explorer plate boundary that is of relevance to local earthquake975.3. Double-difference Tomography of Northern Vancouver Islandhazard assessment.98BibliographyAki, K., and W.H.K. Lee (1976). Determination of Three-dimensional Anomalies under a Seis-mic Array Using First P Arrival Times from Local Earthquakes. 1. A Homogeneous InitialModel, J. Geophys. Res., 81,4381-4399.Abers, G. A., MacKenzie, L. S., Rondenay, S., Zhang, Z., Wech, A. G. and K. C. Creager,(2009). Imaging the source region of Cascadia tremor and intermediate-depth earthquakes,Geology, 37, 1119-1122.Abers, G. A., Nakajima, J., van Keken, P. E., Kita, S., and Hacker, B. R. (2013). Thermal-petrological controls on the location of earthquakes within subducting plates. Earth andPlanetary Science Letters, 369–370, 178–187., A. C. and G. C. Beroza (2014). PageRank for Earthquakes, Seismol. Res. Lett., 85, 2,344-350.Ahrens, T. J. and Schubert, G. (1975). Gabbro-eclogite reaction rate and its geophysical sig-nificance. Reviews of Geophysics Space Physics, 13, 383-400.Allen, R. V. (1978). Automatic earthquake recognition and timing from single traces.Bulletin of the Seismological Society of America, 68(5), 1521–1532. Retrieved from, C. and B. W. Eakins, ETOPO1 1 Arc-Minute Global Relief Model: Procedures, DataSources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24, 19 pp, March2009., J., Kim, W., and A. M. Rubin, (2014). Accurate tremor locations from coherentS and P waves, J. Geophys. Res. Solid Earth, 119, 6, 5000-5013.Armstrong, R. L., Muller, J. E., Harakal, J. E., and Muehlenbachs, K. (1985). The NeogeneAlert Bay Volcanic Belt of northern Vancouver Island, Canada: Descending-plate-edge vol-canism in the arc-trench gap. Journal of Volcanology and Geothermal Research, 26(1–2),75–97., I. G., Husen, S., Flueh, E. R., Gossler, J., Kissling, E., and Alvarado, G. E.(2009). Three-dimensional P-wave velocity structure on the shallow part of the CentralCosta Rican Pacific margin from local earthquake tomography using off- and onshore net-works. Geophysical Journal International, 179(2), 827–849., B.F., Musumi-Rokkaku, S., Satake, K., Tsuji, Y., Ueda, K., and Yamaguchi, D.K.,2015, The orphan tsunami of 1700—Japanese clues to a parent earthquake in North America,2nd ed.: Seattle, University of Washington Press, U.S. Geological Survey Professional Paper1707, 135 p.Audet, P., Bostock, M. G., Mercier, J.-P., and Cassidy, J. F. (2008). Morphology of theExplorer–Juan de Fuca slab edge in northern Cascadia: Imaging plate capture at a ridge-trench-transform triple junction. Geology, 36(11), 895., P., Bostock, M. G., Boyarko, D. C., Brudzinski, M. R. and R. M. Allen, (2010). Slabmorphology in the Cascadia forearc and its relation to episodic tremor and slip, J. Geophys.Res. 115, B00A16.Audet, P., Bostock, M. G., Christensen, N. I., and Peacock, S. M. (2009). Seismic evidencefor overpressured subducted oceanic crust and megathrust fault sealing. Nature, 457, 76–78., N. J., Cassidy, J. F., Dosso, S. E., and Mazotti, S. (2011). Mapping crustal stressand strain in southwest British Columbia, Journal of Geophysical Research, 116, B03314,doi:10.1029/2010JB008003.Balfour, N. J., Cassidy, J. F., and Dosso, S. E. (2012). 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., G. C. and S. Ide, (2011). Slow earthquakes and nonvolcanic tremor, Annu. Rev. EarthPlanet Sci., 39, 1, 271-296.Booth, D. C. and S. Crampin, (1985). Shear wave polarization on a curved wavefront at anisotropic free surface, Geophys. J. R. Astron. Soc., 83, 31-45.Bostock, M. G., Thomas, A. M., Savard, G., Chuang, L., and Rubin, A. M. (2015).Magnitudes and moment-duration scaling of low-frequency earthquakes beneath south-ern Vancouver Island. Journal of Geophysical Research: Solid Earth, 120(9), 6329–6350., M. G., Hyndman, R. D., Rondenay, S., and Peacock, S. M. (2002). An in-verted continental Moho and serpentinization of the forearc mantle. Nature, 417, 536–538., M. G. (2012). The Moho in subduction zones, Tectonophysics, 609, 547-557.Bostock, M. G. and N. I. Christensen, (2012). Split from slip and schist: Crustal anisotropybeneath northern Cascadia from non-volcanic tremor, J. Geophys. Res. 117, B08303.Bostock, M. G., Royer, A. A., Hearn, E. H. and S. M. Peacock, (2012). Low frequency earth-quakes below southern Vancouver Island, G-cubed, 13, 11, 11007.Botev, Z. I., Grotowski, J. F., and D. P. Kroese, (2010). Kernel density estimation via diffusion,Ann. Stat., 38, 5, 2916-2957.100BIBLIOGRAPHYBoyarko, D. C., and Brudzinski, M. R. (2010). Spatial and temporal patterns of non-volcanictremor along the southern Cascadia subduction zone. Journal of Geophysical Research, 115,B00A22, doi:10.1029/2008JB006064.Braunmiller, J., and Nábělek, J. (2002). Seismotectonics of the Explorer region.Journal of Geophysical Research: Solid Earth, 107(B10), ETG 1-1-ETG 1-25., C. M. and Ague, J. J. (2002). Slab-derived fluids and quartz-vein formation in anaccretionary prism, Otago Schist, New Zealand. Geology, 30, 499-502.Brocher, T. M., et al.(2001), Upper crustal structure in Puget Sound, Washington: Resultsfrom 1998 Seismic Hazards Investigation of Puget Sound. Journal of Geophysical Research,106, 13,541–13,564.Brocher, T. M., Parson, T., Trehu, A. M., Snelson, C. M. and Fischer, M. A. (2003). Seis-mic evidence for widespread serpentinized forearc upper mantle along the Cascadia margin.Geology, 31, 267-270Brown, J. R., Beroza, G. C. and D. R. Shelly, (2008). An autocorrelation method to detect lowfrequency earthquakes within tremor, Geophys. Res. Lett., 35, 16, L16305.Brudzinski, M. R. and R. M. Allen, (2007). Segmentation in episodic tremor and slip all alongCascadia, Geology, 35, 10, 907.Calkins, J. A., Abers, G. A., Ekström, G., Creager, K. C. and S. Rondenay, (2011). Shallowstructure of the Cascadia subduction zone beneath western Washington from spectral ambientnoise correlation, J. Geophys. Res., 116, B07302.Canales, J., Carbotte, S. M., Nedimovic, M. R., and Carton, H. (2017). Dry Juan de Fuca slabrevealed by quantification of water entering Cascadia subduction zone, Nature Geoscience,10, 864-870, http://doi: 10.1038/NGEO3050Cassidy, J. F., Ellis, R. M., and Rogers, G. C. (1988). The 1918 and 1957 Vancouver Islandearthquakes. Bull. Seismol. Soc. Am., 78(2), 617–635.Cassidy, J. F., Ellis, R. M., Karavas, C., and Rogers, G. C. (1998). The northern limit ofthe subducted Juan de Fuca Plate system. Journal of Geophysical Research: Solid Earth,103(B11), 26949–26961., J. S. and T. I. Melbourne, (2009). Future Cascadia megathrust rupture delineatedby episodic tremor and slip, Geophys. Res. Lett., 36, 22, L22301.Chatterjee, S.N., Pitt, A.M., Iyer, H.M., (1985). Vp/Vs ratios in the Yellowstone National Parkregion, Wyoming. Journal of Volcanology and Geothermal Research 26, 213–230.Chiarabba, C., Amato, A., and Evans, J. R. (1995). Variations on the NeHT high-resolution tomography method: a test of technique and results for Medicine Lakevolcano, northern California. Journal of Geophysical Research, 100(B3), 4035–4052., N. I. (1984), Pore pressure and oceanic crustal seismic structure, GeophysicalJounal of the Royal Astronomical Society, 79, 411-424.101BIBLIOGRAPHYChristensen, N. I. (1996). Poisson’s ratio and crustal seismology. Journal of Geophysical Re-search, 101, 3139-3156.Chuang, L., Bostock, M., Wech, A., and Plourde, A. (2017). Plateau subduction,intraslab seismicity, and the Denali (Alaska) volcanic gap. Geology, 45, 647–, R.M., Zelt, C. A., Amor, J. A. and Ellis, R. M. (1995). Lithospheric structure in thesouthern Canadian Cordillera from a network of seismic refraction lines. Canadian Journalof Earth Sciences, 32, 1485-1513.Cui, Y., Miller, D., Schiarizza, P., and Diakow, L.J., 2017. British Columbia digital geology.British Columbia Ministry of Energy, Mines and Petroleum Resources, British ColumbiaGeological Survey Open File 2017-8, 9p.Dahm, T., Fischer, T. (2014). Velocity ratio variations in the source region of earthquakeswarms in NW Bohemia obtained from arrival time double differences. Geophysical JournalInternational, 196, 2, pp. 957—970., E. E., and Riddihough, R. P. (1982). The Winona Basin: structure and tectonics.Canadian Journal of Earth Sciences, 19(4), 767–788., E. E., and R. D. Hyndman (1989), Accretion and recent deformation of sedimentsalong the northern Cascadia subduction zone, Geol. Soc. Am. Bull., 101, 1465–1480,doi:10.1130/0016-7606(1989)101<1465:AARDOS>2.3.CO;2.Dragert, H., Wang, K., and James, T. S. (2001). A silent slip event on the deeper Cascadiasubduction interface. Science, 292, 1525-1528.Dragert, H., Wang, K. and G. Rogers, (2004). Geodetic and seismic signatures of episodic tremorand slip in the northern Cascadia subduction zone, Earth Planets Space, 56, 1143-1150.Drew, J. J., and Clowes, R. M. (1990). A re-interpretation of the seismic structure across theactive subduction zone of Western Canada in studies of laterally heterogeneous structuresusing seismic refraction and reflection data, edited by A.G. Green. Geol. Surv. Can. Pap.,89–13, 115–132.Dziak, R. P. (2006). Explorer deformation zone: Evidence of a large shear zone and reorgani-zation of the Pacific-Juan de Fuca-North American triple junction. Geology, 34(3), 213–216., A. M., T.-A. Chou and J. H. Woodhouse, Determination of earthquake sourceparameters from waveform data for studies of global and regional seismicity, J. Geophys.Res., 86, 2825-2852, 1981. doi:10.1029/JB086iB04p02825Earle, S. (2016). Physical Geology. Vancouver, Bristish Columbia: CreateSpace IndependentPublishing Platform. (Retrieved online, D. (1986). Three-dimensional velocity structure in northern California CoastRanges from inversion of local earthquakes arrival times. Bull. Seismol.Soc. Am. 76 : 1 025-32Eberhart-Phillips, D. (1990). Three-dimensional P and S velocity structure in theCoalinga region, California. Journal of Geophysical Research, 95(90), 15343–15363., D., Local earthquake tomography: earthquake source regions, in: SeismicTomography: Theory and Practice, eds. H. M. Iyer and K. Hirahara, Chapman and Hall,London, 613−643, 1993.Ekström, G., M. Nettles, and A. M. Dziewonski, The global CMT project 2004-2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth Planet. Inter., 200-201, 1-9, 2012.doi:10.1016/j.pepi.2012.04.002Ellsworth, W. L., Three-dimensional structure of the island of Hawaii, Ph.D. Thesis, Mas-sachusetts Institute of Technology, Cambridge, MA, 327 p., 1977.Ellsworth, W. L., D. Eberhart-Phillips and E. Kissling, A test of local earthquake tomography,Seis. Res. Lett., 62, 31, 1991.Evans, J. R., Eberhart-Phillips, D., and Thurber, C. H. (1994). User’s manual for SIMULPS12 for imaging Vp and Vp/Vs: A derivative of the Thurber tomographic in version SIMUL3for local earthquakes and explosion. Open-File Report 94-431, 101.Fairchild, L. H., and Cowan, D. S. (1982). Structure, petrology, and tectonic history of theLeech River complex northwest of Victoria, Vancouver Island. Canadian Journal of EarthSciences, 19, 1817-1835.Fisher, D. M., Brantley, S. L., Everett, M., and Dzvonik, J. (1995). Cyclic fluid flow througha regionally extensive fracture network within the Kodiak accretionary prism, Journal ofGeophysical Research, 100, 12881-12894.Fisher, M. A., et al.(1999). Seismic survey probes urban earthquake hazards in Pacific North-west. Eos, Transactions of the AGU, 80(2), 13–17.Foulger, G. R., Miller, A. D., Julian, B. R., and Evans, J. R. (1995). Three‐dimensionalvp and vp/vs structure of the Hengill Triple Junction and Geothermal Area, Ice-land, and the repeatability of tomographic inversion. Geophysical Research Letters., C.G., Dziak, R.P., Matsumoto, H., Schreiner, A.E., 1994. Potential for monitoring low-level seismicity on the Juan de Fuca Ridge using fixed hydrophone arrays. Mar. Tech. Soc.J. 27, 22–30.Frank, W. B. and N. M. Shapiro, (2014). Automatic detection of low-frequency earthquakes(LFEs) based on a beamformed network response, Geophys. J. Int., 197, 2, 1215-1223.Fréchet, J. (1985). Sismogenèse et doublets sismiques, Thèse d’État, Université Scientifique etMédicale de Grenoble, 206 pp.Gao, D., K. Wang, E. E. Davis, Y. Jiang, T. L. Insua, and J. He (2017), Thermal state ofthe Explorer segment of the Cascadia subduction zone: Implications for seismic and tsunamihazards, Geochem. Geophys. Geosyst., 18, doi:10.1002/2017GC006838Geiger, L. (1912). Probability method for the determination of earthquake epicenters from thearrival time only, Bull. St. Louis Univ. 8, 60-71.Ghosh, A., Vidale, J. E., Sweet, J., Creager, K. C. and A. Wech, (2009). Tremor patches inCascadia revealed by seismic array analysis, Geophys. Res. Lett., 36, L17316.103BIBLIOGRAPHYGhosh, A., Vidale, J. E., Sweet, J., Creager, K. C., Wech, A., Houston, H. and E. E. Brodsky,(2010). Rapid, continuous streaking of tremor in Cascadia, G-cubed, 11, 12, 12010.Ghosh, A., Vidale, J. E. and K. C. Creager, (2012). Tremor asperities in the transition zonecontrol evolution of slow earthquakes, J. Geophys. Res., 117(B10):B10301.Gibbons, S. J. and F. Ringdal, (2006). The detection of low magnitude seismic events usingarray-based waveform correlation, Geophys. J. Int., 165, 149-166.Goldfarb, R. J., and Groves, D. I. (2015). Orogenic gold: Common or evolving fluid and metalsources through time. Lithos, 233, 2-26.Goldfinger, C., Ikeda, Y., Yeats, R.S., and Ren, J., 2013, Superquakes and Supercycles, Seis-mological Research Letters, v. 84, no. 1, p. 24-32Greene, A. R., Scoates, J. S., Weis, D., Katvala, E. C., Israel, S., and Nixon, G. T. (2010).The architecture of oceanic plateaus revealed by the volcanic stratigraphy of the accretedWrangellia oceanic plateau. Geosphere, 6(1), 47–73., B. (1996). Eclogite formation and the rheology, buoyancy, seismicity, and H2O contentof oceanic crust, in Subduction: Top to Bottom, Geophys. Monogr. Ser., vol. 96, edited byG. E. Bebout et al., pp. 337-346, AGU, Washington, D.C., doi:10.1029/GM096p0337Hansen, R. T. J., Bostock, M. G. and N. I. Christensen, (2012). Nature of the low velocity zonein Cascadia from receiver function waveform inversion, Earth Planet. Sci. Lett., 337-338,25-38.Hauksson, E. (2011). Crustal geophysics and seismicity in southern California. GeophysicalJournal International, 186(1), 82–98., A. and Nakajima, J. (2017). Seismic imaging of slab metamorphism and genesisof intermediate-depth intraslab earthquakes. Progress in Earth and Planetary Science, 4:12,1-31., J. C., M. G. Bostock, A. A. Royer, and A. M. Thomas (2016), Variations in slowslip moment rate associated with rapid tremor reversals in Cascadia, Geochem. Geophys.Geosyst., 17, 4899–4919, doi: 10.1002/2016GC006489.Horning, G., Canales, J. P., Carbotte, S. M., Han, S., Carton, H., Nedimović, M. R., andvan Keken, P. E. (2016). A 2-D tomographic model of the Juan de Fuca plate from accre-tion at axial seamount to subduction at the Cascadia margin from an active source oceanbottom seismometer survey. Journal of Geophysical Research: Solid Earth, 121, 5859–5879., H. (2015). Low friction and fault weakening revealed by rising sensitivity of tremorto tidal stress. Nature Geoscience, 8, 409–416., H., Delbridge, B. G., Wech, A. G. and K. C. Creager, (2011). Rapid tremor reversalsin Cascadia generated by a weakened plate interface, Nature Geoscience, 4, 6, 404-409.Husen, S., Kissling, E., Flueh, E. R., and Asch, G. (1999). Accurate hypocenter determinationin the shallow part of the Nazca subduction one in Northern Chile using a combined on-/offshore network. Geophysics Journal Internal, 138(May), 687–701.104BIBLIOGRAPHYHusen, S., Kissling, E., and Clinton, J. F. (2011). Local and regional minimum 1D models forearthquake location and data quality assessment in complex tectonic regions: Application toSwitzerland. Swiss Journal of Geosciences, 104(3), 455–469., R. D., Riddihough, R. P., and Herzer, R. (1979). The Nootka Fault Zone - anew plate boundary off western Canada. Geophysical Journal International, 58(3), 667–683., R.D. (1988). Dipping seismic reflectors, electrically conductive zones and trappedwater in the crust over a conducting plate. Journal of Geophysical Research, 93, 13391–13405.Hyndman, R. D., McCrory, P. A., Wech, A. G., Kao, H., and Ague, J. (2015). Cascadia sub-ducting plate fluids channelled to forearc mantle corner: ETS and silica deposition. Journalof Geophysical Research: Solid Earth, 120, 4344-4358., R. D., Yorath, C. J., Clowes, R. M., and Davis, E. E. (1990). The northern Casca-dia subduction zone at Vancouver Island: seismic structure and tectonic history. CanadianJournal of Earth Sciences, 27, 313-329.Hyndman, R. D. and K. Wang, (1993). Thermal constraints on the zone of major thrust earth-quake failure: The Cascadia subduction zone, J. Geophys. Res., 98, 2039-2060.Hyndman, R. D., and K. Wang (1995), The rupture zone of Cascadia great earthquakesfrom current deformation and the thermal regime, J. Geophys. Res., 100, 22,133–22,154,doi:10.1029/95JB01970.Ide, S., Shelly, D. R. and G. C. Beroza, (2007). Mechanism of deep low frequency earthquakes:Further evidence that deep non-volcanic tremor is generated by shear slip on the plate inter-face, Geophys. Res. Lett., 34, L03308.Ide, S., (2010). Striations, duration, migration and tidal response in deep tremor, Nature, 466,356-359.International Seismological Centre, On-line Bulletin,, Internatl. Seismol.Cent., Thatcham, United Kingdom, 2015.Ito, K. (1999). Seismogenic layer, reflective lower crust, surface heat flow and large inlandearthquakes. Tectonophysics, 306, 423-433.Janiszewski, H. A., and Abers, G. A. (2015). Imaging the plate interface in the Cascadiaseismogenic zone: new constraints from offshore receiver functions. Seismological ResearchLetters, 86, 1261-1269.Johns, M. J., Trotter, J. A., Barnes, C. R., Narayan, Y. R., and Dix, G. R. (2012). Bios-tratigraphic, strontium isotopic, and geologic constraints on the landward movement andfragmentation of terranes within the Tofino Basin, British Columbia. Canadian Journal ofEarth Sciences, 49(7), 819–856., D.L., Silberling, N.J., and Hillhouse, J., 1977, Wrangellia—A displaced terrane in north-western North America: Canadian Journal of Earth Sciences, v. 14, p. 2565–2577, doi:10.1139/e77-222.105BIBLIOGRAPHYKao, H., Shan, S. J., Dragert, H., Rogers, G., Cassidy, J. F. and K. Ramachandran, (2005). Awide depth distribution of seismic tremors along the northern Cascadia margin, Nature, 436,841-844.Kao, H., Shan, S., Dragert, H., Rogers, G., Cassidy, J. F., Wang, K., James, T.S. and K.Ramachandran, (2006). Spatial-temporal patterns of seismic tremors in northern Cascadia,J. Geophys. Res., 111(B3):B03309.Kao, H., Wang, K., Chen, R. Y., Wada, I., He, J., and Malone, S. D. (2008). Identifying therupture plane of the 2001 Nisqually Washington, earthquake. Bulletin of the SeismologicalSociety of America, 98, 1546–1558., H., Shan, S. J., Dragert, H. and G. Rogers, (2009). Northern Cascadia episodic tremorand slip: A decade of tremor observations from 1997 to 2007, J. Geophys. Res., 114, B11,B00A12.Katayama, I., Terada, T., Okazaki, K., and Tanikawa, W. (2012). Episodic tremor and slowslip potentially linked to permeability contrasts at the Moho. Nature Geoscience, 5, 731–, S., Katayama, I., and Okazaki, K. (2011). Permeability anisotropy of serpentine andfluid pathways in a subduction zone, Geology, 39, 939-942.Kirby, S., Wang, K., and Dunlop, S. (2002). The Cascadia Subduction Zone and Related Sub-duction Systems-Seismic Structure, Intraslab Earthquakes and Processes, and EarthquakeHazards. Usgs., S., Engdahl, E.R., and Denlinger, R. (1996). Intermediate-depth intraslab earthquakesand arc volcanism as physical expressions of crustal and uppermost mantle metamorphismin subducting slabs, in Subduction: Top to Bottom, Geophys. Monogr. Ser., vol. 96, editedby G. E. Bebout et al., pp. 195-214, AGU, Washington, D.C., doi:10.1029/GM096p0195Kissling, E. (1988). Geotomography with local earthquake data.Kissling, E., Ellsworth, W. L., Eberhart-Phillips, D., and Kradolfer, U. (1994). Initial refer-ence models in local earthquake tomography. Journal of Geophysical Research: Solid Earth,99(B10), 19635–19646., F. W., (2002). User’s Guide to HYPOINVERSE-2000, a Fortran Program to Solve forEarthquake Locations and Magnitudes, U.S. Geol. Surv. Open File Rep. 02-171.Kreemer, C., Govers, R., Furlong, K. P., and Holt, W. E. (1998). Plate boundary deformationbetween the Pacific and North America in the Explorer region. Tectonophysics, 293(3–4),225–238. Rocca, M., McCausland, W., Galluzzo, D., Malone, S., Saccorotti, G., and Del Pezzo,E. (2005). Array measurements of deep tremor signals in the Cascadia subduction zone.textitGeophysical Research Letters, 32, 1–4., T.J., W.H. Bentkowski, E.E. Davis, R.D. Hyndman, J.G. Souther, and J.A. Wright.(1988). Subduction of the Juan de Fuca Plate: thermal consequences, J. Geophys. Res., 93,15,207-15,225.106BIBLIOGRAPHYLewis, T. J., Lowe, C., and Hamilton, T. S. (1997). Continental signature of a ridge-trench-triple junction: Northern Vancouver Island. Journal of Geophysical Research, 102(B4), 7767., C. E. (2004). The chemistry of subduction-zone fluids. textitEarth and PlanetaryScience Letters, 223(1–2), 1–16., P. A., Blair, J. L., Waldhauser, F., and Oppenheimer, D. H. (2012). Juan de Fucaslab geometry and its relation to Wadati-Benioff zone seismicity. textitJournal of GeophysicalResearch, 117, B09306., P. A., Constantz, J. E., Hunt, A. G., and Blair, J. L. (2016). Helium as a tracer for flu-ids released from Juan de Fuca lithosphere beneath the Cascadia forearc. textitGeochemistryGeophysics Geosystems, 17, 2825–2834., P. A., Hyndman, R. D., and Blair, J. L. (2014). Relationship between the Cascadiaforearc mantle wedge, nonvolcanic tremor, and the downdip limit of seismogenic rupture.textitGeochemistry Geophysics Geosystems, 15, 1071-1095,, C. E. (2004). The chemistry of subduction-zone fluids. Earth and Planetary ScienceLetters, 223, 1-16.Monger, J.W.H and Brown, E.H. (2016). Tectonic evolution of the Southern Coast-CascadeOrogen, Northwestern Washington and Southwestern British Columbia, in Ed. E.S. Cheney,The Geology of Washington and Beyond: From Laurentia to Cascadia, University of Wash-ington Press, Seattle U.S.A.Monger, J.W.H., and Journeay, J.M., 1994, Basement geology and tectonic evolution of theVancouver region, in Monger, J.W.H., ed., Geology and geological hazards of the Vancouverregion, southwestern British Columbia: Geological Survey of Canada Bulletin 481, p. 3–25.Morishige, M. and van Keken, P. E. (2017). Along-arc variation in short-term slow slip eventscaused by 3-D fluid migration in subduction zones. Journal of Geophysical Research, 122,1434–1448., V., Kissling, E., Husen, S., and Quintero, R. (2010). Detection of systematic errors intravel-time data using a minimum 1D model: application to Costa Rica seismic tomography.Bulletin of the Seismological Society of America, 100, 629–639.Morell, K. D., Regalla, C., Leonard, L. J., Amos, C., and Levson, V. (2017). QuaternaryRupture of a Crustal Fault beneath Victoria, British Columbia, Canada. GSA Today, 27,4–10., G., Bostock, M. G., Christensen, N. I. and J. Tromp, (2014). Crustal anisotropyin a subduction zone forearc : Northern Cascadia, J. Geophys. Res. Solid Earth, 119, 7,2169-9356.McMechan, G. A., and Spence, G. D. (1983). P-wave velocity structure of the Earth’scrust beneath Vancouver Island. Canadian Journal of Earth Sciences, 20(5), 742–752., A. and E. E. Brodsky, (2008). Deep low-frequency tremor that correlates with pass-ing surface waves, J. Geophys. Res., 113(B1):B01307.107BIBLIOGRAPHYMuller, J E, Northcote, K E and Carlisle, D. (1973). Geological Survey of Canada, Open File170, 1973, 110 pages (1 sheet),, J. W., and Tiffin, D. L. (1974). Patterns of deformation sedimentation and tectonism,southwestern Canadian continental margin. Annales de la Societe Geologique de Belgique,97, pp. 169-183.Nakajima, J., and Hasegawa, A. (2016). Tremor activity inhibited by well-drained conditionsabove a megathrust. Nature Communications, 7, 13863.ć, M. R., Hyndman, R. D., Kumar, R., and Spence, G. D. (2003). Reflection signatureof seismic and aseismic slip on the northern cascadia subduction interface. Nature, 424, 416-420.Nedimović, M. R., Bohnenstiehl, D. R., Carbotte, S. M., Pablo Canales, J., and Dziak, R.P. (2009). Faulting and hydration of the Juan de Fuca plate system. Earth and PlanetaryScience Letters, 284, 94–102., T., Bostock, M. G. and J. F. Cassidy, (2005). New constraints on subduction zonestructure in northern Cascadia, Geophys. J. Int., 161, 3, 849-859.Nixon, G.T., Hammack, J.L., Hamilton, J.V., Jennings, H., Larocque, J.P., Orr, A.J., Friedman,R.M., Archibald, D.A., Creaser, R.A., Orchard, M.J., Haggart, J.W., Tipper, H.W., Tozer,E.T., Cordey, F., and McRoberts, C.A., 2011. Geology, geochronology, lithogeochemistryand metamorphism of the Mahatta Creek area, northern Vancouver Island (NTS 92L/05)1:50 000 scale. British Columbia Ministry of Energy and Mines, British Columbia GeologicalSurvey Geoscience Map 2011-03, 1:50,000 scale.Nowack, R. L., and M. G. Bostock (2013). Scattered waves from low-frequency earthquakes andplate boundary structure in northern Cascadia, Geophys. Res. Lett., 40, 16, 4238-4243.Obana, K., Scherwath, M., Yamamoto, Y., Kodaira, S., Wang, K., Spence, G. D., and Kao,H. (2015). Earthquake activity in northern Cascadia subduction zone off vancouver islandrevealed by ocean-bottom seismograph observations. Bulletin of the Seismological Society ofAmerica, 105(1), 489–495., K., (2002). Nonvolcanic deep tremor associated with subduction in southwest Japan,Science, 296, 1679-1681.Obara, K., Tanaka, S., Maeda, T., and T. Matsuzawa (2010). Depth-dependent activity ofnon-volcanic tremor in southwest Japan, Geophys. Res. Lett., 37, L13306.Okazaki, K., and Katayama, I. (2015). Slow stick slip of antigorite serpentinite under hydrother-mal conditions as a possible mechanism for slow earthquakes. Geophysical Research Letters,42, 1099–1104., C. C., and M. A. Saunders (1982). LSQR: Sparse linear equations and least squaresproblems, ACM Transactions on Mathematical Soft- ware 8/2, 195–209.Patrick Paitz, Alexey Gokhberg and Andreas Fichtner. (2018). A neural network fornoise correlation classification, Geophysical Journal International, 212, 2, p. 1468–1474,, S.M. (2001). Are the lower planes of double seismic zones caused by serpentine dehy-dration in subducting oceanic mantle? Geology, 29, 299-302.Peacock, S. M., Christensen, N. I., Bostock, M. G., and Audet, P. (2011). High pore pres-sures and porosity at 35 km depth in the Cascadia subduction zone. Geology, 39, 471–, S.M., and Wang, K. (1999). Seismic consequences of warm versus cool subductionmetamorphism: Examples from Southwest and Northeast Japan. Science, 286, 937–939Phillips, G. N., and Powell, R. (2010). Formation of gold deposits: a metamorphic devolatiliza-tion model. Journal of Metamorphic Geology, 28, 689-718.Plourde A. P., Bostock, M. G., Audet, P., and Thomas, A. M. (2015). Low-frequency earth-quakes at the southern Cascadia margin. Geophysical Research Letters, 42, 4849-4855.Pommier, A., and Evans, R. L. (2017). Constraints on fluids in subduction zones from electro-magnetic data. Geosphere, 13, 1026–1049., L. A., Creager, K. C., Crosson, R.S., Brocher, T. M., and Trehu, A. M. (2003).Intraslab earthquakes: dehydration of the Cascadia slab. Science, 302, 1197-1200.Peng, Z. and K. Chao, (2008). Non-volcanic tremor beneath the Central Range in Taiwantriggered by the 2001 Mw 7.8 Kunlun earthquake, Geophys. J. Int., 175, 2, 825-829.Ramachandran, K., Dosso, S.E., Spence, G.D., Hyndman, R.D., and Brocher,T.M. (2005). Forearc structure beneath southwestern British Columbia: A three-dimensional tomographic velocity model. Journal of Geophysical Research, 110,, K., and Hyndman, R. D. (2012). The fate of fluids released from subductingslab in northern Cascadia. Solid Earth, 3, 121–129., N., Fichtner, A., Sambridge, M., and Young, M. K. (2014). Seismic To-mography and the Assessment of Uncertainty. Advances in Geophysics, 55, 1–76. Reynen and Pascal Audet. (2017). Supervised machine learning on a network scale:application to seismic event classification and detection, Geophysical Journal International,210, 3, p. 1394–1409,, R. P. (1977). A model for recent plate interactions off Canada’s west coast. Cana-dian Journal of Earth Sciences, 14(3), 384–396., R. P., Currie, R. G., and Hyndman, R. D. (1980). The Dellwood Knolls and theirrole in triple junction tectonics off northern Vancouver Island. Canadian Journal of EarthSciences, 17, pp. 577-593.Robert, F., Boullier, A.-M., Firdaous, K. (1995). Gold–quartz veins in metamorphic terranesand their bearing on the role of fluids in faulting. Journal of Geophysical Research, 100,12,861–12,879.Roecker, S. W., Seismicity and Tectonics of the Pamir-Hindu Kush Region of Central Asia,Ph.D. Thesis, Mas- sachusetts Institute of Technology, Cambridge, MA, 294 p., 1981.109BIBLIOGRAPHYRoecker, S. W., Tomography in zones of collision: practical considerations and examples, in:Seismic Tomography: Theory and Practice, eds. H. M. Iyer and K. Hirahara, Chapman andHall, London, 584−612, 1993.Rogers, G., and Dragert, H. (2003). Episodic tremor and slip on the Cascadia subduction zone:the chatter of silent slip. Science, 300, 1942–1943., G. C., and Hasegawa, H. S. (1978). A second look at the British Columbia earthquakeof June 23, 1946. Bulletin of the Seismological Society of America, 68(3), 653–676. Retrievedfrom, K. M. M., and Furlong, K. P. (1995). Ephemeral plate tectonics at the QueenCharlotte triple junction. Geology, 23(11), 1035–1038.<1035:EPTATQ>2.3.CO;2Rondenay, S., Abers, G. A., and van Keken, P. E. (2010). Seismic imaging of subduction zonemetamorphism. Geology, 36, 275-278.Royer, A. A. and M. G. Bostock, (2013). A comparative study of low frequency earthquaketemplates in northern Cascadia, Earth Planet. Sci. Lett., 402, 247-256.Royer, A. A., Thomas, A. M., and Bostock, M. G. (2015). Tidal modulation and triggeringof low- frequency earthquakes in northern Cascadia. Journal of Geophysical Research, 120,384–405., A. M. and J. Armbruster, (2013). Imaging Slow Slip Fronts in Cascadia With HighPrecision Cross-Station Tremor Locations, G-cubed, 14, 12, 5371-5392.Rubinstein, J. L., La Rocca, M., Vidale, J. E., Creager, K. C., and Wech, A. G. (2008). Tidalmodulation of nonvolcanic tremor. Science, 319, 186–189., J. L., Vidale, J. E., Gomberg, J., Bodin, P., Creager, K. C. and S. D. Malone,(2007). Non-volcanic tremor driven by large transient shear stresses, Nature, 448, 579-82.Shelly, D. R. Beroza, G. C. and S. Ide, (2007). Non-volcanic tremor and low frequency earth-quake swarms, Nature, 446, 305-7.Shelly, D. R., Beroza, G. C., Ide, S. and Nakumula S. (2006). Low-frequency earthquakes inShikoku, Japan, and their relationship to episodic tremor and slip. Nature, 442,, R. H. (2009). Rupturing in overpressured during compressional inversion – the casefrom NE Honshu, Japan, Tectonophysics, 473, 404-416.Sibson, R. H., Robert, F., and Poulsen, H. K. (1988). High-angle reverse faults, fluid-pressurecycling and mesothermal gold-quartz deposits, Geology, 16, 551-555.Sibson, R. H. and Scott, J. (1998). Stress/fault controls on the containment and release ofoverpressured fluids: Examples from gold-quartz vein systems in Juneau, Alaska; Victoria,Australia and Otago, New Zealand, Ore Geology Reviews, 13, 293-306.Spence, G. D., Clowes, R. M., and Ellis, R. M. (1985). Seismic structure across the activesubduction zone of western canada. Journal of Geophysical Research: Solid Earth, 90(B8),6754–6772., C., G. C. Rogers, and J. F. Cassidy (1997). The 1978 Brooks peninsula, VancouverIsland, earthquakes, Bull. Seismol. Soc. Am., 87, 1011–1023.Theunissen, T., Chevrot, S., Sylvander, M., Monteiller, V., Calvet, M., Villaseñor, A., andGrimaud, F. (2018). Absolute earthquake locations using 3-D versus 1-D velocity modelsbelow a local seismic network: example from the Pyrenees. Geophysical Journal International,212(3), 1806-1828. M. Tréhu, Richard J. Blakely, Mark C. Williams; Subducted seamounts and re-cent earthquakes beneath the central Cascadia forearc. Geology ; 40 (2): 103–106. doi:, A. (2012). Understanding and parameter setting of STA-LTA trigger algorithm. NewManual of Seismological Observatory Practice 2 (NMSOP-2),1-80.Trugman, D. T., and Shearer, P. M. (2017). GrowClust: A hierarchical clustering algo-rithm for relative earthquake relocation, with application to the Spanish Springs andSheldon, Nevada, earthquake sequences, Seismological Research Letters, 88,, A., and Linde, N. (2006). Local earthquake (LE) tomography with joint inversionfor P- and S-wave velocities using structural constraints, Geophysical Research Letters, 33,L07303, doi:10.1029/2005GL025485.Thomas, A. M., Bürgmann, R., Shelly, D. R., Beeler, N. M. and M. L. Rudolph, (2012). Tidaltriggering of low frequency earthquakes near Parkfield, California: Implications for faultmechanics within the brittle-ductile transition, J. Geophys. Res., 117(B05301).Thurber, C.H., and Eberhart-Phillips, D. (1999), Local earthquake tomography with flexiblegridding, Computers and Geosciences 25, 809–818.Um, J. and Thurber, C.H. (1987), A fast algorithm for two-point seismic ray tracing, Bull.Seism. Soc. Am. 77, 972–986.U.S. Geological Survey, Earthquake Hazards Program, 2017, Preliminary Determination ofEpicenters (PDE) Bulletin: U.S. Geological Survey,, J. C. and R. S. Crosson, (1990). Determination of teleseismic relative phase arrivaltimes using the multi-channel cross-correlation and least squares, Bull. Seismol. Soc. Am.,80, 1, 150-169.Vidale, J. E., and Shearer, P. M. (2006). A survey of 71 earthquake bursts across southernCalifornia: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers,Journal of Geophysical Research, 111, B0531, doi:10.1029/2005JB004034Waldhauser, F. (2001). hypoDD – A Program to Compute Double-Difference Hypocenter Lo-cations. U.S. Geol. Survey Open File Report 01-113, pp 25.Waldhauser, F. and Ellsworth, W. L (2000). A double-difference earthquake location algorithm:Method and application to the northern Hayward fault. Bulletin of the Seismological Societyof America, 90, 1353-1368.111BIBLIOGRAPHYWang, K. (2000) Stress–strain ‘paradox’, plate coupling, and forearc seismicity at the Cascadiaand Nankai subduction zones. Tectonophysics, 319, 321-338.Wang, K., R. Wells, S. Mazzotti, R. D. Hyndman, and T. Sagiya (2003), A revised disloca-tion model of interseismic deformation of the Casca- dia subduction zone, J. Geophys. Res.,108(B1), 2026, doi:10.1029/2001JB001227.Wang, K., Mulder, T., Rogers, G. C., and Hyndman, R. D. (1995). Case for very low couplingstress on the Cascadia Ssubduction Fault. Journal of Geophysical Research: Solid Earth,100(B7), 12907–12918., A. G. and K. C. Creager, (2008). Automated detection and location of Cascadia tremor,Geophys. Res. Lett., 35, 20, L20302.Wech, A. G., (2010). Interactive tremor monitoring, Seismol. Res. Lett., 81, 4, 664-669.Wech, A. G. (2016). Extending Alaska’s plate boundary: Tectonic tremor generated by Yakutatsubduction, Geology, 44, 597-590.Wech, A. G. and K. C. Creager, (2011). A continuum of stress, strength and slip in the Cascadiasubduction zone, Nature Geoscience, 4, 9, 624-628.Wells, R. E., Blakely, R. J., Wech, A. G., McCrory, P. A. and Michael, A. (2017). Cascadiasubduction tremor muted by crustal faults. Geology, 45, 515–518., D. S. (1993), Confidence intervals for motion and deformation of the Juan de FucaPlate, J. Geophys. Res., 98(B9), 16053–16071, doi: 10.1029/93JB01227.Wilson, D. S. (2002), The Juan de Fuca plate and slab: Isochron structure and Cenozoic platemotions, in The Cascadia Subduction Zone and Related Subduction Systems, edited by S.Kirby, K. Wang, and S. Dunlop, U.S. Geol. Surv. Open File Rep., 02-328, 9–12.Zelt, B. C., Ellis, R. M., Clowes, R. M., and Hole, J. A. (1996). Inversion of three-dimensionalwide-angle seismic data. Journal of Geophysical Research, 101, 8503–8529.Zhang, H., and Thurber, C. H. (2003). Double-difference tomography: The method and itsapplication to the Hayward Fault, California. Bulletin of the Seismological Society of America,93, 1875–1889., H., and Thurber, C. H. (2007). Estimating the model resolution matrix for large seis-mic tomography problems based on Lanczos bidiagonalization with partial reorthogonal-ization. Geophysical Journal International, 170(1), 337–345. ASupplementary material for Chapter 2A.1 Effect of Station GeometryTo illustrate the effect of station geometry on LFE detection sensitivity between different ETSyears, we include the station density and maximum delay time maps for 2003, 2005 and 2013below. Station geometry for 2003 shows decreased coverage near the west coast of VancouverIsland, consistent with epicenter density observations.Figure A.1: (a) Map of the number of stations (inverted triangles) with pre-critical arrivalsfor hypocenters located on the plate interface model of Audet et al. (2010) undersouthern Vancouver Island for the stations available during the 2003 ETS episode.The dashed lines indicate the 20, 30 and 40 km depth contours of the plate interface.LFE detections are more likely in regions of high station density statisfying the pre-critical arrival criterion. (b) Same as (a) but for the maximum delay between S-wavearrival times among all pairs of stations.113A.2. Epicenter Locations for Northern WashingtonFigure A.2: Same as figure A.1, but for the 2005 ETS episode.Figure A.3: Same as figure A.1, but for the 2013 ETS episode.A.2 Epicenter Locations for Northern WashingtonWe analysed data from the Washington state CAFE and Array of Arrays (Ghosh et al., 2012)experiments from ETS episodes in 2010 and 2011. Episodes at these stations lasted approx-imately 18 and 20 days, and resulted in 2176 and 1568 detections. LFEs from Washingtonstate exhibit lower SNR and waveform similarity than those in southern Vancouver Island,with higher noise levels between 1-2 Hz, in particular. To mitigate this shortcoming, we have114A.3. Time-distance Plots for the 2003 and 2005 Episodesexperimented with different bandpass filters, but such filters often reduce the impulsiveness ofthe waveforms for a lower corner frequency higher than about 1.5 Hz. In addition, LFE wave-forms in Washington state appear to have longer coda suggesting a stronger structural imprintthat reduces waveform similarity at nearby stations. This factor decreases the performance ofour approach, as correlation coefficients are smaller and timing uncertainties are larger. LFEepicenters for both the 2010 and 2011 episodes are shown in figure S4 and compared with thetremor locations of Wech and Creager (2008). The main propagation front is clearly evidentfor both episodes though we obtain far fewer detections into the Strait of Juan de Fuca.For northern Washington, the calculated density of epicenters for the ETS episodes of 2010and 2011 is displayed in figure S5. The two asperities observed correspond geographicallyto tremor patches in the catalogue of Wech (2010). Ghosh et al. (2009) have also reportedasperities in the same areas for the 2008 ETS episode.It must be noted that because of station separation constraints and computational efficiency,we separated the stations in southern Vancouver Island and northern Washington in 2 differentsets for which we run the algorithm separately, such that we cannot monitor LFEs continouslyacross a broad region including the Juan de Fuca strait as in the Wech tremor catalogue.A.3 Time-distance Plots for the 2003 and 2005 EpisodesFor completeness, we include below the time-distance plots with comparison to the templatedetections as in figure 8 for the 2003 and 2005 ETS episodes.115A.3. Time-distance Plots for the 2003 and 2005 Episodesa) 124oW  40’  20’  123oW  40’  20’   47oN  20’  40’   48oN  20’ 0 20 40 kmb) 124oW  40’  20’  123oW  40’  20’ 0 20 40 km  08−Aug 10−Aug 12−Aug 14−Aug 15−Aug 17−Aug 19−Aug 21−Aug 23−Aug 24−Augc) 124oW  40’  20’  123oW  40’  20’   47oN  20’  40’   48oN  20’ 0 20 40 kmc) 124oW  40’  20’  123oW  40’  20’ 0 20 40 km  05−Aug 07−Aug 09−Aug 11−Aug 13−Aug 15−Aug 17−Aug 19−Aug 21−Aug 23−AugFigure A.4: LFE epicenters for the 2010 (a) and 2011 (b) ETS episodes in northern Wash-ington. Left panels are from the present study, and right panels show the tremordetections of Wech and Creager (2008).116A.3. Time-distance Plots for the 2003 and 2005 Episodes 124oW  40’  20’  123oW  40’  20’   47oN  20’  40’   48oN  20’ 0 20 40 km  ETS 201001234567 124oW  40’  20’  123oW  40’  20’   47oN  20’  40’   48oN  20’ 0 20 40 km  ETS 2011012345630 km40 kmFigure A.5: Epicenter density maps of northern Washington state computed using thekernel density estimator of Botev et al. (2010) for the 2010 (a) and 2011 (b) ETSepisodes. Inverted triangles are the stations used to locate the events for each re-spective period.117A.3.Time-distancePlotsforthe2003and2005Episodes0204060a)Event count−50050b)Distance along strike [km]Feb−26 Feb−28 Mar−02 Mar−04 Mar−06 Mar−08 Mar−10−50050c)Distance along strike [km]  Distance along dip [km]−50 −40 −30 −20 −10 0  10 20 30 40 Figure A.6: Cross-station detections compared with template detections for the 2003 ETS episode. (a) Number of cross-stationdetections with respect to time. (b) Distance along-strike is plotted against time for our cross-station detections, and datapoints are color-coded for the distance along dip. The profiles used to compute along-strike and along-dip distances are shownin figure 10 (Dark colors are updip, clear colors are downdip). (c) Same as (b) but for template detections. RTRs identifiedby Royer and co-workers (Royer and Bostock, 2013; Royer et al., 2015) are shown in black outlines.118A.3.Time-distancePlotsforthe2003and2005Episodes07−Sep−2005 09−Sep−2005 11−Sep−2005 13−Sep−2005 15−Sep−2005 17−Sep−2005 19−Sep−2005 21−Sep−2005 23−Sep−2005050100150a)Event count−50050b)Distance along strike [km]Sep−07 Sep−09 Sep−11 Sep−13 Sep−15 Sep−17 Sep−19 Sep−21 Sep−23−50050c)Distance along strike [km]  Distance along dip [km]−50 −40 −30 −20 −10 0  10 20 30 40 Figure A.7: Same as figure A.6, but for the 2005 ETS episode.119Appendix BSupplementary material for Chapter 3The following figures, B.1 to B.1 complement Chapter 3.Figure B.1 provides the S-wave velocity model corresponding to sections AA’, BB’ and CC’.We also present two different data resolution tests.First, we performed a checkerboard resolution test with a +/- 5% anomaly in P-wave velocityon 2x2x2 grid cells, which correspond to cells of 24 x 24 x 6 km3. Figures B.2 and B.3 displayour results for sections AA’, BB’ and CC’.Second, we constructed a synthetic model comprising a background 1D Vp model with asuperimposed dipping layer 8 km thick with a Vp/Vs ratio of 2.35 approximately following theslab interface models of Audet et al. (2010) and McCrory et al. (2012). Figures B.4 and B.5display our results for sections AA’, BB’ and CC’.Figure B.6 provides average forearc crust values for Poisson’s ratio and Vp for profiles AA’,BB’ and CC’. The average Poisson’s ratio is computed from the ts/tp traveltimes for a vertical1D column sampled from our velocity models between 10 km depth and the top of the slab ora nominal continental Moho of 36 km landward the mantle wedge corner.120Figure B.1: Vs model for profiles AA’, BB’ and CC’. Symbols have the same significationas in figures 3.2 and 3.3 in the manuscript.121Figure B.2: Synthetic checkerboard model. Note that profiles cross checkerboard gridobliquely leading to some fading for certain sections.122Figure B.3: Recovered checkerboard model. Relocated seismicity in the checkerboardmodel is only shown to illustrate data coverage.123Figure B.4: Synthetic high Vp/Vs LVZ model. We interpolated an 8km dipping layer witha Vp/Vs ratio of 2.35 approximately following the slab interface models of Audet etal. (2010) and McCrory et al. (2012) onto our model grid.124Figure B.5: Recovered high Vp/Vs dipping layer. Relocated seismicity in the syntheticmodel is shown to illustrate data coverage.125Figure B.6: Average Poisson’s ratio and P-wave velocity along profiles AA’, BB’ and CC’.The Poisson’s ratio is computed from the ts/tp traveltimes for a vertical 1D columnsampled from our velocity models between 10 km depth and the top of the slab or anominal continental Moho of 36 km landward the mantle wedge corner. The verticallines delineate the approximate lateral extent of prominent Poisson’s ratio anomalyin figure 3.2.126


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items