Seismic and mechanical attributes of lithospheric deformation and subduction in western Canada by Pascal Audet B.Sc., Université de Montréal, 2002 M.Sc., Université du Québec à Montréal, 2004 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Geophysics) The University Of British Columbia (Vancouver) October, 2008 c© Pascal Audet 2008 Abstract Convergent continental margins are regions of intense deformation caused by the interaction of oceanic plates with continents. The spatial extent of deformation is broadly commensurate with the specific time scale of the causative phenomenon. For example, subduction-related short-term defor- mation is limited to <200 km from the margin, whereas long-term plate con- vergence cause deformation over ∼1000 km landward. Deformation is thus manifested in multiple ways, with attributes depending on the scale of mea- surement. In this thesis we investigate the use of two geophysical approaches in the study of deformation: 1) The analysis of potential-field anomalies to derive estimates of the elastic thickness (Te) of the lithosphere, and 2) The structural study of past and present subduction systems using seismic observations and modelling. Both approaches involve the development of appropriate methodologies for data analysis and modelling, and their ap- plication to the western Canadian landmass. Our findings are summarized as follows: 1) We develop a wavelet-based technique to map variations in Te and its anisotropy; 2) We show how a step-wise transition in Te and its anisotropy from the Cordillera to the Craton is a major factor influencing lithospheric deformation; 3) We implement a waveform modelling tool that includes the effects of structural heterogeneity and anisotropy for teleseismic applications, and use it to model the signature of a fossil subduction zone in a Paleoproterozoic terrane; 4) We use teleseismic recordings to map slab edge morphology in northern Cascadia and show how slab window tectonism and slab stretching led to the creation of the oceanic Explorer plate; 5) We use seismic signals from the subducting oceanic crust to calculate elevated Poisson’s ratio and infer high pore-fluid pressures and a low-permeability plate boundary within the forearc region of northern Cascadia. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xvi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Theme and objectives . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Wavelet analysis of Bouguer gravity and topography . . . 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Wavelet analysis . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Wavelet transform . . . . . . . . . . . . . . . . . . . . 14 2.2.2 Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.3 Cross spectral analysis . . . . . . . . . . . . . . . . . 18 2.2.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 19 iii Table of Contents 2.2.5 Limitations . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.6 Isotropic coherence inversion . . . . . . . . . . . . . . 21 2.2.7 Anisotropic coherence inversion . . . . . . . . . . . . 23 2.3 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.1 Uniform rigidity pattern . . . . . . . . . . . . . . . . 25 2.3.2 Spatially varying rigidity pattern . . . . . . . . . . . 28 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 31 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3 Deformation of continents at convergent margins . . . . . 36 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Te estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.1 Coherence method . . . . . . . . . . . . . . . . . . . . 38 3.2.2 Multitaper analysis . . . . . . . . . . . . . . . . . . . 40 3.2.3 Wavelet transform . . . . . . . . . . . . . . . . . . . . 42 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.1 Isotropic Te . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.2 Anisotropic Te . . . . . . . . . . . . . . . . . . . . . . 49 3.3.3 Geophysical correlations . . . . . . . . . . . . . . . . 52 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4.1 Te, Tc, and heat flow . . . . . . . . . . . . . . . . . . 54 3.4.2 Anisotropy and stress: A primer . . . . . . . . . . . . 57 3.4.3 Conceptual model of deformation . . . . . . . . . . . 59 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.6 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 63 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4 Teleseismic waveform modelling . . . . . . . . . . . . . . . . 70 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2.1 Acoustic one-way operators . . . . . . . . . . . . . . . 72 iv Table of Contents 4.2.2 Elastic one-way operators . . . . . . . . . . . . . . . . 75 4.2.3 Practical considerations . . . . . . . . . . . . . . . . . 77 4.3 Synthetic examples . . . . . . . . . . . . . . . . . . . . . . . 81 4.3.1 1D crustal structure . . . . . . . . . . . . . . . . . . . 81 4.3.2 2D subduction model . . . . . . . . . . . . . . . . . . 83 4.3.3 Remarks on computational efficiency . . . . . . . . . 89 4.4 Application to data from NW Canada . . . . . . . . . . . . . 89 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.6 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 95 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5 Morphology of slab edge in northern Cascadia . . . . . . . 100 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.2 Shallow structure from receiver functions . . . . . . . . . . . 102 5.3 Deep structure from P-wave tomography . . . . . . . . . . . 105 5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.4.1 Further evidence of a sharp slab edge . . . . . . . . . 107 5.4.2 Explorer tectonics: Slab stretching? . . . . . . . . . . 108 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.6 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 111 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6 Overpressured subducted oceanic crust . . . . . . . . . . . . 115 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.5 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 123 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 v Table of Contents 7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 7.1 Main contributions . . . . . . . . . . . . . . . . . . . . . . . 127 7.1.1 Wavelet analysis of potential field data . . . . . . . . 127 7.1.2 Large-scale deformation of convergent margins . . . . 128 7.1.3 Teleseismic waveform modelling . . . . . . . . . . . . 129 7.1.4 Slab morphology in northern Cascadia . . . . . . . . 129 7.1.5 Subducted water and megathrust permeability . . . . 130 7.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.2.1 Mechanical models of the lithosphere . . . . . . . . . 130 7.2.2 Teleseismic attributes of subducted slabs . . . . . . . 132 7.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.3.1 Spherical wavelet transform . . . . . . . . . . . . . . 133 7.3.2 Deformation at triple junction margins . . . . . . . . 133 7.3.3 Gorda slab morphology . . . . . . . . . . . . . . . . . 134 7.3.4 Cascadia forearc structure and ETS . . . . . . . . . . 135 7.3.5 Water transport and seismogenesis . . . . . . . . . . 136 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Master Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . 141 Appendices A Approximate amplitude normalizing factor . . . . . . . . . 159 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 B Supplementary Information for Chapter 5 . . . . . . . . . . 165 B.1 Receiver functions . . . . . . . . . . . . . . . . . . . . . . . . 165 B.2 Phase stacking . . . . . . . . . . . . . . . . . . . . . . . . . . 168 B.3 Tomography cross-sections . . . . . . . . . . . . . . . . . . . 170 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 vi Table of Contents C Supplementary Information for Chapter 6 . . . . . . . . . . 174 C.1 Receiver function analysis . . . . . . . . . . . . . . . . . . . . 174 C.2 Phase weighted stacking . . . . . . . . . . . . . . . . . . . . . 174 C.3 Auto-correlation stacking . . . . . . . . . . . . . . . . . . . . 176 C.4 Plane-wave travel-times . . . . . . . . . . . . . . . . . . . . . 178 C.5 Limitating cases . . . . . . . . . . . . . . . . . . . . . . . . . 182 C.6 Frequency bias . . . . . . . . . . . . . . . . . . . . . . . . . . 183 C.7 Bootstrapping . . . . . . . . . . . . . . . . . . . . . . . . . . 183 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 vii List of Tables 4.1 Properties of subduction zone elements . . . . . . . . . . . . . 84 viii List of Figures 2.1 Results from uniform and isotropic synthetic data . . . . . . 27 2.2 Results using uniform Te=20 km . . . . . . . . . . . . . . . . 28 2.3 Results from uniform anisotropic synthetic data . . . . . . . . 29 2.4 Results from spatially varying isotropic Te . . . . . . . . . . . 30 3.1 Examples of Te from the multitaper method . . . . . . . . . . 44 3.2 Examples of Te from the wavelet transform metho . . . . . . 45 3.3 Te structure of western Canada . . . . . . . . . . . . . . . . . 46 3.4 Anisotropic Te using the multitaper method . . . . . . . . . . 48 3.5 Anisotropic Te using the wavelet transform method . . . . . . 50 3.6 Anisotropic Te structure of western Canada . . . . . . . . . . 51 3.7 Geophysical correlations . . . . . . . . . . . . . . . . . . . . . 53 3.8 Anisotropy correlations . . . . . . . . . . . . . . . . . . . . . 55 3.9 Statistics of anisotropy . . . . . . . . . . . . . . . . . . . . . . 56 3.10 Geodynamic interpretation . . . . . . . . . . . . . . . . . . . 61 4.1 1D modelling results . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Accuracy of energy-flux normalization term . . . . . . . . . . 84 4.3 2D subduction zone model . . . . . . . . . . . . . . . . . . . . 85 4.4 2D modelling results . . . . . . . . . . . . . . . . . . . . . . . 86 4.5 Wiggle plots of 2D results . . . . . . . . . . . . . . . . . . . . 88 4.6 Data and station distribution in NW Canada . . . . . . . . . 90 4.7 Modelled and observed receiver functions . . . . . . . . . . . 92 4.8 Subsurface model of NW Canada . . . . . . . . . . . . . . . . 93 5.1 Tectonic map offshore western Canada . . . . . . . . . . . . . 101 ix List of Figures 5.2 Receiver function results . . . . . . . . . . . . . . . . . . . . . 104 5.3 Deep slab morphology . . . . . . . . . . . . . . . . . . . . . . 106 5.4 Tectonic interpretation of slab morphology . . . . . . . . . . . 110 6.1 Geographical map of northern Cascadia subduction zone . . . 117 6.2 Receiver function results in northern Cascadia . . . . . . . . . 119 6.3 Interpretation of receiver function results . . . . . . . . . . . 122 B.1 Ray diagram for scattered phases . . . . . . . . . . . . . . . . 166 B.2 Receiver functions at station PGC . . . . . . . . . . . . . . . 167 B.3 Polarity reversals . . . . . . . . . . . . . . . . . . . . . . . . . 169 B.4 Phase stacking for station PGC . . . . . . . . . . . . . . . . . 170 B.5 Tomography profile along A1-A2 . . . . . . . . . . . . . . . . 171 B.6 Tomography profile along B1-B2 . . . . . . . . . . . . . . . . 172 C.1 Filtered receiver functions . . . . . . . . . . . . . . . . . . . . 175 C.2 Depth and VP /VS results . . . . . . . . . . . . . . . . . . . . 177 C.3 Window extraction and auto-correlation for Ps . . . . . . . . 179 C.4 Window extraction and auto-correlation for Pps . . . . . . . . 180 C.5 Auto-correlation stacking at PGC . . . . . . . . . . . . . . . . 181 C.6 Frequency bias . . . . . . . . . . . . . . . . . . . . . . . . . . 184 x Preface This thesis is an amalgamation of published and submitted manuscripts jointly written with many co-authors from different institutions. The format of each chapter follows the regulations of the Faculty of Graduate Studies. The content was unaltered to conform with journal copyrights. Most figures were prepared using the Generic Mapping Tool (GMT) software. xi Acknowledgments First and foremost, I would like to express my gratitude to my thesis advisor, Michael Bostock. Michael has been a terrific supervisor from the first day I arrived at UBC with a vague notion that I wanted to study seismology. His insight on relevant scientific problems and his thorough knowledge of the field made this thesis an inspiring scientific endeavour. His availability and support are standard-setting in academia, and a great example to follow in my future career. I’ve had the pleasure to work closely with Mark Jellinek, whose scientific curiosity and taste for important problems got him interested in my previous graduate work and allowed me to gain experience in student supervision. Mark is also the key person who nominated me for the prestigious Miller Fellowship at UC Berkeley, of which he was himself a recipient. I regard Mark as a mentor, but also as a friend outside academia. I would also like to thank Ron Clowes, who made sure UBC was on the top of my list of potential institutions for pursuing my PhD. I also thanks him, along with Liz Hearn, James Scoates, and Simon Peacock, for serving on my supervisory committee. I am grateful to Liz Hearn for showing early interest in my work and for sending me to a conference. Other faculty who influenced my scientific experience during my PhD are my co-authors John Cassidy at PGC, Nikolas Christensen and Simon Peacock at UBC, and Jean-Claude Mareschal at UQAM. Jean-Claude hired me as an summer undergraduate intern and introduced me to the wonderful world of scientific research. The work we did together inspired me to further my education in geophysics. He subsequently supervised my MSc thesis and I am grateful for his constant encouragments and support regarding my ideas and abilities in research. xii Acknowledgments Having been involved in field deployment of seismic stations on Vancou- ver Island, I appreciate how difficult obtaining good quality data can be. For this reason I would like to express my gratitude to Issam Al-Koubbi who was chief technician for POLARIS in BC, to Jean-Philippe Mercier, Julien Chaput and Adam Baig for assistance in field work, and to all those other people who rarely get any recognition for doing all the hard work, providing us with our primary source of precious data. Life at UBC would not have been the same without the friendship of Jean-Philippe Mercier and Andrew Schaeffer who shared both the occasional frustrations and the frequent rewards of scientific research. They were not only great colleagues and officemates but also great friends to hang out with. Jean-Philippe made me appreciate the wonders of outdoor activities in BC and beyond, and Andrew was a great hockey mate. Other significant peo- ple who made my stay in Vancouver memorable are Rob Eso, Daniel Ross, Julien Chaput, and Etienne Poulin. I am proud of having rocked the EOS department as a member of Mr. Ross and Band. Of course several other people contributed to make the EOS department an enjoyable place, espe- cially the grad students and postdocs participating in Mark and Catherine’s weekly group meeting. In particular I would like to thank Hideharu Uno who was willing to suffer through complicated spectral analysis problems and my own inexperience as a supervisor. Outside of UBC I was fortunate to have an amazing group of friends who kept me (thankfully!) distracted from work. Many thanks to my parents, France and Denis, my sister Nancy, and my brother René. Not one word of this thesis would have been written without their never-ending love and support. I am especially honored to be the first member of the Audet and Desroches families up to my grand-parents to submit a PhD thesis. This is a tribute to the many sacrifices made by those people in their lives to make this acheivement possible. Finally, and most importantly, I owe a great deal to my wonderful wife Marie-Eve. We met during the course of this PhD and the many obstacles we had to cross brought us closer together, and recently got us married. I am well aware of the enormous sacrifices she had to make and I am forever xiii Acknowledgments grateful for her love and emotional support. These few words cannot do justice to describe the way she forever changed my life. Financial support was provided by the Natural Sciences and Engineering Research Council (NSERC, Canada), the Fonds Québécois de la Recherche sur la Nature et les Technologies (FQRNT, Québec), and the Cordula and Gunter Paetzold University Graduate Fellowship (UGF, UBC). xiv À Marie-Ève, France et Denis. xv Statement of Co-Authorship Chapter 2 was published with Jean-Claude Mareschal, my former M.Sc graduate advisor at the Université du Québec à Montréal. I identified the research program and performed the research while at UBC. I wrote the manuscript with inputs from Jean-Claude. Chapter 3 was published with Mark Jellinek and Hideharu Uno at UBC. I identified the research program and co-supervised (with Mark) Honours student Hideharu who performed the multitaper analysis. I performed the wavelet analysis and the comparative analysis. I wrote the manuscript in close collaboration with Mark. Chapter 4 was published with Michael Bostock and Jean-Philippe Mercier at UBC. Michael identified the research program. I performed the research, algorithm development, and synthetic modelling. Jean-Philippe contributed the analysis of real data. I wrote the manuscript in collaboration with Michael. Chapter 5 was published with Michael Bostock and Jean-Philippe Mercier at UBC, and John Cassidy at the Pacific Geoscience Centre of the GSC. I performed the major part of research and the interpretation of results, with inputs from the co-authors. Jean-Philippe provided the tomography results. I wrote the manuscript in collaboration with all co-authors. Chapter 6 was published with Michael Bostock, Nikolas Christensen, and Simon Peacock at UBC. Michael and I identified the research program. I performed the research in collaboration with Michael. Nikolas and Simon participated in the interpretation of results. I wrote the manuscript in close collaboration with Michael and with inputs from Nikolas and Simon. xvi Chapter 1 Introduction Plate tectonics (i.e. relative motion of rigid lithospheric plates) is the sur- face expression of deep mantle convection driven by the secular cooling of the Earth and radioactive heating of rocks [Turcotte and Schubert , 2002]. This recognition revolutionized geology and geophysics and indeed consti- tuted a paradigm shift in the Earth sciences [Kuhn, 1962]. In particular, interactions between adjacent plates across narrow boundaries (∼100 km) explain the occurence of earthquakes, distribution of volcanoes, and oroge- nesis. Among the three types of plate boundaries (convergent, divergent, transforming), convergent margins are regions of plate consumption by sub- duction that provide the most effective force driving plate motion [Forsyth and Uyeda, 1975]. Convergent continental settings, characterized by the subduction of oceanic plates underneath a continent, are of particular so- cietal interest due to potential earthquake and volcanic hazards affecting densely populated coastal regions, and the economic impacts of precious ore-body extraction produced in subduction zones. Repeated episodes of subduction have also contributed to the stabilization of continental litho- sphere [Jordan, 1978]. Characterizing structure and deformation that affect past and present active continental margins is thus central in understanding tectonic and geodynamic processes that shape the Earth surface. 1.1 Theme and objectives Convergent continental margins can be broadly separated into several zones of deformation from the coast to the continental interior. On a regional scale (>1000 km) we can define the (often allochthonous) Cordillera associated with fold and thrust belts, the adjacent foreland basin where thick sedimen- 1 Chapter 1. Introduction tary layers overlay an older (autochthonous) basement, and the Cratonic Shield which forms the backstop of deformation. These structures are formed over several geological eras marking alternating episodes of subduction with passive margin rifting [Thomas, 2006]. At a smaller scale (∼200 km), de- formation associated with subduction processes is confined to the arc and forearc, and sometimes backarc, depending on the type of convergent mar- gin (extensional or compressional) [Stern, 2002]. In a compressive setting, subduction-related short-term deformation occurs mostly in the forearc due to locking on the plate interface. The spatial configuration of deformation zones is thus commensurate with the time scale of the causative phenomenon. This character can be exploited to make inferences on the structure and mechanics of the litho- sphere at a convergent margin at different scales, which is the main theme of this thesis. Several geophysical approaches exist to address such questions. Ultimately, the study of plate tectonics relies on rheological models of the lithosphere, for which there exists a range of classes (elastic, viscous, plastic, brittle, mixed classes, etc.) that describe the response of the lithosphere to geological loads, depending on the time scale of deformation and loading [Ranalli , 1995]. The first approach employed in this thesis makes use of the fact that at time scales >1 Ma, the lithosphere behaves as an elastic plate in re- sponse to loading [Turcotte and Schubert , 2002; Watts, 2001]. Elasticity is the rheological class linearly relating stress and strain in solids. In other words, the amount of instantaneous deformation undergone by an elastic body is proportional to the constraint imposed on it. Such a description of the lithosphere implies that the topography of the Earth surface (i.e. ob- served deformation) is a direct reflection of forces acting at its boundaries and at depth. The main source of loading at time scales >1 Ma affect- ing stable plates is provided by the gravitational body force [Turcotte and Schubert , 2002]. An elastic lithospheric plate will respond to surface (topog- raphy) or internal (density contrasts) gravitational loading by flexure. The characteristic wavelength of flexure gives an indication of the plate rigidity, which, given reasonable values of Young’s (elastic) modulus, is related to an 2 Chapter 1. Introduction equivalent elastic thickness (Te). Mapping variations in Te can give valuable information on the tectonic evolution of continents. The technique generally used to estimate Te relies on the spectral rela- tionships between topography and gravity anomalies [Forsyth, 1985]. This approach has generated conflicting results on continents which has triggered a vigorous debate over the factors contributing to the strength of the litho- sphere [Burov and Diament , 1995; Burov and Watts, 2006; Jackson, 2002; Maggi et al., 2000]. Part of the problem resides in the fact that estimat- ing Te is a complex spatio-spectral problem that requires advanced spectral methods [Pérez-Gussinyé et al., 2004]. In the first part of this thesis the objectives are to develop a high-resolution spectral method based on the wavelet transform to improve the mapping of Te and apply it to data from western Canada to infer the mechanical controls of the lithosphere on con- tinental deformation at an active margin. The second approach utilized also makes use of the elasticity of rocks but at a much shorter time scale (seconds to hours) [Ranalli , 1995]. This behavior is inferred mainly from seismic waves generated by earthquakes and travelling through the Earth. By analyzing the speed at which seismic waves are propagating and the scattering pattern generated by their inter- action with structural heterogeneity, we can obtain information on elastic properties of the subsurface. This information can in turn be related to tectonic processes operating at inaccessible depths. One of the most popular techniques to image lithospheric structure is based on the scattering of teleseismic P-waves [Bostock , 2007]. Teleseismic waves are generated by earthquakes at epicentral distances >30◦ that travel through the mantle and crust and are recorded at the surface. When a P- wave interacts with a discontinuity in structure on the receiver side (i.e. an interface separating material with different seismic velocities), a secondary S wave is generated that arrives shortly after the direct P wave. Because the free-surface is a quasi-perfect reflector of seismic energy, the direct P wave is reflected back to the subsurface and further interacts with the disconti- nuity, generating reflected S waves that arrive later at the surface. This set of secondary pulses (i.e. receiver function) is recorded in the coda of the 3 Chapter 1. Introduction direct P wave and is recovered by deconvolution of the P wave train with radial and transverse S-component seismograms. Receiver functions are sen- sitive mostly to S velocity structure. In the second part of this thesis the objectives are to calculate and model receiver functions at several arrays of broadband seismic stations across western Canada in order to image litho- spheric structure of past and present subduction zones. The ultimate goal is to map morphology of the subducting plate and to use this information to resolve tectonic and geodynamic problems. 1.2 Outline In Chapter 2 we first give an overview of the problem of estimating Te with standard spectral methods due to the dependence of the calculation on the wavelength of flexure. We then propose and describe the use of a wavelet transform-based method to estimate Te and its anisotropy. We illustrate the performance of this method using synthetic data with varying patterns of Te and Te anisotropy. In Chapter 3 we are interested in the factors controlling lithospheric de- formation at a continental convergent margin from the Cordillera to the Cra- ton. We address this problem by using and comparing a standard window- based spectral technique (multitaper analysis) and the wavelet-based method on data from western Canada. Both methods also calculate the anisotropy in Te, a useful parameter that is largely ignored in Te studies. Anisotropy in Te measures a preferred direction of deformation and can be compared with different geophysical indicators of anisotropy. The increased resolution obtained with the wavelet method enables us to further compare our results with a set of geophysical observables from the region and formulate a con- ceptual model of deformation where spatial variations in Te have a leading order influence on the mechanics of continental convergent margins. In Chapter 4 we describe the implementation of a waveform modelling tool based on the one-way approximation of teleseismic wavefield propa- gation. The one-way wave equation is obtained from the full (two-way) wave equation by neglecting wavefield propagation in the reverse direction. 4 Chapter 1. Introduction Forward wave-mode coupling due to 1D structure and anisotropy is fully described by this approach, while the coupling of modes due to lateral het- erogeneity is approximated by a technique similar to the “phase-screen” method where the scattering coefficients are calculated locally and linearly interpolated between heterogeneous regions. We show how free-surface mul- tiples can be generated by making multiple passes through the model and reversing the sense of propagation. We demonstrate the application of the method using synthetic wavefields propagating through 1D and 2D mod- els. The technique is further employed to model waveforms from an array of broadband seismograph stations deployed in northwestern Canada as part of the CANOE seismic experiment. We interpret receiver function images that show a dipping layer of anisotropic material at a Paleo-Proterozoic suture zone to be the fossil signature of a subducted slab. Chapter 5 presents results from a recently deployed POLARIS array of broadband seismic stations on northern Vancouver Island, Canada. Offshore of Vancouver Island the oceanic Explorer micro-plate is thought to be an ephemeral adjustment to the accommodation of Pacific-Juan de Fuca-North America plate motions. Because plate kinematics result from a balance of forces acting at their boundaries, it is necessary to understand how structure of the descending slab affects plate reconfiguration in this region. We use receiver functions and P-wave travel-time tomography to map the morphol- ogy of the subducting slab in the forearc region beneath northern Vancouver Island and the interior of British Columbia. The results allow us to formu- late a tectonic model for the Explorer plate supported by geophysical and geological data that describes mechanisms that led to its creation and its ultimate fate. In Chapter 6 we are interested in the seismic signature of water trans- port within subducted oceanic crust of the Juan de Fuca slab in northern Cascadia. This is an important problem to investigate because water plays a key role in a variety of geodynamic processes, and especially in the origin and evolution of plate tectonics on Earth. We calculate and stack receiver functions over a POLARIS portable array of seismometers across southern Vancouver Island to estimate the Poisson’s ratio of subducted oceanic crust. 5 Chapter 1. Introduction This parameter is sensitive to the ratio between P and S velocities, and indirectly to pore-fluid pressure. Our results allow us to infer the presence of pore-fluids at elevated pressure and to estimate the permeability of the plate interface. This information is crucial in geodynamic modelling of fault slip and the mechanism of episodic tremor and slip. 6 Bibliography Bostock, M. G., Regional methods, in Seismology and structure of the Earth, Treatise on Geophysics, vol. 1, edited by G. Schubert, pp. 73–78, Elsevier, Amsterdam, The Netherlands, 2007. Burov, E. B., and M. Diament, The effective elastic thickness (Te) of con- tinental lithosphere: What does it really mean?, J. Geophys. Res., 100, 3905–3927, 1995. Burov, E. B., and A. B. Watts, The long-term strength of continental lithosphere: “jelly sandwich” or “crème brûlée”?, GSA Today, 16, 4–10, 2006. Forsyth, D., and S. Uyeda, On the relative importance of the driving forces of plate motion, Geophys. J. R. Astron. Soc., 43, 163–200, 1975. Forsyth, D. W., Subsurface loading and estimates of the flexural rigidity of continental lithosphere, J. Geophys. Res., 90, 12,623–12,632, 1985. Jackson, J., Strength of the continental lithosphere, GSA Today, 12, 4–10, 2002. Jordan, T. H., Composition and development of the continental tecto- sphere, Nature, 274, 544–548, 1978. Kuhn, T. S., The structure of scientific revolutions, University of Chicago Press, Chicago, USA, 1962. Maggi, A., J. A. Jackson, D. McKenzie, and K. Priestley, Earthquake fo- cal depths, effective elastic thickness, and the strength of the continental lithosphere, Geology, 28, 495–498, 2000. 7 Bibliography Pérez-Gussinyé, M., A. Lowry, A. B. Watts, and I. Velicogna, On the recovery of effective elastic thickness using spectral methods: Examples from synthetic data and from the Fennoscandian Shield, J. Geophys. Res., 109, B10,409, doi:10.1029/2003JB002,788, 2004. Ranalli, G., Rheology of the Earth, 2nd ed., Chapman and Hall, London, UK, 1995. Stern, R. J., Subduction zones, Rev. Geophys., 40, 1012, doi:10.1029/2001RG000,108, 2002. Thomas, W. A., Tectonic inheritance at a continental margin, GSA Today, 16, doi:10.1130/1052–5173, 2006. Turcotte, D. L., and G. Schubert, Geodynamics, Cambridge Univ. Press, Cambridge, UK, 2002. Watts, A. B., Isostasy and Flexure of the Lithosphere, Cambridge Univ. Press, Cambridge, UK, 2001. 8 Chapter 2 Wavelet analysis of the coherence between Bouguer gravity and topography: Application to the elastic thickness of the lithosphere 2.1 Introduction The correlations between the topography and gravity anomalies provide im- portant information on the level of isostatic compensation of the lithosphere at the geological time scale, and reflect its thermo-mechanical state [Watts, 2001]. The response of the lithosphere to surface (e.g. mountain belts, sedimentary basins) and internal loadings (e.g. Moho undulations) is mod- eled by assuming that regional isostasy is achieved by the flexure of a thin elastic plate overlying an inviscid fluid. The effective elastic thickness Te of the lithosphere is defined as the thickness of an equivalent elastic plate that would give the same response under the known tectonic loading. It is obtained from the flexural rigidity parameter, D, used in the equation of A version of this chapter has been published. P. Audet, and J.-C. Mareschal. Wavelet analysis of the coherence between Bouguer gravity and topography: Application to the elastic thickness anisotropy in the Canadian Shield. Geophysical Journal International, 168, 287-298, 2007. c©2007 Blackwell Publishing 9 Chapter 2. Wavelet analysis of Bouguer gravity and topography flexure of a thin elastic plate D = ET 3e 12(1 − ν2) , (2.1) where E is the Young modulus and ν Poisson’s ratio (assumed to be 1011 Pa and 0.25 throughout this work). Elastic thickness depends on many fac- tors, including the density structure, the thermal and stress state, and the mechanical properties of the lithosphere [Burov and Diament , 1995]. In the oceans, the rheology is relatively simple and the estimates of the flexural rigidity correlate well with the depth to the 450-600 oC isotherm, calcu- lated from a cooling plate model [Watts, 2001]. In the continents, where the rheological properties of the lithosphere are vertically and laterally hetero- geneous, Te does not correspond to an isotherm or to a physical boundary. In general, Te is low in young and tectonically active regions and increases in the stable continent [Lowry and Smith, 1994; Flück et al., 2003]. An ap- proximate correlation between the long term strength of the lithosphere and its age has been suggested in a few regions, e.g. in Europe [Poudjom Djo- mani et al., 1999; Pérez-Gussinyé and Watts, 2005] , and Australia [Simons and van der Hilst , 2002]. This simple model does not apply everywhere: for example, little correlation of Te with heat flow or the geology is found in some shields, e.g. in the Siberian craton [Poudjom-Djomani et al., 2003] and the Canadian Shield [Audet and Mareschal , 2004a] . There are two main approaches in the estimation of Te: the direct and inverse approaches. In the former, forward modeling of the gravity anoma- lies computed from the assumed tectonic loading can be compared with the observed gravity field to infer the mechanical properties of the litho- sphere. This method is useful for certain geological settings such as moun- tain belts, seamounts and sedimentary basins, where the loading structure is well known [e.g. Karner et al., 1983; Stewart and Watts, 1997]. More of- ten, however, the estimates of Te are inverted from the spectral relationships between topography and gravity anomalies, assuming that the lithosphere behaves as a thin elastic plate. Following Forsyth (1985), the majority of researchers use the coherence between Bouguer gravity and topography to 10 Chapter 2. Wavelet analysis of Bouguer gravity and topography estimate Te because it allows the decomposition of the loads into surface and internal components and is less sensitive to short wavelength noise in the data than the admittance between free air gravity and topography. The coherence between two fields F and G is defined in the Fourier transform domain as: γ20(k) = | < FG∗ > |2 < FF ∗ >< GG∗ > , (2.2) where <> denotes some averaging in the two-dimensional (2D) wavevector space k, and the asterisk indicates complex conjugation. The coherence is a measure of the phase relationships between two fields. If no averaging were done, the coherence would always be 1. The most common averaging method consists of binning over different annuli of wavenumber bands, but this destroys the azimuthal information in the spectra. For a discussion on the different ways of averaging the spectra, see Simons et al. [2000] . For uncorrelated fields, the phases of the cross-spectra at a given wavevector are randomly distributed and averaging cancels the coherence. For cor- related fields, the phases of the cross-spectra interfere constructively and averaging yields a high coherence. When surface or internal loads are fully compensated by the deflection of the plate at long wavelengths, the Bouguer anomaly is negatively correlated with the topography, resulting in a high co- herence. At shorter wavelengths, the loads are supported by the strength of the plate, and, if there is no correlation between initial surface and internal loading, the Bouguer anomaly is incoherent with the topography. The tran- sition wavelength from low to high coherence depends both on the rigidity and the loading structure of the plate. For estimating Te, the observed co- herence is compared with the coherence predicted for a thin elastic plate with some assumptions about the loading scheme [Forsyth, 1985]. Major differences in the estimation of Te result from the use of various spectral estimators when calculating the isotropic coherence. The fact that high flexural rigidity implies long transition wavelengths imposes a lower limit on the size of the windows used to calculate the spectra. Tests with synthetic data have shown that the Te estimates based on the modified (win- dowed or mirrored) Fourier periodogram and multitaper methods are highly 11 Chapter 2. Wavelet analysis of Bouguer gravity and topography sensitive to window size [Ojeda and Whitman, 2002; Audet and Mareschal , 2004a; Pérez-Gussinyé et al., 2004] . The multitaper method calculates the spectra with multiple orthogonal windows used as data tapers to reduce the variance of the estimates, and averaging is done over different, (ap- proximately) independent, subsets of the data [Simons et al., 2000] . The resolution degrades with the number of tapers used. Pérez-Gussinyé et al. [2004] showed that serious discrepancies occur when comparing the multi- taper coherence with theoretical curves, due to the large bias introduced near the transition wavelength by the tapering procedure. Parametric spec- tral estimators (e.g., maximum entropy) have been found to perform better on synthetic data [Lowry and Smith, 1994; Audet and Mareschal , 2004a] . However, the basic assumption of the parametric estimation, that the data were produced by an auto-regressive stochastic process, while likely to be met by numerically generated fractal surfaces, remains to be verified in the case of heterogeneous and anisotropic data, such as continental topography and gravity fields. The flexural rigidity of the lithosphere is usually assumed to be isotropic. This is a convenient assumption because it allows the reduction of the prob- lem to one dimension (1D) by averaging the azimuthal information in the spectra. It has been shown by several recent studies [e.g. Lowry and Smith, 1995; Simons et al., 2000, 2003; Rajesh et al., 2003; Swain and Kirby , 2003b; Audet and Mareschal , 2004b] that the coherence increases in one direction compared to the azimuthal average, and this anisotropy reflects the preferred direction of isostatic compensation where the lithosphere is weaker. The re- trieval of anisotropy in the coherence is hampered by the lack of a suitable averaging method. Multiple windowing techniques are better suited to this task than either the modified periodogram or the maximum entropy method because the coherence can be calculated at each wavevector, enabling the detection of anisotropy. However, all the methods mentioned above return a single estimate of Te or its anisotropy at each window, thus limiting the spatial resolution. In the course of quantifying the lateral and azimuthal variations of the coherence between topography and Bouguer anomaly, the knowledge of the 12 Chapter 2. Wavelet analysis of Bouguer gravity and topography local phase information of these fields is necessary. A multiple windowing technique that uses orthonormal Hermite polynomials in 2D, providing the necessary condition for both spatial and spectral localizations and enabling the retrieval of anisotropy in the coherence, has already been developed [Simons et al., 2003] . However, the use of many windows decreases the spectral resolution and direct comparison of the coherence with predictions is inaccurate, unless the same bias is carried in the calculation of the predicted coherence [Pérez-Gussinyé et al., 2004] . In recent years, there have been new developments in the application of wavelet analysis to geophysical data. The main advantage of the wavelet analysis over the Fourier transform approach is that it is applicable to non- stationary time series. We believe that the study of isostasy would benefit from a wavelet point of view since the wavelet transform uses optimally sized windows to obtain the spatial localization and does not require multi- windowing. So far, the majority of applications using the wavelet trans- form in isostatic studies have been restricted to the calculation of the local isotropic coherence. Kido et al. [2003] constructed an isotropic wavelet- like kernel by azimuthal averaging a 2D Gabor function and correcting for spherical geometry. Stark et al. [2003] derived a tensor-like wavelet that makes use of the multiple derivatives of the real valued Gaussian function in different directions. Recently, Kirby [2005] showed that only the complex valued Morlet wavelet and its relatives are able to reproduce the Fourier spectrum accurately because of the similarity between the basis functions of the Fourier transform and the Morlet wavelet. This property is important for the comparison of the observed 1D coherence curves with theoretical predictions to yield Te estimates. Moreover, the directional selectivity of the Morlet wavelet in the spatial domain naturally allows for the detection of azimuthal information in the spectra. Kirby [2005] developed a quasi-isotropic wavelet dubbed the ”fan” wavelet by superposing optimally spaced Morlet wavelets in a given azimuthal range. This method is further exploited in a recent paper by Kirby and Swain [2006] where they use a restricted fan wavelet to detect anisotropy in the coher- ence. This paper will first describe how the wavelet transform based on this 13 Chapter 2. Wavelet analysis of Bouguer gravity and topography directional wavelet can be used in the study of lateral and azimuthal varia- tions in the 2D coherence. The method is then demonstrated on synthetic gravity/topography with homogeneous, heterogeneous isotropic, and with homogeneous anisotropic elastic thickness. We apply this method to data from the Canadian Shield and compare the results with relevant geological and geophysical information from this area. 2.2 Wavelet analysis 2.2.1 Wavelet transform The literature on wavelets being quite extensive, we will only briefly re- view here the basic theory necessary for our particular application and refer the interested reader to Foufoula-Georgiou and Kumar [1994], Torrence and Compo [1998], or Holschneider [1999]. The 2D wavelet transform is defined as the convolution of a signal with a scaled and rotated kernel ψθa,b called a wavelet [Foufoula-Georgiou and Kumar , 1994]. Manipulating the wavelet dilation parameter (or scale) a and azimuth θ has the effect of analyzing a signal at location b at different scales and directions, where a can be re- lated to an equivalent Fourier wavenumber. The wavelet transform is thus a multi-resolution operation that allows the detection of non-stationarity and a localized characterization of the signal at different scales and azimuths. The shape of the wavelet determines its localization in both space and wavenum- ber domains. We shall use the more general terminology physical space and spectral space for the representation of a signal in the two reciprocal domains. A wavelet must satisfy two conditions: (1) zero mean, to insure a wave- like behavior, and (2) compact support in both physical and spectral spaces. We construct the family of wavelets by translating and dilating the argument of a mother wavelet ψθa,b(r) = 1 a ψ ( 1 a Cθ(r− b) ) , (2.3) where r is the location in the 2D physical plane (rx, ry), and Cθ is the 14 Chapter 2. Wavelet analysis of Bouguer gravity and topography rotation operator Cθ(r) = (rx cos(θ) + ry sin(θ),−rx sin(θ) + ry cos(θ)), (2.4) 0 ≤ θ < 2π. There are two classes of wavelet transform: (1) orthogonal and (2) non- orthogonal. The orthogonal wavelet transform allows the decomposition of a signal into a set of orthogonal basis functions, like a Fourier transform, and its validity is restricted to discrete signals. This property is appealing for the archiving of large data sets, digital image compression and filtering. It also provides a minimally redundant representation of a signal. The orthogonal wavelet kernels are usually complicated functions to express analytically in both spaces. The non-orthogonal wavelet transform is defined at arbitrary wavelengths, a useful property for the spectral analysis of continuous signals. In this study we use a discrete version of the non-orthogonal continuous wavelet transform (conventionally termed CWT) Wf of f(r), defined as Wf (a, θ,b) = ∫ R2 f(r)ψθ∗a,b(r), d 2r, (2.5) where the asterisk denotes complex conjugation. The kernel functions ψθa,b determine the resolution of the wavelet transform in both physical and spec- tral spaces. At large scales, the wavelet is spread in physical space and the wavelet transform picks up information at long wavelengths. At smaller scales, the wavelet is well localized in physical space and the wavelet trans- form picks up information at short wavelengths. The wavelet resolution is bounded by the uncertainty principle, which states that physical and spec- tral localizations cannot be simultaneously measured with arbitrarily high precision. Hence the spatial localization of the wavelet comes at the expense of the spectral localization. In Fourier domain, the wavelet transform ŴF becomes ŴF (a, θ,k) = F̂ (k)ψ̂ ∗ a(Cθ(k)), (2.6) 15 Chapter 2. Wavelet analysis of Bouguer gravity and topography where k is the two dimensional wavevector (kx, ky), F̂ (k) is the 2D Fourier transform of f(r) and ψ̂∗a(Cθ(k)) is the complex conjugate of the Fourier transform of the wavelet ψ at scale a and azimuth θ. The multiplication of the spectra in the Fourier domain is faster to perform than the convolution in equation (2.5), and one can take advantage of the Fast Fourier Transform (FFT) algorithm for transforming back and forth from space to Fourier domain. 2.2.2 Wavelets The choice of a wavelet depends on several factors, including the need for a real or complex valued function, shape and resolution of the desired wavelet, and ultimately the capacity of resolving anisotropic features in 2D fields. In our case where we want to calculate the coherence between two fields, it is important to use a complex valued wavelet to retrieve the phase information. To compare the coherence with Fourier-derived predictions, it is also natural to use a wavelet constructed from complex exponentials modulated by a smoothing function. The Morlet wavelet is thus an ideal candidate for this task and exhibits a natural property of directional selectivity which allows the detection of azimuthal information in 2D fields. Morlet wavelet In 2D, the Morlet wavelet is an oscillating function of wavevector k0 = (k0x, k 0 y) modulated by a Gaussian smoothing function ψ(r) = e−ik0·re− 1 2 |Ar|2, (2.7) where |k0| ≈ 5.336 in order to satisfy the zero mean condition, and A is the 2 × 2 anisotropic diagonal matrix where the first non-zero element is ǫ, with 0 < ǫ ≤ 1, and the second element is 1 [Antoine et al., 1993; Kumar , 1995; Antoine et al., 1996] . By setting ǫ = 1, the Morlet wavelet has an isotropic envelope in the spectral space. This property will be convenient in the construction of the fan wavelet to be described later. We refer to this 16 Chapter 2. Wavelet analysis of Bouguer gravity and topography wavelet as the isotropic Morlet wavelet. In the spectral space the Morlet wavelet is a Gaussian band-pass filter centered on the wavevector k0 ψ̂(k) = e− 1 2 |k−k0| 2 . (2.8) The Morlet wavelet in the spectral space is almost entirely supported in the positive domain, subject to the restriction |k| > |k0|. The wavevector k0 can be arbitrarily chosen and will filter information in the direction given by tan(θ0) = k0y k0x . In Fourier domain, the peak representation of the Morlet wavelet is at wavenumber kF = |k0| a , (2.9) which is considered as the equivalent Fourier wavenumber. “Fan” wavelet In the context of finding a suitable isotropic wavelet that can be interpreted as a Fourier spectrum, Kirby [2005] showed that a controlled superposition of directionally adjacent isotropic Morlet wavelets, dubbed the “fan” wavelet because of its shape in spectral space, provided better results than any of the commonly used wavelets (DoG (Derivative of Gaussian), Paul, Perrier and Poisson). The fan wavelet is expressed as ψ̂F (k) = 1 Nθ Nθ∑ i=1 ψ̂(Cθi(k)). (2.10) The geometry of the fan wavelet is achieved by averaging adjacent isotropic Morlet wavelets separated by an azimuthal sampling of δθ = 2 √−2 ln p/|k0|, where the optimal value of p was found to be 0.75 [Kirby , 2005]. The fan wavelet is obtained by setting a single range of azimuthal increments Nθ = int(∆θ/δθ), where ∆θ is the azimuthal extent. For ∆θ = π, the fan wavelet includes 11 adjacent Morlet wavelets and is quasi-isotropic. For 17 Chapter 2. Wavelet analysis of Bouguer gravity and topography smaller values of ∆θ, the fan wavelet is anisotropic. 2.2.3 Cross spectral analysis Equation (2.5) is used to compute a wavelet scalogram, which represents the energy density of a function f(r) at scale a, location b and direction θ Sff (a, θ,b) = |Wf (a, θ,b)|2. (2.11) The wavelet cross-scalogram of any aritrary 2D functions f and g is calculated in the same way as (2.11) Sfg(a, θ,b) =Wf (a, θ,b) ·W ∗g (a, θ,b). (2.12) Equations (2.11) and (2.12) are combined to estimate the local wavelet com- plex admittance Qfg(a, θ,b) = Sfg(a, θ,b) Sff (a, θ,b) (2.13) and the complex impedance Q ′ fg(a, θ,b) = Qgf (a, θ,b). The local wavelet coherence is the product of admittance with impedance: γ2fg(a, θ,b) = |Sfg(a, θ,b)|2 Sff (a, θ,b)Sgg(a, θ,b) . (2.14) At this point, if no averaging is performed, the coherence is 1 at all scales and azimuths. In the wavelet implementation using the fan wavelet, the averaging is done over distinct azimuths γ2fg(a, θ,b) = < |Sfg(a, θ,b)| >2θ′ < Sff (a, θ,b) >θ′< Sgg(a, θ,b) >θ′ , (2.15) where θ′ = θ ± ∆θ. Details on the azimuthal averaging are discussed in the next section. One can also calculate a global estimate of the coherence from the wavelet scalograms by averaging the local wavelet spectra over the 18 Chapter 2. Wavelet analysis of Bouguer gravity and topography spatial coordinates b at each scale and azimuth. γ̄2fg(a, θ) = | < Sfg(a, θ,b) >b |2 < Sff (a, θ,b) >b< Sgg(a, θ,b) >b . (2.16) 2.2.4 Algorithm We follow the standard practice and calculate the convolution in Fourier domain using the fast Fourier transform (FFT) algorithm to transform back and forth from space to spectral domains. We now describe in more de- tails the steps involved in the calculation of the cross-spectral quantities (equations 2.11 to 2.14). We start by Fourier transforming the Bouguer and topography fields. In the second step, we define (1) the range of scales to be computed, depend- ing on the size of the data window and the size of the grid, and (2) the discrete angle increments for the azimuthal selectivity of the fan wavelet (≈ 16.3o). Second we loop over all scales and azimuths up to 180o and (i) cal- culate the daughter wavelet in spectral space for this set of parameters; (ii) compute the Fourier representation of the wavelet transform by multiplying the daughter wavelet with the transformed fields (equation 2.6); (iii) inverse Fourier transform the result to obtain the local wavelet transform in phys- ical space (equation 2.5); (iv) calculate the scalograms and cross-scalogram (equations 2.11 and 2.12); (v) average the scalograms and cross- scalogram depending on the type of analysis; and finally (vi) combine the resulting quantities to obtain the local admittance and coherence (equations 2.13 and 2.14). Step (v) is the key procedure in the algorithm, as it controls the az- imuthal resolution and the variance of the coherence. In the isotropic case, the averaging is performed over individual wavelet spectra calculated from a range of Morlet wavelets spanning 180o in the spectral space. By aver- aging the local spectra over a restricted range of azimuths, the anisotropic coherence is obtained at the azimuth corresponding to the middling value of the range. The width of the range controls the variance by virtue of the number of selected spectra in the averaging, while the overlap between 19 Chapter 2. Wavelet analysis of Bouguer gravity and topography adjacent ranges controls the azimuthal resolution. Note that no averaging in wavenumber (or scale) is performed. We follow Kirby and Swain [2006] and set the azimuthal range to 90o in the anisotropic analysis. The az- imuthal distance between overlapping ranges is set to 5o, which is smaller than the azimuthal separation between adjacent Morlet wavelets.This in- troduces some redundancy and increases computational time but improves upon the azimuthal resolution in the inversion. 2.2.5 Limitations While the use of the FFT to compute the convolution in the spectral domain speeds up the calculations, it can also bias the scalograms at the larger scales. Because we are dealing with finite size 2D fields, errors will occur at the edges of the grids, as the FFT assumes that the data are periodic. This means that signals at one edge of the field will get wrapped around the other end [Torrence and Compo, 1998]. This effect is more pronounced at larger scales as the influence of each wavelet extends further in space. There is thus a cone of influence beyond which edge effects become important. In our application, it means that for grid points close to any edge of the grid which have a large transition wavelength in the coherence (large Te), the estimated elastic thickness and the anisotropy are affected by data from the opposite side of the grid. This effect increases with increasing contrast in Te at opposing edges of the study area. There exists many pre-processing techniques aimed at removing the edge effect in spectral analysis. The most obvious one that does not require any data manipulation is to select a data window that is much wider than the region of interest with the smallest grid sample possible. This procedure is ideal but not applicable in all cases as the data availability and resolution are usually limited. Another possibility is to mirror the edges of the win- dow, calculate the spectra on the larger data set, and retain only the initial quadrant in the final analysis. However it is well known that this procedure introduces artificially long wavelengths and can greatly affect the value of Te [Lowry and Smith, 1994; Ojeda and Whitman, 2002; Swain and Kirby , 20 Chapter 2. Wavelet analysis of Bouguer gravity and topography 2006]. For example, small wavelength coherent signals at one edge of the initial data window will be mapped at longer wavelengths in the enlarged grid due to the mirroring, and the coherence curve will be shifted accord- ingly, resulting in the elastic thickness being over-estimated at the edges. The anisotropic coherence might suffer even more from the mirroring since the weak direction will also be mirrored at grid points close to any edge of the data window. The preferred method used in this study pads the 2D fields with ze- ros to the next power of 2 in each direction, applies the FFT on the larger grid, and retains only the relevant quadrant in the computation of the scalo- grams. Spectral leakage due to discontinuities in the data are dealt with by the application of a Hanning taper to all sides of the initial data set. This procedure down-weights the wavelet spectral power at the edges, but the phase of the cross-spectra is confidently recovered and no artificial informa- tion is introduced in the data. 2.2.6 Isotropic coherence inversion In the inversion step, Te is estimated by minimizing the misfit between the observed coherence and predictions from an ideal elastic lithosphere, pro- vided a certain amount of assumptions about the loading scheme, depth to the surface of compensation, and the elastic parameters. One approach compares the observed coherence with theoretical (analytical) curves, that is, with the solution of the elastic plate for an infinitely large window and a uniform load ratio f [Kirby and Swain, 2004]. However, as outlined by Pérez-Gussinyé et al. [2004] and Simons et al. [2000], there is a bias in the spectral estimates of the coherence. The comparison of biased estimates with theoretical (unbiased) predictions results in an inaccurate estimate for the elastic thickness. The wavelet transform does not escape this problem because the resolution in physical and spectral spaces is bounded by the un- certainty principle. Wavelets are band-pass functions, hence the transform at a single scale will include information from neighboring wavenumbers. As the scales increase, the wavelets are more localized in the spectral space 21 Chapter 2. Wavelet analysis of Bouguer gravity and topography and the transform is closer to the Fourier spectrum, at the cost of spatial resolution. At wavenumbers near the transition from low to high coherence, the estimated coherence is biased towards higher values because the contri- bution to the coherence of the wavelet at neighboring lower wavenumbers will smear out the long wavelength information. The fractal nature of the topography and Bouguer gravity fields will also put more weight on the long wavelengths and together these effects will shift the transition wavelength toward lower values, resulting in the elastic thickness being underestimated. The assumption of uniform load ratio f is also a crude approximation be- cause it is a wavenumber dependent parameter calculated from the data that affects the shape of the coherence, albeit in a weaker manner than Te [Forsyth, 1985]. Pérez-Gussinyé et al. [2004] suggest to use the same spectral estimator for calculating the predicted and the observed coherence. In this case, the same biases affect both quantities and their effects balance each other in the inversion for the elastic thickness. For the wavelet spectral estimator, the Fourier-based deconvolution method of Forsyth [1985] can be adapted to fit the purpose of Te calculation by assuming that adjacent local spectra are decoupled, as was done by Stark et al. [2003] and Swain and Kirby [2006]. The method proceeds as follows: (1) using the wavelet transforms of the Bouguer gravity WB and topography WH at each scale and azimuth, we select a particular spatial point and solve for the initial surface WHi and internal WBi loads for a given value of Te; (2) the initial loads are decomposed into the four components of the surface (WHt and WBt) and subsurface (WHb and WBb) deflection that they cause; and (3) the auto and cross spectra are averaged and combined to form the local predicted coherence, assuming that the initial loads are statistically uncorrelated: γ2p = | <WHtW ∗Bt +WHbW ∗Bb>θ |2 (<WHtW ∗ Ht >θ + <WHbW ∗ Hb >θ)(<WBtW ∗ Bt >θ + <WBbW ∗ Bb >θ) , (2.17) where the brackets indicate averaging over all azimuths. Thus by minimizing the misfit between the observed and predicted coherence, we are also able 22 Chapter 2. Wavelet analysis of Bouguer gravity and topography to calculate the local subsurface to surface loading ratio f , which is given by f2 = (ρm − ρc)2 ρ2c < WBiW ∗ Bi >θ < WHiW ∗ Hi >θ , (2.18) where WBi and WHi are the components of the best fit solution, and ρc and ρm are the crust and mantle density, respectively. Some authors have recently suggested to use the transition wavelength from low to high coherence (typically a coherence of 0.5) as a proxy for the mechanical strength, because it does not require making many uncer- tain assumptions as those needed for inverting Te, such as the depth and magnitude of the density contrast [Simons and van der Hilst , 2003; Rajesh and Mishra, 2004]. The inversion of the observed coherence from synthetic data using theoretical curves should then be equivalent to calculating the transition wavelength since the density structure and the uniform load ratio are known. It is expected that both inversion methods should yield similar results on synthetic data when the transition wavelength is much smaller than the window size. At longer transition wavelengths (higher Te) the lim- ited wavelet resolution introduces bias in the estimated coherence near the transition wavelength and the Fourier deconvolution method should be more accurate than the inversion based on theoretical solutions. In the inversion, the misfit is calculated by the L2 norm weighted by the variance of the observed coherence. In the next section we compare the abil- ity of the different inversion methods torecover Te through the application of the wavelet method on synthetic grids. 2.2.7 Anisotropic coherence inversion The full inversion of the 2D coherence for anisotropic Te necessitates the modeling of the lithosphere as a thin plate with anisotropic mechanical properties. Although this method has already been developed and used in isostatic studies [Swain and Kirby , 2003b; Kirby and Swain, 2006] , its application in a laterally varying lithosphere remains ambiguous because of the complexity observed in the 2D coherence. The method compares the 2D 23 Chapter 2. Wavelet analysis of Bouguer gravity and topography wavelet coherence with theoretical solutions to the equation of the flexure of an orthotropic elastic plate given a large number of assumptions, such as intrinsic orthotropy, and isotropy in the subsurface to surface load ratio, which is kept constant at all wavelengths. The discussion about the bias in the estimates of the wavelet coherence and how it affects the recovery of Te also applies in this case. These arguments, coupled with the observed com- plexity in the observed wavelet coherence, render the interpretation in terms of the elastic thickness in two perpendicular directions debatable. However the response of the lithosphere should in theory be close to that of a thin elastic plate, and the orthotropic elastic plate model has the advantage of determining physical estimates of the anisotropy, whether the assumptions are valid or not. We thus choose to invert the 2D coherence using the the- oretical solutions as in Kirby and Swain [2006] . Relative comparisons can be made between the recovered isotropic and anisotropic elastic thicknesses, while the estimates of the weak direction are regarded as absolute. In the 2D coherence inversion, the misfit is weighted by the inverse of wavenumber to down-weight the contribution of short wavelength coherent features that are not associated with the response of an elastic plate [Kirby and Swain, 2006]. 2.3 Numerical examples Synthetic gravity and topography were generated as two non- correlated fractal fields with a spectral slope of β = −3 representing the initial surface and subsurface loads that are applied to a thin elastic plate with known elastic parameters and density structure [Macario et al., 1995]. The elastic plate is assumed to consist of a crust of uniform thickness and density. The subsurface load is placed at the Moho 40 km deep, with a 500 kg/m3 density contrast at the interface and a crustal density of 2670 kg/m3. The ampli- tudes of the loads were scaled by imposing a constant ratio f = 1 between surface and subsurface loads. The flexure produced by both loads was cal- culated by solving the thin elastic plate equilibrium equation in the spectral domain for uniform rigidity plates and with a finite-difference solution in 24 Chapter 2. Wavelet analysis of Bouguer gravity and topography the spatial domain for the case with a varying rigidity pattern. 2.3.1 Uniform rigidity pattern We first verified the robustness of the wavelet transform method on syn- thetic data with spatially uniform flexural rigidity. We treat the isotropic and anisotropic cases independently and compare the ability of the wavelet transform and the different inversion methods in recovering the input Te. In all examples, synthetic data were generated on 256 × 256 grids with a sampling interval of 8 km. The spectra were calculated on a smaller sub-grid (128 × 128) to avoid using periodic data sets. Isotropic plate We first tested the homogeneous isotropic plate. Applying the procedure outlined above, we generated 100 data sets for each isotropic elastic plate thickness of 10, 20, 40, and 80 km. Various tests were performed on these synthetic grids. First we compared the global wavelet coherence with the- oretical predictions to evaluate the ability of the wavelet transform to re- produce the Fourier coherence. The theoretical curves are calculated from equation (10) of Pérez-Gussinyé et al. [2004], corrected for the numerator (variable φ not squared). We also compared the global wavelet coherence with the predicted global wavelet coherence as described in section 2.6. The inversion for the best fitting elastic thickness is effectively done for each of the 100 runs and each input value of Te, a procedure that is now rou- tinely carried out in studies of the recovery of Te with spectral methods [e.g. Macario et al., 1995; Swain and Kirby , 2003a; Pérez-Gussinyé et al., 2004; Audet and Mareschal , 2004a]. In Fig. 2.1 the histograms show the distributions, mean values and standard deviations for each input Te for both inversion methods. The relative difference between the observed and predicted solutions and the standard deviations increase when Te is > 40 km. When the observed wavelet coherence is compared with the predicted wavelet coherence rather than the theoretical function, the solutions are closer to the input values except for the lowest Te, where the estimates are 25 Chapter 2. Wavelet analysis of Bouguer gravity and topography similar. This could be because the theoretical solutions are as good as the deconvolution solutions for low Te’s, where the transition wavelength is much less than the area dimensions, as argued before. We also calculated the local coherence at each spatial point for a sin- gle synthetic data set. Fig. 2.2 shows example results with both inversion methods for input Te = 20 km: (a-b) maps of the Te variations; (c) observed global isotropic coherence curve and the best-fitting predicted and theoret- ical coherence. The Te recovered from the global wavelet coherence is 19.0 km when Te is inverted from the theoretical curves, and the mean and stan- dard deviations of the local estimates on the whole grid is 19.5 ± 3.7 km; estimates for Te are respectively 22.0 km for the global coherence, and 24.5± 4.1 km for the mean and standard deviations when Te is inverted from the predicted coherence. The dominant wavelength in the spatial variations of the recovered Te is about 500 km, similar to what was found by Kirby and Swain [2004]. With controlled synthetic data, the spatial variations in the estimated Te are on the order of 15% the input value. Anisotropic plate Anisotropy in the flexural rigidity parameters Dx and Dy, where x and y are two perpendicular directions in a Cartesian coordinate system, produces a decrease in the transition wavelength in the direction where the rigidity is lowest. The rigidities are converted to equivalent elastic thicknesses in two perpendicular directions using constant and isotropic Poisson’s ratio andY- oung’s modulus. Following Swain and Kirby [2003b] method for the genera- tion of synthetic gravity and topography with homogeneous anisotropic elas- tic thickness, and using the same procedure as before we generated synthetic data sets with uniform elastic thickness pairs of (Tex, Tey, φ) = (5, 10, 0), (20, 10, 0), (50, 10, 0) km, where φ is the weak direction calculated anti- clockwise from the x-axis. The global wavelet anisotropic coherence inver- sion gives (Tex, Tey, φ) = (5, 10, 0), (18, 10, 5), (40, 12, 0). Fig. 2.3 shows the wavelet coherence at each scale and azimuth and the best fit theoretical coherence at a randomly chosen point on the grid. The local coherence in- 26 Chapter 2. Wavelet analysis of Bouguer gravity and topography 0 10 20 30 40 50 60 Pe rc en t 0 5 10 15 20 25 30 Te = 10 km <Te> = 9.4 ± 0.3 km a 0 10 20 30 40 50 Pe rc en t 0 20 40 60 80 Te = 20 km <Te> = 18.5 ± 1.1 km b 0 10 20 30 40 50 Pe rc en t 0 30 60 90 120 Te = 40 km <Te> = 33.9 ± 5.7 km c 0 10 20 30 40 Pe rc en t 0 30 60 90 120 Estimated Te Te = 80 km <Te> = 59.1 ± 14.6 km d 0 10 20 30 40 50 60 Pe rc en t 0 5 10 15 20 25 30 Te = 10 km <Te> = 11.3 ± 1.7 km e 0 10 20 30 40 50 Pe rc en t 0 20 40 60 80 Te = 20 km <Te> = 19.6 ± 3.1 km f 0 10 20 30 40 50 Pe rc en t 0 30 60 90 120 Te = 40 km <Te> = 40.9 ± 8.0 km g 0 10 20 30 40 Pe rc en t 0 30 60 90 120 Estimated Te Te = 80 km <Te> = 64.3 ± 13.9 km h Figure 2.1: Results from synthetic data with isotropic and spatially uniform Te for four different values of Te: (a,e) 10 km, (b,f) 20 km, (c,g) 40 km, and (d,h) 80 km. The observed global coherence is calculated using the isotropic fan wavelet. (a-h) distributions of recovered Te from the inversion of the coherence using the theoretical (a-d), and the deconvolution solutions (e-h). 27 Chapter 2. Wavelet analysis of Bouguer gravity and topography 20 200 400 600 800 1000 200 400 600 8001000 a 200 400 600 8001000 b 10 12 14 16 18 20 22 24 26 Elastic Thickness (km) 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e 10-2 10-1 Wavenumber (1/km) c Figure 2.2: Example results on synthetic data with uniform Te = 20 km. (a-b) maps of the Te variations calculated on a coarser grid (1/16 of total grid points). (a) is obtained with the inversion from theoretical curves; (b) is inverted from the predicted coherences. (c) Observed (filled circles), theoretical (dotted line), and predicted (solid line) global coherence curves, version results are sensibly the same as for the global results and show that we can recover the direction of highest coherence with a good accuracy (± 10o). For large Te contrasts, the recovered Te values are inaccurate but the coherence anisotropy is well marked (Fig 2.3c,f). We must also keep in mind that the inversion with theoretical solutions underestimates Te, particularly when Te >40 km for this size of grid. 2.3.2 Spatially varying rigidity pattern The motivation for choosing a wavelet method over traditional spectral methods is an improved spatial resolution for recovering variations in Te. To address the issue of resolution, we calculated the local wavelet coherence on synthetic Bouguer and topography grids with a Te structure consisting of a circular core with diameter of 700 km and Te = 50 km, with radially decreasing values towards the edges of the window to Te = 0 km [Pérez- Gussinyé et al., 2004]. The data were generated with a uniform load ratio of f = 1, and the depth to Moho = 40 km. The study area size is 2048×2048 km2, with a sampling interval of 8 km. In this case the whole grid is retained 28 Chapter 2. Wavelet analysis of Bouguer gravity and topography 1722 33 49 73 109 163 242 360 536 798 1187 1766 W av e le ng th (k m) 1722 33 49 73 109 163 242 360 536 798 1187 1766 W av e le ng th (k m) 30 60 90 120150180 Azimuth (deg) 1722 33 49 73 109 163 242 360 536 798 1187 1766 W av e le ng th (k m) 1722 33 49 73 109 163 242 360 536 798 1187 1766 1722 33 49 73 109 163 242 360 536 798 1187 1766 30 60 90 120150180 Azimuth (deg) 1722 33 49 73 109 163 242 360 536 798 1187 1766 a d eb fc 0.0 0.2 0.4 0.6 0.8 1.0 Coherence Figure 2.3: Results of the wavelet method using the directional Morlet wavelet on synthetic data with spatially uniform and anisotropic Te. The azimuth is measured anticlockwise from the x axis. (a-c) observed wavelet coherence; (d-f) best fit theoretical coherence. (a,d) Tex = 5, Tey = 10 km; (b,e) Tex = 20, Tey = 10 km; and (c,f) Tex = 50, Tey = 10 km. 29 Chapter 2. Wavelet analysis of Bouguer gravity and topography 20 20 30 400 800 1200 1600 2000 400 800 120016002000 a 20 30 40 400 800 120016002000 b 0 10 20 30 40 50 Elastic Thickness (km) 10 10 20 3040 50 400 800 120016002000 c Figure 2.4: Recovered Te pattern from the input model shown in (c). (a) is obtained by inverting the local coherence with theoretical curves. Esti- mates of (b) are calculated by inverting the local coherence with the adapted deconvolution method described in section 2.2.6. in the calculations and the zero-padding/tapering procedure is not applied since the data are periodic in both directions. The results in Fig. 2.4 show the recovered Te pattern obtained with both inversion methods along with the input structure. Both inversion methods show a roughly circular Te structure with values increasing from ≈ 0 km at the edges to ≈ 50 km towards the center. This figure should be compared with Fig. 4a-h of Pérez-Gussinyé et al. [2004] . While their coherence results show a similar pattern when they use a 800× 800 km2 window, the resolution is limited to a 1200 × 1200 km2 square region at the center of the grid. The wavelet method, on the other hand, is able to recover the low values at the edges with good accuracy. However this might not be the case for different Te models as discussed in section 2.5. The RMS errors calculated on both grids are 9.1 and 8.4, respectively for (a) and (b). The deconvolution method also produces higher estimates of Te in the central region of the study area. The larger RMS error than in Kirby and Swain [2004] and Stark et al. [2003] can be partly explained by the sharpness in the structure of our model, and by the wide range of input Te values. This rather extreme example illustrates well the limits of 30 Chapter 2. Wavelet analysis of Bouguer gravity and topography spectral methods in recovering spatial variations in Te: structures with high values of Te are well estimated only if their wavelength is longer than the transition wavelength from low to high coherence. When this condition is not verified, the spatial variations of Te represent a weighted average of the true distribution. In light of these results on synthetic data, it is clear that one must be careful when interpreting maps of Te values, since both the processing and the inversion steps rely on a large number of assumptions, and even with known input parameters the inversion is flawed. The accuracy of recovered Te values strongly depends on the size of the study area compared to the transition wavelength. 2.4 Conclusion We have tested and used a wavelet based method to determine the coherence between Bouguer gravity and topography and to obtain estimates of the effective elastic thickness of the lithosphere and its anisotropy. Tests with synthetic data show that effective elastic thickness can be estimated within ± 15% with the wavelet coherence. Tests also show that the wavelet method can detect anisotropy and determine the directional variations in effective elastic thickness. 2.5 Acknowledgments We are grateful to Marta Pérez-Gussinyé for kindly providing some of the synthetic data used in this study. This paper was much improved by thor- ough and constructive reviews by Jon Kirby and Marta Pérez-Gussinyé. 31 Bibliography Antoine, J.-P., P. Carrette, R. Murenzi, and B. Piette, Image analysis with two-dimensional continuous wavelet transform, Signal Processing, 31, 241– 272, 1993. Antoine, J.-P., R. Murenzi, and P. Vandergheynst, Two-dimensional direc- tional wavelets in image processing, J. Imag. Sys. Technol., 7, 152–165, 1996. Audet, P., and J.-C. Mareschal, Anisotropy of the flexural response of the lithosphere in the Canadian Shield, Geophys. Res. Lett., 31, L20,601, doi:10.1029/2004GL021,080, 2004a. Audet, P., and J.-C. Mareschal, Variations in elastic thickness in the Cana- dian Shield, Earth Planet. Sci. Lett., 226, 17–31, 2004b. Burov, E. B., and M. Diament, The effective elastic thickness (Te) of con- tinental lithosphere: What does it really mean?, J. Geophys. Res., 100, 3905–3927, 1995. Flück, P., R. D. Hyndman, and C. Lowe, Effective elastic thickness Te of the lithosphere in western Canada, J. Geophys. Res., 108, 2430, doi:10.1029/2002JB002,201, 2003. Forsyth, D. W., Subsurface loading and estimates of the flexural rigidity of continental lithosphere, J. Geophys. Res., 90, 12,623–12,632, 1985. Foufoula-Georgiou, E., and P. Kumar (Eds.),Wavelets in Geophysics, Aca- demic Press, San Diego, Calif., 1994. 32 Bibliography Holschneider, M., Wavelets: An Analysis Tool, Oxford University Press, Oxford, 1999. Karner, G. D., M. S. Steckler, and J. A. Thorne, Long-term thermo- mechanical properties of the continental lithosphere, Nature, 304, 250–253, 1983. Kido, M., D. A. Yuen, and A. P. Vincent, Continuous wavelet-like filter for a spherical surface and its application to localized admittance function on Mars, Phys. Earth Planet. Inter., 135, 1–14, 2003. Kirby, J. F., Which wavelet best reproduces the Fourier power spectrum?, Comput. Geosc, 31, 846–864, 2005. Kirby, J. F., and C. J. Swain, Global and local isostatic coherence from the wavelet transform, Geophys. Res. Lett., 31, doi:10.1029/2004GL021,569, 2004. Kirby, J. F., and C. J. Swain, Mapping the mechanical anisotropy of the lithosphere using a 2-D wavelet coherence, and its application to Australia, Phys. Earth Planet. Int., 158, 122–138, 2006. Kumar, P., A wavelet-based methodology for scale-space anisotropic anal- ysis, Geophys. Res. Lett., 22, 2777–2780, 1995. Lowry, A. R., and R. B. Smith, Flexural rigidity of the Basin and Range- Colorado Plateau-Rocky Mountain transition from coherence analysis of gravity and topography, J. Geophys. Res., 99, 20,123–20,140, 1994. Lowry, A. R., and R. B. Smith, Strength and rheology of the western U. S. Cordillera, J. Geophys. Res., 100, 17,947–17,963, 1995. Macario, A., A. Malinverno, and W. F. Haxby, On the robustness of elastic thickness estimates obtained using the coherence method, J. Geophys. Res., 100, 15,163–15,172, 1995. Ojeda, G. Y., and D. Whitman, Effect of windowing on lithosphere elas- tic thickness estimates obatined via the coherence method: Results from 33 Bibliography northern South America, JGR, 107, 2275, doi:10.1029/2000JB000,114, 2002. Pérez-Gussinyé, M., and A. B. Watts, The long-term strength of Europe and its implications for plate- forming processes, Nature, 436, 381–384, 2005. Pérez-Gussinyé, M., A. Lowry, A. B. Watts, and I. Velicogna, On the recovery of effective elastic thickness using spectral methods: Examples from synthetic data and from the Fennoscandian Shield, J. Geophys. Res., 109, B10,409, doi:10.1029/2003JB002,788, 2004. Poudjom Djomani, Y. H., J. D. Fairhead, and W. L. Griffin, The flexu- ral rigidity of Fennoscandia: Reflection of the tectonothermal age of the lithospheric mantle, Earth Planet. Sci. Lett., 174, 139–154, 1999. Poudjom-Djomani, Y. H., S. Y. O’Reilly, W. L. Griffin, L. M. Nat- apov, Y. Erinchek, and J. Hronsky, Upper mantle structure beneath eastern Siberia: Evidence from gravity modelling and mantle petrology, Geochem. Geophys. Geosys., 4, 1066, doi:10.1029/2002GC000,420., 2003. Rajesh, R. S., and D. C. Mishra, Lithospheric thickness and mechanical strength of the indian shield, Earth Planet. Sci. Lett., 225, 319–328, 2004. Rajesh, R. S., J. Stephen, and D. C. Mishra, Isostatic response and anisotropy of the Eastsern Himalayan-Tibetan Plateau: A reap- praisal using multitaper spectral analysis, Geophys. Res. Lett., 30, 1060, doi:10.1029/2002GL016,104, 2003. Simons, F. J., and R. D. van der Hilst, Age-dependent seismic thickness and mechanical strength of the Australian lithosphere, Geophys. Res. Lett., 29, 1529, doi:10.1029/2002GL014,962, 2002. Simons, F. J., and R. D. van der Hilst, Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere, Earth Planet. Sci. Lett., 211, 271–286, doi:10.1016/S0012–821X(03)00,198– 5, 2003. 34 Bibliography Simons, F. J., M. T. Zuber, and J. Korenaga, Isostatic response of the Aus- tralian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper spectral analysis, J. Geophys. Res., 105, 19,163–19,184, 2000. Simons, F. J., R. D. van der Hilst, and M. T. Zuber, Spatio-spectral lo- calization of isostatic coherence anisotropy in Australia and its relation to seismic anisotropy: Implications for lithospheric deformation, J. Geo- phys. Res., 108, 2250, doi:10.1029/2001JB000,704, 2003. Stark, C. P., J. Stewart, and C. J. Ebinger, Wavelet transform mapping of effective elastic thickness and plate loading: Validation using synthetic data and application to the study of southern African tectonics, J. Geo- phys. Res., 108, 2258, doi:10.1029/2001JB000,609, 2003. Stewart, J., and A. B. Watts, Gravity anomalies and spatial variations of flexural rigidity at mountain ranges, J. Geophys. Res., 102, 5327–5352, 1997. Swain, C. J., and J. F. Kirby, The effect of ”noise” on estimates of the elastic thickness of the continental lithopshere by the coherence method, Geophys. Res. Lett., 30, 1574, doi:10.1029/2003GL017,070, 2003a. Swain, C. J., and J. F. Kirby, The coherence method using a thin anisotropic plate model, Geophys. Res. Lett., 30, 2014, doi:10.1029/2003GL018,350, 2003b. Swain, C. J., and J. F. Kirby, An effective elastic thickness map of Australia from wavelet transforms of gravity and topography using Forsyth’s method, Geophys. Res. Lett., 33, L02,314, doi:10.1029/2005GL025,090, 2006. Torrence, C., and G. P. Compo, A practical guide to wavelet analysis, Bull. Am. Meter. Soc., 79, 61–78, 1998. Watts, A. B., Isostasy and Flexure of the Lithosphere, Cambridge Univ. Press, Cambridge, UK, 2001. 35 Chapter 3 Mechanical controls on the deformation of continents at convergent margins 3.1 Introduction At a continental convergent boundary long-term tectonic loading may de- form a lithospheric plate permanently by two mechanisms: brittle failure in the cool upper part of the lithosphere and by effectively viscous flow in the hotter deeper regions. The maximum stress that can be supported by the lithosphere is thought to occur at the intersection of these regimes, the so- called brittle-ductile transition [Kohlstedt et al., 1995]. Under a present-day stress regime brittle failure is apparent from seismicity. A number of studies have shown that the depth to which seismicity occurs (the seismogenic zone) plausibly corresponds to a lower bound to the elastic thickness, Te, which for a lithostatic state of stress is the vertical fraction of the lithosphere that behaves elastically over geological time scales [Maggi et al., 2000; Watts and Burov , 2003]. The flexural rigidity D ∝ ET 3e , where E is Young’s modulus, is a rheological property that governs the resistance of the lithosphere to flexure by plate boundary forces [Watts, 2001]. The strong dependence of D on Te implies that the magnitude and spatial variations of Te can have a significant influence on the degree and style of deformation due to long-term tectonic loads. In particular, it is expected that spatial variations in Te can A version of this chapter has been published. P. Audet, A. M. Jellinek, and H. Uno. Mechanical controls on the deformation of continents at convergent margins. Earth and Planetary Science letters, 264, 151-166, 2007. c©2007 Elsevier B.V. 36 Chapter 3. Deformation of continents at convergent margins prescribe where stress concentration may occur and consequently determine the location of earthquakes, as well as mechanical anisotropy, due to the localized brittle failure of crustal rocks [Lowry and Smith, 1995]. The mechanics governing the tectonic evolution of convergent plate bound- aries remain a major focus of both observational [Peacock , 1990; Stern, 2002] and modelling [Becker et al., 1999; Billen et al., 2003; Stegman et al., 2006] studies of subduction zones. However, although potentially important, the influences of Te variations on the rheology, dynamics and evolution of con- vergent margins are poorly understood [Burov and Diament , 1995; Gaspar- Escribano et al., 2001]. A central reason is the difficulty in constraining spatial variations in Te. Only a few studies [e.g. Flück et al., 2003; Rajesh et al., 2003] have focused on estimating Te in active convergent tectonic settings using window-based methods. More recently, Tassara et al. [2007] used the wavelet method to calculate Te in South America. In this study we calculate Te and Te anisotropy in western Canada by using and compar- ing the multitaper and wavelet transform methods. Western Canada is an ideal target for Te studies because it has a rich and diverse tectonic history spanning ∼3 Ga, where rigidity variations may have played a major role in continental evolution. This region encompasses three main tectonic do- mains: the young and tectonically active Phanerozoic Cordillera to the west, the Paleo-Proterozoic Interior Plains covered by a thick layer of Phanero- zoic sediments, and the stable Archean Craton to the east. The tectonic evolution of the plate margin has been dominated by convergence since the Mezosoic [Price, 1986]. The plate boundary at the western edge of the con- tinent is currently defined by the subduction of the Juan de Fuca plate to the south, and the strike-slip Queen Charlotte fault system to the north. To gain a sound understanding of the mechanical controls on the de- formation of the lithosphere at convergent margins, Te results should be compared with a range of geophysical indicators of lithospheric deformation [Eaton and Jones, 2006]. In western Canada several experiments conducted by Lithoprobe and Polaris provide a wealth of seismic, magnetotelluric, and heat flow data. In this paper we start by using and comparing the multitaper and wavelet methods for estimating Te and Te anisotropy. We 37 Chapter 3. Deformation of continents at convergent margins then compare the results with a range of geophysical observables from the same region. The geophysical correlations allow us to formulate a concep- tual model of deformation at convergent margins that can be tested in future studies of lithospheric dynamics. 3.2 Te estimation 3.2.1 Coherence method Te is estimated by minimizing the misfit (L2 norm) between the observed spectral coherence between topography and Bouguer gravity anomalies and the coherence predicted for the flexure of a uniform elastic plate bending under the weights of topography at the surface and internal loads related to density variations [Watts, 2001]. The plate response is modelled either as isotropic or anisotropic, and the coherence is inverted for a single parame- ter, Te, or the three parameters of an assumed orthotropic elastic plate (i.e. having different rigidities in two perpendicular directions), Tmin, Tmax, and φ, the direction of weakest rigidity [Swain and Kirby , 2003]. The coher- ence is an averaging property which measures the consistency of the phase dependence between two signals or fields at distinct wavelengths. Applied to topography and gravity data sets, it provides a measure of the degree of flexural compensation of tectonic loads. An important physical parameter for understanding this property is the characteristic flexural wavelength, λF ∼ ( ET 3e ∆ρg )1/4 , (3.1) where ∆ρ is the density contrast at the depth of compensation, g is the gravitational acceleration, and ”∼” means ”scales as”. At short wave- lengths (λ ≪ λF ) the loads are fully supported by elastic stresses and the Bouguer anomaly is not coherent with the topography. In contrast, at long wavelengths (λ ≫ λF ) loads are compensated isostatically and the Bouguer anomaly is coherent with topography. Typical coherence curves vary smoothly from 0 to 1, reflecting the transition from flexurally sup- 38 Chapter 3. Deformation of continents at convergent margins ported to fully compensated loads. In particular, equation 3.1 indicates that the wavelength of transition from compensated to uncompensated loads increases for increasing plate rigidity. On the continents the elastic thickness estimated using the coherence method ranges between very low (Te < 10 km) and very high (Te > 100 km) values, reflecting the wide variety of tectonic settings and the thermal evo- lution of the lithosphere [Pérez-Gussinyé et al., 2004]. Discrepancies occur, however, when the results from different spectral estimation and load de- convolution techniques applied to the same datasets are compared [Pérez- Gussinyé et al., 2004; Audet and Mareschal , 2007]. This problem arises because Te is a wavelength dependent parameter (equation 3.1), and the accuracy of its estimation is controlled by the spatio-spectral resolution of the analysis technique. Window-based methods have good resolution in frequency, but lack spatial resolution, depending on window width. Spatio- spectral (or space-scale) methods, on the other hand, offer a compromise on the resolution in both reciprocal domains, and are thus able to characterize spatial discontinuities. It is thus important to use and compare different spectral methods whenever possible. In contrast to most previous analyses of Te, generally relying on a single spectral technique [e.g. Lowry and Smith, 1995; Simons et al., 2000; Kirby and Swain, 2006; Audet and Mareschal , 2007], we maximize both spatial and wavelength resolutions by using and comparing two spectral techniques for the calculation of the coherence and the estimation of Te: the multitaper [Simons et al., 2000; Pérez-Gussinyé et al., 2004] and wavelet transform [Kirby and Swain, 2006; Audet and Mareschal , 2007] methods. Applied to western Canada, we show that both techniques recover the spatial pattern of Te and the Te anisotropy, but the wavelet method provides both better spatial resolution and higher absolute estimates of Te. In the inversion, the predicted coherence is calculated by deconvolving the topography and Bouguer anomaly into surface and internal components of the initial loading structure [Forsyth, 1985]. See Pérez-Gussinyé et al. [2004] for the details of this calculation. This step requires information on the crustal structure (Moho depth and density contrasts) that we extract 39 Chapter 3. Deformation of continents at convergent margins from the LITH5.0 crustal model [Perry et al., 2002]. We take the subsurface load to be the density contrast at the Moho [Watts, 2001]. The flexural rigidity is converted to Te using a Poisson ratio of 0.25 and Young modulus of 100 GPa. Resolution tests and statistical analysis of estimated Te values from synthetic modelling show that the accuracy of recovered Te strongly depends on the size of the window compared to the flexural wavelength [Audet and Mareschal , 2007], which is expected from equation (3.1). A reasonable estimate of the accuracy of Te estimated from spectral methods is ±15% of the recovered value. Moreover, as Te is sensi- tive to a range of parameters [Burov and Diament , 1995; Lowry and Smith, 1995], additional uncertainties in the choice of material parameters and mod- elling approaches can be described as log-normal and of the order 30% - 50% [Lowry and Smith, 1995]. Since only major trends are of interest, the vari- ance of Te estimates will not affect the interpretation of the results. 3.2.2 Multitaper analysis The multitaper analysis is performed using a set of orthogonal basis func- tions that minimize the variational problem of broadband bias reduction as data tapers [Thomson, 1982], such that multiple estimates of the same spectrum are approximately uncorrelated and can be averaged to reduce the estimation variance. The set of functions is characterized by two num- bers: the time-frequency bandwidth product NW that controls the wave- length resolution and spectral leakage, and the number of significant tapers k that governs the estimation variance. The multitaper spectral estimates are calculated using Discrete Prolate Spheroidal Slepian (DPSS) tapers with NW = 3 and k = 3 throughout this study. See Pérez-Gussinyé et al. [2004] for the description of the general application of the multitaper method in Te studies. In the multitaper method, the load deconvolution is done in a given rect- angular window of the study area and the estimated Te value is assigned to the area covered by the window [Pérez-Gussinyé et al., 2004]. Spatial reso- lution is achieved by moving the windows over the entire data set with some 40 Chapter 3. Deformation of continents at convergent margins overlap. The size d of the data windows is critical in the estimation of Te be- cause small windows provide high spatial resolution but cannot resolve long flexural wavelengths (c.f. equation (3.1), d < λF ), while large windows (i.e. d≫ λF ) provide an average of the true distribution of Te within the window, thus degrading the spatial resolution. Although a popular method [Pérez- Gussinyé et al., 2004], two important limitations of the multitaper method for Te estimation emerge because of the trade-off between window size and frequency resolution and cause: (1) limited spatial and wavelength resolu- tion, and (2) lower estimates of Te. To counter these effects, we map out re- sults using a combination of window sizes of 300 × 300 km2, 500 × 500 km2, and 800 × 800 km2 (shown at lower left corner of Figure 3a), with an over- lap of 70% between adjacent windows. When the misfit does not converge, we assign a value of 60 km to the window, which is the largest resolvable Te value using the 800 × 800 km2 window. Results are extrapolated to the edges of the study area and the maps are averaged. The spatial recovery of Te is limited by the size of the smallest window and is restricted to a 1484 × 1484 km2 square region in the center of the study area. Finally, the resulting map is filtered using a Gaussian function of width 140 km to produce a continuous image of Te that is more easily comparable to the re- sults obtained with the wavelet transform method. Results do not change for Gaussian functions with different window widths. The multitaper method also allows for the determination of directional variations of the coherence because the spectra are effectively calculated and averaged at each wavevector (kx, ky) in the 2D spectral plane. The results of the inversion for the anisotropic Te cannot be averaged spatially as in the isotropic analysis (above), because the weak directions can vary between estimates obtained with different window sizes. However the results can be qualitatively compared with those obtained with the wavelet method (below) to assess the general trends seen in the map of Te anisotropy. 41 Chapter 3. Deformation of continents at convergent margins 3.2.3 Wavelet transform In the wavelet analysis there is no need for fixed spatial window because the wavelet transform, defined as the convolution of a signal with a scaled and rotated kernel called a wavelet, is a space-scale operation that calcu- lates the spectra at each grid point. The convolution is performed over the entire grid of the study area at different scales, which are related to a Fourier wavenumber. In this study we use the directional 2D Morlet wavelet because it is well suited for the study of Te and the anisotropy [Kirby and Swain, 2006; Audet and Mareschal , 2007]. The Fourier load deconvolution of Forsyth [1985] is performed at each grid point by making the assumption of spatial decoupling between adjacent local wavelet spectra [Stark et al., 2003; Kirby and Swain, 2004; Audet and Mareschal , 2007]. In other words, we assume that the wavelet coefficients at one grid point are statistically independent from neighboring coefficients, and apply a Fourier-based de- convolution to the wavelet spectra. In reality, however, spatial and spectral (space-conjugate) resolutions trade-off in the wavelet analysis as implied by the uncertainty principle in information theory, which states that physical and spectral localizations cannot be measured simultaneously with arbitrar- ily high precision. Short-scale wavelets are well localized in space but broad in the spectral domain. Large-scale wavelets are well localized in the spectral domain but broad in space. This results in a smoothed version of Te where it is large, but allows for larger estimates of Te than the multitaper method due to the presence of longer wavelengths contributing to the spectra. In the anisotropic analysis, the coherence is calculated at distinct az- imuths using a directional Morlet wavelet [Kirby and Swain, 2006; Audet and Mareschal , 2007] on a coarse grid for a better visualization of the varia- tions of the weak axes. See Audet and Mareschal [2007] or Kirby and Swain [2006] for the description of the application of the wavelet transform method in Te studies. 42 Chapter 3. Deformation of continents at convergent margins 3.3 Results 3.3.1 Isotropic Te Figure 3.1 shows example results of the 1D (azimuthally averaged) multita- per coherence for the three window sizes at different locations (N. Cordillera, S-E. Cordillera, and Craton). Misfit curves show a clear minimum for small values of Te even for the smallest window used. As the size of the window increases, more wavelengh information from neighboring regions is retained in the coherence and the estimate can vary by as much as a factor of 5. In the Craton the maximum of the coherence curves remains small even for the largest window size, and the misfits do not converge to a stable Te value, indicating that Te is large in this region. In comparison, the wavelet coher- ence shows the characteristic roll-over from uncompensated (0 coherence) to fully compensated (coherence of 1) loads as well as a clear minimum misfit in all cases (Figure 3.2). Multitaper and wavelet maps of Te are shown in Figure 3.3. Regardless of the technical differences between the two methods, qualitatively, the results agree encouragingly well and correlate with major physiographic provinces in western Canada. Te varies from <30 km in the Phanerozoic Cordillera to >60 km in the Paleoproterozoic Craton over a distance of a few hundred kilometers. This trend is in qualitative agreement with a previous study of Te in western Canada that used the maximum entropy spectral estimator [Flück et al., 2003], although the maximum estimates of Te in the Craton are generally a factor 2 to 3 smaller. Despite a good overall qualitative agreement between the Te maps, there is a pronounced difference in the Craton where multitaper estimates are lower than 30 km in some instances. The discrepancies stem in part from the ad hoc methodology adopted in section 2.2 for the mapping of the win- dowed multitaper estimates, and from the difference in spatial and spectral resolutions between the two methods. In the Craton Te is presumably large, such that small windows would not resolve the transition wavelength from low to high coherence. However, high coherence (>0.3, c.f. Figure 3.1) at short wavelengths (the ”spurious” coherence mentioned in Swain and Kirby 43 Chapter 3. Deformation of continents at convergent margins N. Cordillera S-E. Cordillera Craton Te = 4 km Te = 6 km Te = 10 km Te = 4 km Te = 8 km Te = 14 km Te > 60 km Te > 60 km Te > 60 km 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e a 0.0 0.5 1.0 1.5 2.0 M is fit d 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e b 0.0 0.5 1.0 1.5 2.0 M is fit e 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e 102 101 Wavenumber (1/km) c 0.0 0.5 1.0 1.5 2.0 2.5 3.0 M is fit 10 20 30 40 50 60 70 80 Elastic thickness (km) f 300 x 300 km2 800 x 800 km2 500 x 500 km2 Figure 3.1: Examples of Te estimated with the multitaper method for the three window sizes: a-c) Observed (solid lines) and predicted (dashed lines) coherence curves for three locations in western Canada ((a) North Cordillera, (b) South-East Cordillera, and (c) Craton) ; d-f) L2 norm (root mean square) error between observed and predicted coherence curves for a range of Te values and for each window used. Sizes of windows are 300 × 300 km2 (red curves), 500 × 500 km2 (blue curves), and 800 × 800 km2 (green curves). 44 Chapter 3. Deformation of continents at convergent margins N. Cordillera S-E. Cordillera Craton 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e a 0.0 0.1 0.2 0.3 0.4 0.5 M is fit d Te = 24 km 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e b 0.0 0.1 0.2 0.3 0.4 0.5 M is fit e Te = 16 km 0.0 0.2 0.4 0.6 0.8 1.0 Co he re nc e 102 101 Wavenumber (1/km) c 0.0 0.1 0.2 0.3 0.4 0.5 M is fit 20 40 60 80 100 120 140 Elastic thickness (km) f Te = 51 km Figure 3.2: Examples of Te estimated with the wavelet transform method for the three locations indicated on the plots. The coherence curves show the characteristic roll-over from uncompensated (0 coherence) to fully com- pensated (coherence of 1) loads as well as a clear minimum misfit in all cases 45 Chapter 3. Deformation of continents at convergent margins 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ a 0 10 20 30 40 50 60 70 80 90 Elastic thickness (km) 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ b Cordillera Interior Plains Craton Figure 3.3: Te structure of western Canada calculated from the coherence between Bouguer gravity and topography using the multitaper (a) and the wavelet transform (b) methods. Data are projected onto a Transverse Mer- cator grid. Although the values differ in (a) and (b), the qualitative Te structure is approximately equivalent, in particular the location of the Te gradient. The Te results from the multitaper method are lower than those obtained with the wavelet method due to the use of moving windows of dif- ferent sizes (shown by the boxes at the lower left corner in (a)) which limits the spatial and wavelength resolutions [Simons et al., 2000]. Thus, we con- sider the wavelet results to be more reliable. Major physiographic provinces are bounded by solid white lines. The results show a well defined transition from low values in the Cordillera, intermediate values in the buried Interior Plains, and high values in the exposed Craton. Inset shows the location of the study area with respect to North America. 46 Chapter 3. Deformation of continents at convergent margins [2006] and Audet and Mareschal [2007]) not associated with the flexure of an elastic plate can be mapped as Te and thus lower the Te average over the region. The extent to which this short-wavelength anomaly affects the mapping of Te is thus controlled by the spatial resolution: in the multitaper estimates, it will dominate the Te mapping until a large enough window is able to resolve the flexural wavelength (c.f. equation 3.1, d ∼ λF ). This situation is amplified by the fact that our study area is small, and thus only a small number of the largest windows is used. In the wavelet method the deconvolution is performed at each grid point and thus incorporates less information from neighboring grid points. Because this high coherence is observed at short wavelengths, the mapping of the wavelet estimates of Te is less affected by the spurious coherence. Although a more careful picking of Te estimates according to the Te sen- sitivity for each window size [see for e.g. Flück et al., 2003; Pérez-Gussinyé et al., 2004] could perhaps erase this low-Te anomaly, we prefer to avoid mak- ing any a priori assumption about the underlying Te structure because this introduces an additional observational bias. Alternatively, Pérez-Gussinyé et al. [2007] calculated multitaper Te estimates using three window sizes and produced three maps of Te in South America that compared favorably with the wavelet method in the same region [Tassara et al., 2007]. Another ap- proach is to weight the misfit by the inverse of wavenumber, as was done in Kirby and Swain [2006]; Swain and Kirby [2006] and Tassara et al. [2007]. This filtering procedure downweights the contribution of the high coherence at small wavenumbers, effectively erasing the low-Te anomaly, at the cost of spatial resolution [Kirby and Swain, 2006; Pérez-Gussinyé et al., 2007]. We note that the wavelet map of Te also shows an intriguing low-Te region at the north-east corner of the study area, with a lesser spatial extent. Such low-Te estimates due to short-wavelength coherence is also found elsewhere in the Canadian Shield [Audet and Mareschal , 2004, 2007]. The origin of those features is beyond the scope of this study. We suggest that a more thorough comparison between the wavelet and multitaper methods applied to Te is needed in order to isolate the contribu- tion of the methods themselves to the differences observed. 47 Chapter 3. Deformation of continents at convergent margins 1a 1b 1c 1d 1e 1f 2a 2b 2c 2d 2e 2f 3a 3b 3c 3d 3e 3f 0.0 0.2 0.4 0.6 0.8 1.0 Coherence 0 pipi/2−pi −pi/2 0 pipi/2−pi −pi/2 0 pipi/2−pi −pi/2 0 pi pi/2 −pi −pi/2 0 pi pi/2 −pi −pi/2 0 pi pi/2 −pi −pi/2 0 pi pi/2 −pi −pi/2 0 pi pi/2 −pi −pi/2 0 pi pi/2 −pi −pi/2 kx kx ky ky ky ky ky ky N. Cordillera S-E. Cordillera Craton Figure 3.4: Anisotropic inversion of Te using the multitaper method. Ob- served (a-c) and theoretical (d-f) 2D coherences using different windows: (1a-f) is 300 × 300 km2; (2a-f) is 500 × 500 km2; and (3a-f) is 800 × 800 km2. Wavenumbers are normalized by each window size. The observed 2D coher- ences correspond to the three locations described in Figure S1: (1-3a) N. Cordillera, (1-3b) S-E. Cordillera, and (1-3c) Craton. The major axis of the high-coherence ellipse points in the direction of lower Te, or Tmin. 48 Chapter 3. Deformation of continents at convergent margins 3.3.2 Anisotropic Te Figure 3.4 shows example results for the anisotropic coherence calculated with the multitaper method at the same three locations as those in Fig- ure 3.1. Also shown are best-fitting 2D coherences using analytical solu- tions of the orthotropic elastic plate [Kirby and Swain, 2006; Audet and Mareschal , 2007]. In a 2D wavenumber plot (kx, ky) of coherence, anisotropic Te typically produces an elliptical shape of high coherence, whereas isotropic Te appears as a circular region of high coherence. Figure 3.4 indicates that the anisotropy is well-developed in the N. Cordillera as the results for the three window sizes are consistent. In the S-E. Cordillera the coherence is quasi isotropic except for the estimate obtained with the largest window. The Craton estimates show an absence of coherence at the longest wave- lengths as in the 1D case, indicating that still larger windows are necessary to resolve Te and the anisotropy in this region. The wavelet results (Figure 3.5) are presented in 2D as wavelength- an- gle plots. The signature of anisotropic Te is then characterized by a 180 ◦ periodic variation of high coherence, whereas the signature of isotropic Te is a flat line. The wavelet results agree with the multitaper results and show that the anisotropy is well developed in the N. Cordillera and in the S-E. Cordillera, but almost absent in the Craton. Maps of Te anisotropy us- ing the multitaper method show considerable scatter with the major trends described above, and we choose not to interpret them. The resulting map of the anisotropic inversion of Te using the wavelet method is shown in Figure 3.6. The magnitude of Te anisotropy is given by the ratio (Tmax−Tmin)/Tmax. Te anisotropy is maximized in the Cordillera and across the Interior Plains where Te < 30 km, closely mimicking the isotropic Te pattern. The mechanical weak (i.e. low rigidity) direction is dominantly SW-NE in this region, approximately perpendicular to the main belts of the Cordillera and parallel to the direction of Cascadia subduction. East of Te = 30 km contour, the Te anisotropy progressively decreases in magnitude. In addition, the trend of the weak direction, φ, rotates to a SE-NW direction in the Craton and in the northernmost part of the Interior 49 Chapter 3. Deformation of continents at convergent margins a 20 31 47 71 109 167 255 389 595 908 1387 1850 W av e le ng th (k m) d (φ,Tmin,Tmax) = (205o, 5km, 30km) b 20 31 47 71 109 167 255 389 595 908 1387 1850 W av e le ng th (k m) e (φ,Tmin,Tmax) = (215o, 10km, 22km) 60 90 120 150 180 210 Azimuth (deg) c 20 31 47 71 109 167 255 389 595 908 1387 1850 W av e le ng th (k m) 60 90 120 150 180 210 Azimuth (deg) f (φ,Tmin,Tmax) = (125o, 50km, 60km) N. Cordillera S-E. Cordillera Craton Figure 3.5: Anisotropic inversion of Te using the wavelet transform method. The coherence is plotted as function of wavelength and azimuth. Observed (a-c) and theoretical (d-f) 2D coherences for the three locations: (a) N. Cordillera, (b) S-E. Cordillera, and (c) Craton. Legend in 2D theoretical plots indicates best-fitting parameters of the orthotropic elastic plate. The anisotropy is well developed in the N. Cordillera and in the S-E. Cordillera, but almost absent in the Craton. The scale bar is the same as in Figure 3.4 50 Chapter 3. Deformation of continents at convergent margins 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ 0.0 0.2 0.4 0.6 0.8 1.0 (Tmax −Tmin)/Tmax a 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ 0 10 20 30 40 50 60 70 80 90 Tmin (km) b 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ 0 10 20 30 40 50 60 70 80 90 Tmax (km) c Figure 3.6: Anisotropic Te structure of western Canada from the wavelet transform method. The 2D coherence is inverted using the thin orthotropic plate model, i.e. having different rigidities in two perpendicular directions, for three parameters: the minimum (b) and maximum (c) values of Te, and φthe orientation of the weak direction calculated anti-clockwise from the x-axis (a). The weak direction of mechanical anisotropy in (a) is plotted as thin solid black bars superposed on a color scale map of the anisotropic magnitude. The patterns of Tmin and Tmax closely mimick the isotropic Te structure of Figure 3.3. In particular, the maximum gradient occurs at the same location, and accompanies a change in both the magnitude and orientation of the weak direction. 51 Chapter 3. Deformation of continents at convergent margins Plains. We note also that the maps of Tmin and Tmax follow closely the isotropic Te structure of Figure 3.3. In particular, the maximum gradient occurs at the same location. 3.3.3 Geophysical correlations In Figure 3.7a we plot the location of highest Te gradient using the wavelet results, along with the location of earthquake epicenters for all events with magnitude larger than M = 3 that occured since 1979. On the map is also shown a 150 × 1800 km2 corridor trending SW-NE within which Te, Te gradient (|∇Te|), and crustal thickness (Tc) were extracted and aver- aged, and plotted in Figure 3.7b. Figure 3.7c shows Te with the percentage of seismicity and heat flow values along the same profile. The highest Te gradient is found in the Interior Plains, 100 to 300 km east of the largest concentration of seismicity, and also where Te < 30 km. The position and magnitude of the Te gradient also broadly correlate with a similar (i.e. fac- tor of about 2) eastward decrease in vertical heat flux [Hyndman and Lewis, 1999], consistent with inferences from xenolith thermobarometry and seis- mic data [Currie and Hyndman, 2006]. The Te gradient corresponds also to a transition in crustal thickness from ∼30 km in the Cordillera to >40 km in the Craton [Perry et al., 2002], and with a gradient of low to high S-wave seismic velocity anomalies imaged by seismic tomography [Frederiksen et al., 2001]. In Figure 3.8a we plot the Te weak directions along with seismic and magnetotelluric measurements of crustal and upper mantle anisotropy and indicators for crustal stresses. Figure 3.8b shows the magnitude of the Te anisotropy with averages of Tmin and Tmax along the same profile as Fig- ure 3.7. In the southern Cordillera and across the Interior Plains the fast axis of seismic anisotropy, as measured from SKS splitting data, show coher- ent fast polarizations in the SW-NE direction [Shragge et al., 2002; Currie et al., 2004]. The geoelectric strikes align in the same direction as the SKS data [Boerner et al., 1999; Ledo and Jones, 2001]. Both anisotropy indi- cators align with the Te weak directions and with the maximum horizontal 52 Chapter 3. Deformation of continents at convergent margins 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ a 0.0 0.5 1.0 1.5 T e gr ad ie nt 0 10 20 30 40 50 60 70 80 90 D ep th [k m] Te Te = 30 km Tc |∇Te| b 0 2 4 6 8 10 % S ei sm ici ty 0 200 400 600 800 1000120014001600 Distance along profile [km] Te Te = 30 km 0 30 60 90 120 H ea t f lo w [m W /m 2 ] c Figure 3.7: Correlations between the Te structure and geophysical data. In (a) the seismicity (with magnitude > 3) of western Canada, represented by black circles at earthquake epicenters recorded between 1979 and 2006, is plotted against the location of the maximum gradient of Te shown in orange. Thick grey lines are the province boundaries. Dashed dark red line is the 30 km contour which roughly follows the maximum gradient limit. (b) shows the average of Te, Te gradient, and Tc calculated within the corridor delimited by the thick gray dashed lines in (a), and projected along the profile shown by the thick gray line. In (c) the Te is plotted against the seismicity (gray bins) and heat flow values (diamonds) extracted within the same corridor and projected along the profile. The west-east decrease in seismicity correlates with the location of the Te gradient, where Te ∼ Tc ∼ 30 km, and where heat flow decreases by a factor of 2. 53 Chapter 3. Deformation of continents at convergent margins compressive stress directions. This striking uniform pattern of anisotropy persists over a wide area to about 57◦ latitude, and correlates with a change in Te from Te < Tc to Te > Tc. This also corresponds to the boundary defined by the largest Te gradient, where Te ∼30 km, and where the Te anisotropy changes in both amplitude and direction. In order to quantitatively compare the anisotropy datasets, we resampled the seismic, magnetotelluric and stress data onto the coarse Te anisotropy grid. Figure 3.9 shows histograms of the difference in angle between the Te weak directions and the seismic, magnetotelluric, and stress directions, along with a histogram of percent alignment to within 15◦ and 22.5◦ [Simons and van der Hilst , 2003]. Stress and magnetotelluric data show well-defined Gaussian curves centered on ∼10◦ to 20◦ difference with Te anisotropy (Fig- ure 3.9b,c), and the alignments are well above the 15◦ and 22.5◦ signif- icance levels (Figure 3.9d). The angle differences between Te weak axes and seismic fast axes show a quasi-bimodal distribution centered on ∼10◦ and ∼75◦ (Figure 3.9a). When binned into three distinct regional clusters (Cascadia backarc, Cordilleran Front, and Craton), the alignment between Te weak axes and the seismic fast axes in the Cordilleran portion emerges (Figure 3.9a,d, light blue). The Cascadia and Craton regions show angle differences larger than 45◦ with a marked proportion of data towards or- thogonality. However more seismic data are required to assess the general trends observed. 3.4 Discussion 3.4.1 Te, Tc, and heat flow The relationship between Te and Tc requires some discussion because it provides insight into the factors that influence Te most strongly [Burov and Diament , 1995; Lowry and Smith, 1995]. That Te < Tc in the Cordillera suggests that the average elastic properties of this region are governed, in part, by the composition and constitution of crustal rocks. In contrast, the elastic properties of the Craton, where Te > Tc, must involve a component 54 Chapter 3. Deformation of continents at convergent margins 230˚ 240˚ 240˚ 250˚ 250˚ 50˚ 60˚ 60˚ a Te weak axis SKS splitting MT aniso Max. comp. stress 0.0 0.5 1.0 T e a n is ot ro py 0 200 400 600 800 1000 1200 1400 1600 Distance along profile [km] 0 10 20 30 40 50 60 70 80 90 D ep th [k m] Tmin Tmax Te = 30 km Tc (Tmax−Tmin)/Tmax b Figure 3.8: Anisotropy correlations in western Canada. In (a) the mechan- ical weak axes (black bars) are plotted along with seismic (SKS splitting - blue bars [Shragge et al., 2002; Eaton et al., 2004; Currie et al., 2004], and electrical (geoelectric strikes - orange bars [Boerner et al., 1999]) anisotropy data. Green and red bars are maximum compressive horizontal stress di- rections inferred from shear-wave splitting above local earthquakes in the Cascadia forearc [Currie et al., 2001] and borehole breakouts, respectively. Dashed dark red line is the Te = 30 km contour. (b) shows the average of Tmin, Tmax, and the magnitude of Te anisotropy calculated within the corridor and projected along the profile as in Figure 3.7b,c. The boundary defined by Te ∼ 30 km correlates with a decrease in Te anisotropy. In the southern Cordillera and across the Interior Plains all indicators align in the direction of plate motion and subduction of Juan de Fuca beneath North america, suggesting a stress control on the deformation of western Canada. 55 Chapter 3. Deformation of continents at convergent margins 0 10 20 30 40 50 Fr eq ue nc y (% ) −90 −60 −30 0 30 60 90 Angular difference <∆θ> = 48.5 ± 29.0 <∆θ> = 7.8 ± 6.1 <∆θ> = 57.7 ± 23.3 a) SKS 0 8 16 24 32 40 Fr eq ue nc y (% ) −90 −60 −30 0 30 60 90 Angular difference <∆θ> = 18.9 ± 22.3 <∆θ> = 31.5 ± 24.7 b) Stress 0 4 8 12 16 20 Fr eq ue nc y (% ) −90 −60 −30 0 30 60 90 Angular difference <∆θ> = 11.1 ± 22.6 c) MT 0 20 40 60 80 100 % A lig nm en t d SK S Co rd ille ra SK S Cr at on SK S Ca sc ad ia M T St re ss Ca sc ad ia Figure 3.9: Statistics of anisotropy to test whether the observed alignment between Te weak directions and other anisotropy indicators is significant. In (a), (b), and (c) the histograms show the angular difference between the Te weak axes and seismic fast axes, maximum horizontal compressive stress directions, and geoelectric strikes respectively. The seismic anisotropy in (a) is further divided into three regions corresponding to the Cordillera (light blue), the Craton (dark blue), and the Cascadia backarc (hatched dark blue). Shown in (d) are the percentage of alignment to within 15◦ (colored bins), and to within 22.5◦ (light gray). Concinuous and dashed black lines indicate the levels at which randomly distributed alignment should fall within 15◦ and 22.5◦. High percentage of alignment suggests a causal link between the observables. 56 Chapter 3. Deformation of continents at convergent margins of the lithospheric mantle, although the nature and extent of the mechanical coupling of the crust to the underlying lithosphere is unknown. In principle the vertical structure of the lithosphere and its influence on plate rigidity is determined, in part, locally and regionally by geology and may be complex. However, the eastward increase in Te is significantly greater than that which is expected if plate rigidity is related to changes in Tc alone (i.e. D is not proportional to T 3c ). Indeed, the roughly factor of 2 variation in Te is better explained by the similar variations in heat flux that will scale approximately as 1/Te if Te is proportional to the lithospheric thickness [Maggi et al., 2000]. This result suggests that although a geologically complex region, Te is predominantly governed by temperature in western Canada, as speculated in previous studies [Hyndman and Lewis, 1999; Flück et al., 2003]. 3.4.2 Anisotropy and stress: A primer Seismic anisotropy refers to variations of seismic wavespeeds along differ- ent directions of propagation or polarization [Silver , 1996; Savage, 1999]. The most popular method for investigating seismic anisotropy is the split- ting (or birefringence) of shear-waves travelling through the mantle. In an anisotropic mantle with a horizontal axis of symmetry, the two perpendicu- lar shear-wave polarizations will split into a fast and a slow component. The time difference between the two phases provides an estimate of the magni- tude of the anisotropy, whereas the fast polarization azimuth measures the direction of the symmetry axis. Time delays and azimuths give an indication of elastic anisotropy within the lithosphere, and can thus image deformation of the rocks in the Earth’s mantle. Under the continents, the seismic anisotropy is attributed either to frozen fabric in the lithospheric mantle [Silver , 1996], or to present-day sublitho- spheric mantle convection [Vinnik et al., 1992]. In western Canada, almost all the seismic anisotropy data are measured by the splitting of SKS phases and are interpreted to result from shear-induced lattice-preferred orienta- tion (LPO) of olivine in the upper mantle due to viscous flow [Bank et al., 2000; Shragge et al., 2002; Currie et al., 2004]. Despite the excellent lateral 57 Chapter 3. Deformation of continents at convergent margins resolution afforded by this method, these measurements are hampered by a poor depth resolution and/or the presence of complex anisotropic layering [Park and Levin, 2002], hindering a robust estimation of the source of the anisotropy. Electrical anisotropy is seen by the azimuthal dependence of the phase difference between the two orthogonal horizontal components (transverse electric - TE - and transverse magnetic - TM) of the time- dependent electro- magnetic (EM) fields of magneto-telluric (MT) measurements [Boerner et al., 1999; Eaton et al., 2004; Jones, 2006]. The magnitude of the electrical anisotropy is generally of a few orders of magnitude, rather than a few per- cents as in shear-wave anisotropy, and hence rules out the possibility of an intrinsic dry crystal or grain electrical anisotropy. Electrical anisotropy may be generated by macroscopic, graphitized (resistive) shear zones within the lithosphere that mark past tectonic events [Mareschal et al., 1995], by grain- scale anisotropy of strain-induced LPO of wet olivine crystals oriented by present-day plate motion or mantle flow [Simpson, 2001], or by enhanced fluid flow in the crust and upper mantle due for example to water being incorporated into the subduction zone [Ledo and Jones, 2001; Soyer and Unsworth, 2006] or originating from a meteoric source. MT strike directions at a particular period represent an integrative effect of the electrical structure at a broad depth range related to the diffusion of the EM fields into a conductor. However the penetration depth of the TE field can differ from that of the TM field by tens of kilometers due to strong structural heterogeneity [Jones, 2006]. Electrical anisotropy data shown on Figure 3.8 show the strike directions at a period of 320 s, which is broadly representative of mid-crustal depths [Boerner et al., 1999]. A reliable indi- cator of the electrical anisotropy magnitude is the phase separation between the impedance estimates for the two transverse modes that is consistent on a regional scale. In the Interior Plains the phase difference is maximum south of 57◦, and diminishes to the north and southeast [Boerner et al., 1999]. Electrical anisotropy in this region requires a hydrothermal alteration event most likely related to the tectonic assembly of Laurentia circa 1850 Ma [Boerner et al., 1999]. Note that the high correlation between present-day 58 Chapter 3. Deformation of continents at convergent margins stress and MT anisotropy rather suggests a causal link and we provide an alternative explanation for the MT anisotropy in the next section. Seismic and electrical anisotropies mostly sample upper mantle/mid- crustal materials, but although they do not measure the same physical properties, they may be measuring the response of the same causative effect. The anisotropy in Te provides an integrated measure of the rheology of the lithosphere at crustal and upper mantle depths and thus it fills the gap be- tween these anisotropy measurements and surface geology features [Lowry and Smith, 1995; Simons and van der Hilst , 2003]. A number of mechanisms can induce mechanical anisotropy of the litho- sphere. The azimuthal dependence of the coherence between gravity anoma- lies and topography indicates that flexural compensation is facilitated in the direction perpendicular to the direction that has accumulated the most deformation per unit of topographic loading, which is accomodated by fold- ing, faulting, or buckling of the crust (i.e. the weak axis is perpendicular to strike) [Simons and van der Hilst , 2003]. If the preserved gravity struc- ture reflects the current state of stress of the lithosphere due for example to plate convergence, then the weak direction aligns in the direction of maxi- mum compressive stress [Lowry and Smith, 1995], and that of shear strain induced seismic anisotropy if flow is parallel to convergence. High ampli- tude Te anisotropy thus reflects the preserved strain field that is likely to be related to crustal stresses and associated brittle processes (e.g. faulting, seismicity) at active convergent margins. 3.4.3 Conceptual model of deformation The diverse geophysical correlations of Figures 3.7 and 3.8, and the dis- cussion in section 4.2, suggest a simple conceptual model for the current state of stress and the associated brittle deformation in western Canada (Figure 3.10). The alignment of the seismic anisotropy data with crustal stress indicators suggests that elastic stresses, expressed as both flexurally- supported topography, and brittle failure, indicated by earthquakes, arise in response to driving plate boundary forces and retarding mantle tractions 59 Chapter 3. Deformation of continents at convergent margins related to Cascadia subduction and associated mantle flow. The sharp and nearly factor of 10 increase in the thickness of the elastic lithosphere from the Cordillera to the Craton is interpreted to cause bending related to a compressional component of the stress field to be concentrated a couple of hundred kilometers west of the Te gradient. The location of this stress concentration is consistent with expectations from simple calculations of the flexural characteristics of beams and plates of variable rigidity. Elastic bending stresses scale as ( E/(1 − ν2)) (Te/R), where ν and R are Poisson’s ratio and the radius of curvature of the deformation, and will exceed the ∼30 MPa yield strength of the crust for even very small vertical deflections of the cordilleran lithosphere (i.e. for R < O(104) km, say). That seismicity and mechanical weakening occurs where bending stresses are maximized is thus not surprising. If the bending is simple (approximately cylindrical, Figure 3.10) failure is expected to occur in tension (leading to normal faults) and compression (leading to reverse faults) in the upper and lower halves of the elastic lithosphere, respectively. Results from geologic mapping [Carr et al., 1988] and regional GPS studies are consistent with this picture and, moreover, indicate that significant upper crustal extension and normal faulting has occurred in this region over an extended period of time [Carr et al., 1988]. It is worth noting that our conceptual elastic model predicts the location of failure, regardless of the amount of finite-amplitude deformation produced thereafter, which involves no elasticity. The enhanced west-east electrical conductivity is sufficiently large to require a fluid phase [Boerner et al., 1999; Ledo and Jones, 2001; Soyer and Unsworth, 2006], where the fluids may come either from the mantle or meteoric sources. If this electrical anisotropy is explained by a conductivity that is enhanced by charge being carried by the flow of crustal fluids it is reasonable to speculate that mechanical failure of the rocks in this region leads also to fracture permeability, with enhanced flow along pre-existing grain boundaries, and that this mechanism dominates over macroscopic fluid flow along faults, which would create a NW-SE trending MT anisotropy. Interestingly, the maximum compressive direction of crustal stress indi- cators in the Cascadia forearc, inferred from the splitting of direct shear- 60 Chapter 3. Deformation of continents at convergent margins Figure 3.10: Cartoon showing a geodynamic interpretation of the elastic thickness results. The driving forces for lithospheric deformation at the convergent margin of western Canada are due to the subduction of the Juan de Fuca plate beneath North America. The elastic plate flexes in response to the loads applied at the plate boundary, and the stresses concentrate where rigidity is lowest, west of the maximum gradient of Te. Brittle failure and seismicity occur in the upper and lower halves of the plate as a result of displacement along normal and reverse faults, respectively. Note that the amplitude of the flexure is significantly exaggerated (to finite-amplitude) to highlight the underlying mechanics of the flexure and failure processes. 61 Chapter 3. Deformation of continents at convergent margins waves generated by local earthquakes [Currie et al., 2001] and in-situ mea- surements, are oriented in a NW-SE direction (Figure 3.7b). These obser- vations are in agreement with the forearc sliver model of [Wang , 1996] who predicts large arc-parallel compressive stresses due to oblique subduction. The anisotropy in Te in this region is in agreement with the predictions and observations of crustal stresses in the Cascadia forearc. The northern Cordillera represents a remarkably tectonically and seis- mically active region where the collision of the Yakutat Block drives the deformation at least 1000 km farther inland [Hyndman et al., 2005]. Our limited coverage of this region hampers a more quantitative comparison of the mechanical controls on the deformation, although the very steep gradi- ent in Te north of 60 ◦ suggests that the same mechanisms operate in this region (see Figure 3.7). More Te data north of 60 ◦ latitude along with con- straints from other geophysical observables are needed to verify the validity of our mechanical model in the northern Cordillera. The conceptual model explains the stress-induced deformation where the lithosphere is weakest and where deformation is most likely controlled by the rheology of the crust, but fails to explain all the features observed on the anisotropy correlation map (Figure 3.8), notably in the Craton. In this re- gion Te is larger than Tc but much lower than the ∼300 km thick lithosphere, suggesting that at least part of the lithospheric mantle contributes to the rigidity. The fact that seismic fast axes align in the direction of absolute plate motion indicates that the anisotropy is consistent with shear-induced LPO of olivine in the upper mantle due to viscous flow [Bank et al., 2000]. Because the very rigid Craton largely escapes deformation and is unaffected by present-day stress levels within the lithosphere, the Te anisotropy therein most likely reflects the fossil strain field from past tectonic events and may not correlate with the deformation associated with the present-day absolute plate motion at the base of the lithosphere [Simons and van der Hilst , 2003]. On the other hand, if the seismic anisotropy were due to frozen fabric within the lithosphere, the seismic fast axes should anti-correlate with the Te weak directions [Simons and van der Hilst , 2003]. The absence of constraints from other geophysical observables renders this hypothesis difficult to assert in a 62 Chapter 3. Deformation of continents at convergent margins testable way. 3.5 Conclusion The estimation of the spatial variations of the elastic thickness is a complex spatio-spectral problem because the elastic response of the plate to load- ing is dominated by features with wavelengths much larger than Te. To address this issue, we use the multitaper and wavelet transform methods to obtain detailed maps of Te and its anisotropy in western Canada. Our results are compared to relevant data sets and suggest that under the cur- rent tectonic stress regime the locations of brittle fracturing of crustal rocks, indicated by seismicity, mechanical anisotropy, and enhanced electrical con- ductivity (i.e. fracture permeability), are governed by the detailed spatial variations in plate rigidity. Whether or not this simple conceptual model is validated by more thorough geodynamic modelling of the deformation of western Canada, our work indicates that an accurate determination of the spatial variations in Te is fundamental to understanding the modern pro- cesses governing lithospheric deformation at active continental margins in general. 3.6 Acknowledgments This paper was much improved by thorough and constructive reviews by three anonymous reviewers. 63 Bibliography Audet, P., and J.-C. Mareschal, Anisotropy of the flexural response of the lithosphere in the Canadian Shield, Geophys. Res. Lett., 31, L20,601, doi:10.1029/2004GL021,080, 2004. Audet, P., and J.-C. Mareschal, Wavelet analysis of the coherence be- tween Bouguer gravity and topography: Application to the elastic thickness anisotropy in the Canadian Shield, Geophys. J. Int., 168, 287–298, 2007. Bank, C.-G., M. G. Bostock, R. Ellis, and J. Cassidy, A reconnaissance teleseismic study of the upper mantle and transition zone beneath the Archean Slave craton in NW Canada, Tectonophysics, 319, 151–166, 2000. Becker, T. W., C. Facenna, R. J. O’Connell, and D. Giardini, The develop- ment of subduction in the upper mantle: insights from experimental and laboratory experiments, J. Geophys. Res., 104, 15,207, 1999. Billen, M. I., M. Gurnis, and M. Simons, Multiscale dynamic models of the Tonga-Kermadec subduction zone, Geophys. J. Int., 153, 359–388, 2003. Boerner, D. E., R. D. Kurtz, J. A. Craven, G. M. Ross, and F. W. Jones, A synthesis of electromagnetic studies in the Lithoprobe Alberta Basement Transect: constraints on Paleoproterozoic indentation tectonics, Can. J. Earth Sci., 37, 1509–1534, 1999. Burov, E. B., and M. Diament, The effective elastic thickness (Te) of con- tinental lithosphere: What does it really mean?, J. Geophys. Res., 100, 3905–3927, 1995. 64 Bibliography Carr, S. D., R. Parrish, and D. L. Parkinson, Eocene extensional tectonics and geochronology of the southern Omineca Belt, British Columbia and Washington, Tectonics, 7, 181–212, 1988. Currie, C. A., and R. D. Hyndman, The thermal structure of subduction zone backarcs, J. Geophys. Res., 111, B08,404, doi:10.1029/2005JB004,024, 2006. Currie, C. A., J. F. Cassidy, and R. D. Hyndman, A regional study of shear wave splitting above the Cascadia subduction zone: Margin-parallel crustal stress, Geophys. Res. Lett., 28, 659–662, 2001. Currie, C. A., J. F. Cassidy, R. D. Hyndman, and M. G. Bostock, Shear wave anisotropy beneath the Cascadia subduction zone and western North American craton, Geophys. J. Int., 157, 341–353, 2004. Eaton, D., A. Jones, and I. Ferguson, Lithospheric anisotropy structure in- ferred from collocated teleseismic and magnetotelluric observations: Great Slave Lake shear zone, northern Canada, Geophys. Res. Lett., 31, L19,614, doi:10.1029/2004GL020,939, 2004. Eaton, D. W., and A. G. Jones, Tectonic fabric of the subcontinental litho- sphere: Evidence from seismic, magnetotelluric and mechanical anisotropy, Phys. Earth Planet. Int., 158, 85–91, 2006. Flück, P., R. D. Hyndman, and C. Lowe, Effective elastic thickness Te of the lithosphere in western Canada, J. Geophys. Res., 108, 2430, doi:10.1029/2002JB002,201, 2003. Forsyth, D. W., Subsurface loading and estimates of the flexural rigidity of continental lithosphere, J. Geophys. Res., 90, 12,623–12,632, 1985. Frederiksen, A. W., M. G. Bostock, and J. F. Cassidy, S-wave velocity structure of the Canadian upper mantle, Phys. Earth Planet. Inter., 124, 175–191, 2001. 65 Bibliography Gaspar-Escribano, J. M., J. D. van Wees, M. ter Voorde, S. Cloetingh, E. Roca, L. Cabrera, J. A. Munoz, P. A. Ziegler, and D. Garcia-Castellanos, Three-dimensional flexural modelling of the Ebro Basin (NE Iberia), Geo- phys. J. Int., 145, 349–367, 2001. Hyndman, R. D., and T. J. Lewis, Geophysical consequences of the Cordillera-Craton thermal transition in southwestern Canada, Tectono- physics, 306, 397–422, 1999. Hyndman, R. D., P. Flueck, S. Mazzotti, T. J. Lewis, J. Ristau, and L. Leonard, Current tectonics of the northern Canadian Cordillera, Can. J. Earth Sci., 42, 1117–1136, 2005. Jones, A. G., Electromagnetic interrogation of the anisotropic Earth: Look- ing into the Earth with polarized spectacles, Phys. Earth Planet. Int., 158, 281–291, 2006. Kirby, J. F., and C. J. Swain, Global and local isostatic coherence from the wavelet transform, Geophys. Res. Lett., 31, doi:10.1029/2004GL021,569, 2004. Kirby, J. F., and C. J. Swain, Mapping the mechanical anisotropy of the lithosphere using a 2-D wavelet coherence, and its application to Australia, Phys. Earth Planet. Int., 158, 122–138, 2006. Kohlstedt, D. L., J. B. Evans, and S. J. Mackwell, Strength of the litho- sphere: Constraints imposed by laboratory experiments, J. Geophys. Res., 100, 17,587–17,602, 1995. Ledo, J., and A. G. Jones, Regional electrical resistivity structure of the southern Canadian Cordillera and its physical interpretation, J. Geo- phys. Res., 106, 30,755–30,769, 2001. Lowry, A. R., and R. B. Smith, Strength and rheology of the western U. S. Cordillera, J. Geophys. Res., 100, 17,947–17,963, 1995. 66 Bibliography Maggi, A., J. A. Jackson, D. McKenzie, and K. Priestley, Earthquake fo- cal depths, effective elastic thickness, and the strength of the continental lithosphere, Geology, 28, 495–498, 2000. Mareschal, M., R. L. Kellett, R. D. Kurtz, J. N. Ludden, S. Ji, and R. C. Bailey, Archaean cratonic roots, mantle shear zones and deep electrical anisotropy, Nature, 375, 134–137, 1995. Park, J., and V. Levin, Seismic anisotropy: Tracing plate dynamics in the mantle, Science, 296, 485–489, 2002. Peacock, S. M., Fluid processes in subduction zones, Science, 248, 329–337, 1990. Pérez-Gussinyé, M., A. Lowry, A. B. Watts, and I. Velicogna, On the recovery of effective elastic thickness using spectral methods: Examples from synthetic data and from the Fennoscandian Shield, J. Geophys. Res., 109, B10,409, doi:10.1029/2003JB002,788, 2004. Pérez-Gussinyé, M., A. R. Lowry, and A. B. Watts, Effective elastic thick- ness of South America and its implications for intracontinental deforma- tion, Geochem. Geophys. Geosys., 5, Q05,009, doi:10.1029/2006GC001,511, 2007. Perry, H. K. C., D. W. S. Eaton, and A. M. Forte, LITH5.0: a revised crustal model for Canada based on Lithoprobe results, Geophys. J. Int., 150, 285–294, 2002. Price, R. A., The southeastern Canadian Cordillera: Thrust faulting, tec- tonic wedging, and delamination of the lithosphere, J. Struct. Geol., 8, 239–254, 1986. Rajesh, R. S., J. Stephen, and D. C. Mishra, Isostatic response and anisotropy of the Eastsern Himalayan-Tibetan Plateau: A reap- praisal using multitaper spectral analysis, Geophys. Res. Lett., 30, 1060, doi:10.1029/2002GL016,104, 2003. 67 Bibliography Savage, M. K., Seismic anisotropy and mantle deformation: What have we learned from shear wave splitting?, Rev. Geophys., 37, 65–106, 1999. Shragge, J., M. G. Bostock, C. G. Bank, and R. M. Ellis, Integrated tele- seismic studies of the southern Alberta upper mantle, Can. J. Earth Sci., 39, 399–411, 2002. Silver, P. G., Seismic anisotropy beneath the continents: Probing the depths of geology, Annu. Rev. Earth Planet. Sci., 24, 385–432, 1996. Simons, F. J., and R. D. van der Hilst, Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere, Earth Planet. Sci. Lett., 211, 271–286, doi:10.1016/S0012–821X(03)00,198– 5, 2003. Simons, F. J., M. T. Zuber, and J. Korenaga, Isostatic response of the Aus- tralian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper spectral analysis, J. Geophys. Res., 105, 19,163–19,184, 2000. Simpson, F., Resistance to mantle flow inferred from the electromagnetic strike of the Australian upper mantle, Nature, 412, 632–635, 2001. Soyer, W., and M. Unsworth, Deep electrical structure of the northern Cascadia (British Columbia, Canada) subduction zone: Implications for the distribution of fluids, Geology, 34, 53–56, doi:10.1130/G21,951.1, 2006. Stark, C. P., J. Stewart, and C. J. Ebinger, Wavelet transform mapping of effective elastic thickness and plate loading: Validation using synthetic data and application to the study of southern African tectonics, J. Geo- phys. Res., 108, 2258, doi:10.1029/2001JB000,609, 2003. Stegman, D. R., J. Freeman, W. P. Schellart, L. Moresi, and D. May, Influence of trench width on subduction hinge retreat rates in 3-D models of slab rollback, Geochem. Geophys. Geosys, 7, Q03,012, 2006. Stern, R. J., Subduction zones, Rev. Geophys., 40, 1012, doi:10.1029/2001RG000,108, 2002. 68 Bibliography Swain, C. J., and J. F. Kirby, The coherence method using a thin anisotropic plate model, Geophys. Res. Lett., 30, 2014, doi:10.1029/2003GL018,350, 2003. Swain, C. J., and J. F. Kirby, An effective elastic thickness map of Australia from wavelet transforms of gravity and topography using Forsyth’s method, Geophys. Res. Lett., 33, L02,314, doi:10.1029/2005GL025,090, 2006. Tassara, A., C. Swain, R. Hackney, and J. Kirby, Elastic thickness structure of South America estimated using wavelets and satellite-derived gravity data, Earth Planet. Sci. Lett., 253, 17–36, 2007. Thomson, D. J., Spectrum estimation and harmonic analysis, Proc. IEEE, 70, 1055–1096, 1982. Vinnik, L. P., L. I. Makeyeva, A. Milev, and A. Y. Usenko, Global pat- terns of azimuthal anisotropy and deformations in the continental mantle, Geophys. J. Int., 111, 433–447, 1992. Wang, K., Simplified analysis of horizontal stresses in a buttressed forearc sliver at an oblique subduction zone, Geophys. Res. Lett., 23, 2021–2024, 1996. Watts, A. B., Isostasy and Flexure of the Lithosphere, Cambridge Univ. Press, Cambridge, UK, 2001. Watts, A. B., and E. B. Burov, Lithospheric strength and its relationship to the elastic and seismogenic layer thickness, Earth Planet. Sci. Letters, 213, 113–131, 2003. 69 Chapter 4 Teleseismic waveform modelling with a one-way wave equation 4.1 Introduction The last decade has seen an increase in the quality and number of temporary regional broadband seismic experiments conducted over specific geological targets that have imaged the lithosphere in remarkable detail. Such exper- iments demonstrate that elastic properties of continental crust and the un- derlying mantle are both laterally heterogeneous and anisotropic at a range of scales. The characterization of these properties provides key constraints on lithospheric deformation and evolution. For example, teleseismic imag- ing of the lithosphere below Precambrian cratons has brought new evidence of early-Earth continental stabilization through shallow subduction by the identification of fossil structures within the lithospheric mantle [e.g. Pop- peliers and Pavlis, 2003; Mercier et al., 2008]. Other examples include the characterization of deep tectonic structures with multiple implications for modern lithospheric processes [e.g. Vinnik et al., 1995; Kosarev et al., 1999; Yuan et al., 2000; Bostock et al., 2002; Zandt et al., 2004]. Results such as these have prompted the funding of new experiments (e.g. INDEPTH), infrastructure (e.g. POLARIS), or both (e.g. USArray). In most analyses of teleseismic wavefields for lithospheric structure (e.g. A version of this chapter has been published. P. Audet, M. G. Bostock, and J.-P. Mercier. Teleseismic waveform modelling with a one-way wave equation. Geophysical Journal International, 171, 1212-1225, 2007. c©2007 Blackwell Publishing 70 Chapter 4. Teleseismic waveform modelling receiver functions, SKS splitting, traveltime tomography), one or both of two simplifying assumptions are commonly made. These assumptions are 1) that the wavefield propagation is adequately represented by ray theory, and 2) that structure can be (locally) characterized as one-dimensional. Con- sequently, interpretation of results in regions of complex geology can prove difficult or, worse, may yield erroneous results. In this context, an improved exploitation of high-quality broadband data sets for optimum recovery of the elastic tensor requires the use of comprehensive modelling tools. Unfor- tunately, fully three-dimensional complete techniques (e.g. finite differences [Virieux , 1986; Levander , 1988], spectral elements [Komatitsch and Tromp, 1999]) are still impractical for inverse modelling on desktop computers and may provide too much information, rendering the interpretation difficult. Ray tracing methods are intuitively and computationally appealing but may neglect important wave effects. A computationally attractive yet near-complete representation of the seismic wavefield is based on the solution of the one-way wave equation. In this approach one identifies a preferred direction of propagation and splits the wavefield into constituents travelling in opposite directions. The one- way formulation is a valid approximation in teleseismic studies because the wavefronts are propagating at near-vertical incidence and perpendicular to the predominantly horizontal layering of the Earth. The main advantages of one-way techniques are firstly the explicit separation between the effects of forward wave propagation (e.g. geometrical spreading and forward wave coupling) and back-scattering from variations in material properties, and secondly the promise of computational efficiency. The one-way approach thus renders straightforward the isolation of wave effects considered impor- tant in the analysis without having to model waveform complexities arising in full wavefield modelling, effectively reducing numerical costs. The ability to turn off reflections will be most advantageous in receiver function studies because it can sometimes be difficult to distinguish wave-type conversions arising at deeper interfaces from free-surface reflections arising at shallower ones, and in shear-wave splitting analysis of SKS phases where only the near-vertical forward propagation of the shear phase through anisotropic 71 Chapter 4. Teleseismic waveform modelling media is of interest. In this paper we discuss the implementation of a wide- angle one-way wave equation in a teleseismic context. 4.2 Theory The one-way method is based on the factorization of the generally anisotropic, elastic wave equation into two operators which are first order with respect to a given direction of propagation. One factor gives the one-way equation for waves travelling in the forward direction. In order to fully understand the develpment of the elastic one-way operators, it is helpful to first consider the acoustic case. 4.2.1 Acoustic one-way operators There has been a large number of studies concerning the acoustic approxima- tion of one-way wave equations, partly because the acoustic wave equation is scalar and thus its algebraic approximation is much simpler, and because it reasonably represents the physics of compressional wave (P -wave) propa- gation [McCoy et al., 1986; McCoy and Frazer , 1986; Thomson, 2005]. Consider a Cartesian coordinate system (x1, xα) in which x1 is the pre- ferred direction of wavefield propagation, and α = 2, 3. The acoustic, source- free, wave equation is given by the Helmholtz equation and is written as (neglecting density) ( ∂21 +H2(x1, xα, ∂α;ω) ) φ(x1, xα;ω) = 0, (4.1) where ∂1, ∂α are the spatial derivatives, φ is the acoustic wavefield, ω is the angular frequency, and H2 is the Helmholtz (two-way) operator given by H2(x1, xα, ∂α;ω) = ( ∂2α + ( ω c(x1, xα) )2) , (4.2) where c is the variable wavespeed of the medium. The ω dependency will be omitted through the remainder of this paper for brevity. In the case of a 72 Chapter 4. Teleseismic waveform modelling homogeneous medium the Helmholtz operator has a simple polynomial form in the Fourier domain given by H̃2(pα) = ω 2 ( 1 c2 − p2α ) , (4.3) where H̃2 is the Fourier representation of H2 and pα is the horizontal slow- ness component. H2 is quadratic in pα and is thus a Partial Differential Operator (PDO). The homogeneous, source-free, Helmholtz equation can be factorized into forward- and backward-going components (∂1 + iH1(∂α)) (∂1 − iH1(∂α))φ(x1, xα) = 0, (4.4) where H1 is the square-root operator that has a Fourier representation given by H̃1(pα) = ω √ 1 c2 − p2α. (4.5) The square-root operator is not a polynomial in the Fourier domain and thus it belongs to the class of Pseudo-Differential Operators (PSDO), which forms a generalized class of PDOs. The Fourier representation of a PSDO is given in terms of its symbol. In this form (4.5) H̃1 is the vertical component of the slowness vector. Retaining only one factor in (4.4) gives the acoustic homogeneous one-way wave equation. In a laterally varying medium, the Helmholtz operator is a complicated polynomial of the Fourier variable, and the square-root operator is usually found in an approximate manner. The most widely used is the parabolic approximation based on a Taylor expansion of H2 in the spatial domain [Claerbout , 1985]. The validity of this approach is limited by a series of conditions (in different combinations) concerning the medium parameters (weak inhomogeneities, weak gradients) and wavefield (narrow-angle) prop- agation. Although computationally efficient, this approximation has the ma- jor disadvantage of restricting the propagation to near-vertical angles, even in the case of a homogeneous medium. The phase-screen method decom- poses the square-root operator into a space-dependent term that accounts 73 Chapter 4. Teleseismic waveform modelling for the effects of homogeneous propagation at a spectrum of slownesses, and a wavenumber-dependent term that describes local lateral variations at a single slowness [Wu, 1994; Wild and Hudson, 1998; Hoop et al., 2000]. The method also breaks down for wide-angle propagation. The operators re- main exact in the case of a homogeneous medium. Other implementations include the exact but computationally expensive modal wavefield decompo- sition where H1 is obtained by expanding the Helmholtz operator into its eigenfunctions [Grimbergen et al., 1998], and the wide-angle asymptotic ap- proximation, where the one-way operators are found from a high-frequency expansion of the square-root operator using Fourier analysis [McCoy and Frazer , 1986; Thomson, 1999]. All the ideas presented above can be generalized to the elastic and gen- erally anisotropic one-way (vector) wave equation [e.g. Thomson, 1999]. For example, a parabolic approximation has been successfully implemented to model waveform distortions due to elastic heterogeneity and anisotropy at a low computational cost [Angus et al., 2004; Angus, 2005]. Thomson [1999] also derived an expression for the wide-angle propagation of elastic wave- fields based on plane-wave (Fourier) decomposition and a high-frequency approximation of the square-root operator. The leading-order term in fre- quency has a simple form in terms of the eigen-polarizations and slownesses for plane waves, and describes the forward phase propagation and coupling of the wave components. Higher-order terms describing the scattering due to inhomogeneities, and sometimes referred to as energy-flux normalization terms [e.g. Chapman, 2004], are left unexploited. An acoustic implementa- tion of the wide-angle, one-way method showed that physically motivated approximations greatly improved upon the numerical efficiency [Thomson, 2005]. It was found that higher-order terms are of utmost importance for the proper description of the differential transmission of energy across an inter- face. Here we build upon these findings and develop the elastic analogue of the acoustic implementation. The necessary theoretical background of the elastic wide-angle, one-way wave equation is provided in detail in the paper by Thomson [1999]. A brief summary is presented in the next subsection. 74 Chapter 4. Teleseismic waveform modelling 4.2.2 Elastic one-way operators The source-free elastic one-way wave equation for waves travelling in the forward (x1) direction in a homogeneous medium is given by (∂1 − iωD(∂α))u(x1, xα) = 0, (4.6) where D = GP1G −1 will be termed the propagator matrix, the 3 × 3 matrices G and P1 are respectively the 3-column eigen-displacements of the forward going waves and the diagonal matrix of corresponding eigenvalues of p1, and u is the 3-component displacement vector [Thomson, 1999]. The propagator matrix is the square-root operator, solution to the factorization of the homogeneous wave equation, and is the equivalent of equation (4.5) for an elastic medium. The theory described herein takes into account the full elastic tensor and thus readily allows for waveform modelling through the most general anisotropic media. The propagator matrix properly describes the phase shift propagation and mode coupling between the different wave types composing the forward-going wavefield. No assumptions have been made in the derivation of equation (4.6), and hence it is exact. When inhomogeneities are present, the differential operators are repre- sented in the Fourier domain with respect to the transverse coordinates, where Fourier analysis is used to obtain an approximate solution. The fac- torization of the inhomogeneous wave equation requires the solution of the square-root operator which is a PSDO that has a Fourier representation in terms of its symbol. In the remainder of the paper we will adopt the shorthand sym(D(xα, ∂α)) = D̃(xα, pα) ≡ D(xα, pα), and omit the x1 de- pendency unless otherwise specified. Since an exact representation of the symbol of the root operator is difficult to obtain, a high-frequency asymp- totic expansion is commonly applied [Thomson, 1999]. This approximation still allows for finite-frequency scattering because the signal is band-limited in practice, but restricts the length scales for observing material hetero- geneity [McCoy et al., 1986]. Only the first two terms of the expansion are retained and the one-way equation is written (appendix) 75 Chapter 4. Teleseismic waveform modelling ∂1u(xα) = ( ω 2π )2 ∫ ∫ ( iωD0(xα, pα) + iD1(xα, pα) +O(ω −1) ) · exp[iω(xα − yα)pα]u(yα) dyα dpα. (4.7) Thomson [1999] shows that the leading-order term in frequency is given by the propagator matrix derived in the homogeneous case D0(xα, pα) = G(xα, pα)P1(xα, pα)G(xα, pα) −1, extended by a decomposition/recomposition over transverse plane waves. To recover the wavefield, the leading-order one-way wave equation (c.f. equa- tion (4.7), omitting D1) performs a number of operations: (1) the first integration over yα decomposes u(yα) into local plane waves; (2) the ma- trix G−1 resolves the displacement into quasi-P, and quasi-S1, S2 plane modes; (3) the diagonal matrix P1 defines the rate of progression of each local plane mode along x1; (4) the matrix G reconstitutes the total plane wavefield after a step in that direction; and (5) the integration over pα re- constitutes the wavefronts. In this form the wide-angle propagator correctly models geometric ray theory, Maslov- and Kirchhoff-like representations of the elastic wavefield, because no restriction is made concerning the range of slownesses allowed [Thomson, 1999; Angus, 2007]. However, the solution of (4.7) requires time-consuming integrations in both spatial and wavenumber (or slowness) domains. Thomson [2005] showed that physically motivated assumptions can greatly improve upon numerical efficiency, as described in the next subsection. The second-order term represents the energy-flux normalization term that properly describes the amplitude scaling due to the lateral and depth dependent variations of the medium parameters. The exact form of this term requires the inversion of matrices and an increase in computational 76 Chapter 4. Teleseismic waveform modelling costs. We derive in the appendix an approximate factor that takes the form D1(xα, pα) ≈ − i 2 G(xα, pα)P ′ 1(xα, pα)G(xα, pα) −1 P′1(xα, pα) = ( (C11(xα)P1(xα, pα)) −1 ∂1 (C11(xα)P1(xα, pα)) ) (4.8) where C11 is, in the special case of isotropy, a diagonal matrix given by the inverse wavespeeds of the wave modes (see appendix). The approximate amplitude scaling factor (4.8) neglects the effects due to gradients of the medium parameters in the transverse direction, and is based on the narrow- angle approximation. However the one-way operators still describe lateral gradient effects through the dependence of all the terms on xα. This form was also chosen because it is implemented in a straightforward manner, as described in the following subsections. 4.2.3 Practical considerations Higher order extrapolation The one-way equation (4.7) defines the derivative of the wavefield in the preferred direction of propagation. The leading term can be differentiated with respect to x1 to yield the second derivative ∂ 2 1u accurate to within the one-way approximation [Thomson, 1999]. The extended propagator is then expanded as [eq. 5.9 Thomson, 1999] ǫ∂1u(xα) + 1/2ǫ 2∂21u(xα) = ( ω 2π )2 ∫ ∫ (iωD0(xα, pα)ǫ +1/2iω∂1D0(xα, pα)ǫ 2 − 1/2ω2D20(xα, pα)ǫ2) · exp[iω(xα − yα)pα]u(yα)dyαdpα +O(ǫ3). (4.9) A simplification occurs if the propagator is evaluated at the midpoint x1 + 1/2ǫ ≡ xǫ, and after re-introducing the x1 dependency the one-way 77 Chapter 4. Teleseismic waveform modelling wave equation can be written u(x2ǫ, xα) =( ω 2π )2 ∫ ∫ G(xǫ, xα, pα) exp[iωP1(xǫ, xα, pα)ǫ]G −1(xǫ, xα, pα) · exp[iω(xα − yα)pα]u(x1, yα)dyαdpα +O(ǫ3). (4.10) It can be seen from equation (4.10) that the conversion of wave types across an interface are taken up by a form of boundary condition for displace- ments through the productG−1(x+ǫ)G(x−ǫ). If the material properties vary across the interval [x−ǫ, x+ǫ], this matrix product introduces off-diagonal el- ements and generates converted waves that are then propagated with the exponential phase term. In the derivation of the above expression (4.10), it was assumed that lateral gradients are small, and that the stepsize ǫ is small enough for the Taylor expansion to be valid. However, Thomson [2005] suggests that larger stepsize may be used in practice, because the form (4.10) corresponds to the exact separable solution in a homogeneous medium. The stepsize can then be chosen to densely sample the medium where vertical and lateral inhomogeneities exist, and remain coarse elsewhere. Natural interpolation In the form (4.10), the one-way wave equation still appears to require a pro- hibitive amount of computation, because the Fourier synthesis necessitates (1) the decomposition of the wavefield into Fourier (plane waves) compo- nents over N lateral xα grid points (first use of the FFT with N log2N operations), and (2) the formation of a unique propagator, or the associated eigenvectors/ eigenvalues matrices (G,P1,G −1), and inverse transform at each lateral xα grid point (N use of the FFT, N 2 log2N operations). The propagation of the wavefield thus requires N + 1 FFT calls at each frequency. To reduce the costs of computations, a key observation [Thom- son, 2005] is made concerning the form of the one-way wave equation. For a laterally homogeneous (layered) model, the propagator obviously needs only 78 Chapter 4. Teleseismic waveform modelling be evaluated at a single lateral grid point, and so the modelling requires only two calls to the FFT to describe the propagation. If the model con- sists of only a few homogeneous regions described by M < N special points (or nodes) on the lateral grid, and if in between points are interpolated to give smooth wavespeed gradients, the propagator needs only be evaluated at those M points. This yields M coverings of the entire lateral grid, and the second integral becomes a weighted sum of the M partial contributions to the wavefield, where the weights are the interpolation coefficients calcu- lated for the Earth model. The number of FFT required to describe the propagation of the wavefield is thus reduced to M , and hence necessitates MN log2 N operations, i.e. computational costs are linearly proportional to M . See Thomson [2005] for details on computation. We note that the lat- erally varying regions represent variations of the full elastic tensor, and thus include variations in anisotropic parameters as well as variations in isotropic wavespeeds. Scattering operators Near vertical reflection/transmission (R/T) at horizontal boundaries can be addressed directly in the one-way approach. The transmission operator is simply the forward square-root operator given by (4.7). Explicitly, the transmitted wavefield is calculated by including the second order term D1 in the asymptotic expansion of the square-root operator, D, and takes the form uT (x1 + ǫ, xα) = ( ω 2π )2 ∫ ∫ G(xǫ, xα, pα)T(xǫ, xα, pα) · exp[iωP1(xǫ, xα, pα)ǫ]G−1(xǫ, xα, pα) · exp[iω(xα − yα)pα]u(x1, yα)dyαdpα, (4.11) where the energy-flux term is obtained by the integration of (4.8) over a small step ǫ in the forward direction following the procedure described previously, 79 Chapter 4. Teleseismic waveform modelling and is written T(xǫ, xα, pα) = (C11(x+ǫ, xα)P1(x+ǫ, xα, pα)) −1 · (C11(x−ǫ, xα)P1(x−ǫ, xα, pα))1/2 . (4.12) The reflection operator is not derived explicitly because the 3 × 3 for- mulation neglects the coupling between forward and backward going waves. However to the same order of approximation we can use the local plane-wave coefficients for reflection computed using the material properties across the interface at lateral position xα and slowness pα, because they represent the higher order terms in the asymptotic expansion of the corresponding pseudo- differential operators [McCoy et al., 1986; McCoy and Frazer , 1986]. The natural interpolation idea also applies to the scattering operators since the single-interface R/T is asymptotically local [Thomson, 2005]. The extension to a sequence of laterally inhomogeneous layers relates closely to the theory of invariant embedding [Kennett , 1983], where the exact reflection and transmission response is calculated for a stack of homogeneous layers. From an implementation perspective, the first-order effects of scattering are obtained by making 2 passes through the model; one in the forward direction and a second in the reverse direction. The reflection coefficients are stored during the first pass and added to the wavefield on the second pass [Thomson, 2005]. Free surface multiples can be included in the wavefield by calculating the local plane wave free surface reflection coefficients and applying them to the wavefield at the surface on the first pass. In the second pass the reflected wavefield at the surface is taken as the initial wavefield and marched in the reverse direction. On the third pass the reflection coefficients calculated on the second pass are added to the initial wavefield. 2.5D modelling Although the wide-angle one-way formulation allows for waveform modelling through 3D media, we restrict our applications to 2D models with invari- ant material properties in the perpendicular direction. For planar incident 80 Chapter 4. Teleseismic waveform modelling wavefields there is a single 2D horizontal slowness for any individual source. We can thus model 3D propagation through 2D media by giving a har- monic dependence to the wavenumber component of the incident wavefield in the perpendicular direction because it will scatter independently from the in-plane components. 4.3 Synthetic examples 4.3.1 1D crustal structure This example serves as a benchmark for the one-way implementation of tele- seismic wave propagation with the aim of validating the amplitude correction term (4.12) and first-order scattering by invariant embedding described in the previous section. The 1D model consists of a 36 km layer with typical crustal properties (VP = 6.55 km/s, VS = 3.70 km/s, ρ = 2800 kg/m 3) overlying a mantle-like half-space (VP = 8.10 km/s, VS = 4.50 km/s, ρ = 3500 kg/m3). The model has a total vertical extent of 100 km, and the vertical spacing in the extrapolation is varied from 15 km within the ho- mogeneous regions, to 250 m near the discontinuity. Because the model is laterally homogeneous, the natural interpolation idea states that the model needs only be sampled at one lateral grid point. Wavefield extrapolation is still carried out on the 2D grid and seismograms are obtained at each lateral grid point. The lateral extent is 384 km with 3 km spacing. We selected the seismogram at the center of the grid to compare with the reference solution. The input wavefield is a positively-polarized plane P-wave Ricker pulse coming from below with a dominant wavelength of 10 km (at mantle P ve- locities) and a horizontal slowness of 0.08 s/km, which is the maximum limit usually encountered in teleseismic studies [Bostock , 1998]. Time spacing is taken as 0.25 s and 129 frequencies are used. The propagated one-way wave- field is compared with a reference solution calculated with the reflectivity technique [Kennett , 1983; Thomson, 1996], which calculates the exact R/T response for a vertical stack of plane homogeneous layers. Figure 4.1 shows the horizontal and vertical component seismograms. 81 Chapter 4. Teleseismic waveform modelling −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Ve rti ca l d isp la ce m en t 0 5 10 15 20 25 30 Time (s) One−Way RMatrix PmS PPmP PSmP + PPmS PSmS −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 H or iz on ta l d isp la ce m en t 0 5 10 15 20 25 30 Time (s) One−Way RMatrix PmS PPmP PSmP + PPmS PSmS Figure 4.1: Comparison of the one-way simulations with the reflectivity technique Top: vertical component seismograms of an initial P-wave ricker pulse with a slowness of 0.08 [s/km]. Bottom: horizontal component seismo- grams. We artificially introduced a time offset of 0.25s to visually compare the one-way results with the reference solution. The one-way results prop- erly reproduce the Moho conversion (PmS) and the free surface multiples of the different phases (PPmP ,PPmS ,PSmP ,PSmS) and the exact arrival times, although the amplitudes are not fully corrected by the approximate one-way scattering term. 82 Chapter 4. Teleseismic waveform modelling The one-way results reproduce the P to S wave conversion generated from the velocity contrast at the Moho, termed PmS , and the free-surface re- verberations and conversions of the different phases PPmP , PSmP+PPmS, PSmS . The modelling reproduces the relative arrival times exactly, although the amplitudes are not fully corrected by our approximate amplitude nor- malizing factor (e.g. PSmS). This example demonstrates one of the limits of the wide-angle approach. Increasingly accurate results are found for wavefields propagating at smaller slownesses. To illustrate this effect we calculated the relative difference in amplitude for increasing slowness for both the horizontal and vertical components (Figure 4.2). The amplitude difference increases quasi linearly with increasing slowness. The amplitude difference is lower for the horizontal component and it reaches zero at normal incidence. The amplitude difference reaches 4% for the vertical component at a slowness of 0.08. This fit is reasonable considering the large horizontal slowness of the incident plane wave, the large velocity contrast, and the level of approximation in obtaining the amplitude correction term (see appendix). Note that computation cost is independent of slowness. Computation time: ∼10 s using single precision Intel FC compiler on a 2.2 GHz AMD Athlon(tm) 64 processor. 4.3.2 2D subduction model The simple subduction model consists of a downgoing oceanic lithospheric crust and mantle, overriding continental crust and mantle, and sub-lithospheric oceanic mantle. The geometry and velocities are based on a simplification of the Cascadia subduction zone model of [McNeill , 2005] (Figure 4.3). The oceanic crust is defined with a thickness of approximately 6 km. The conti- nental Moho is taken at a depth of 37 km at the eastern edge of the study area. The subduction zone model is shown as gray shaded contours in Fig- ure 4.3. The model is defined on a 360 × 110 km2 Cartesian grid. The lateral and vertical node spacings are equal to 500 m, or about 20 points per dominant wavelength of the initial P-wave in the sub-slab mantle. Ma- 83 Chapter 4. Teleseismic waveform modelling 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 Ac cu ra cy 0.00 0.02 0.04 0.06 0.08 Slowness [s/km] Vertical Horizontal Figure 4.2: Accuracy of the energy-flux normalization term versus slowness. Gray and black curves show the relative amplitude difference between one- way simulations and reflectivity technique for the vertical and horizontal components, respectively. The accuracy decreases for increasing slowness, which is expected due to the narrow-angle assumption in the derivation of the energy-flux term (appendix A). terial properties are defined in Table 4.1. The same geometry is used in the one-way and reference calculations. The reference solution is calculated using the pseudo-spectral method of Kosloff et al. [1984], which calculates the spatial derivatives of the full wave equation with the use of the Fast Fourier Transform, instead of by differencing (Finite-Difference) or through interpolation functions (Finite- Element). Time derivatives are approximated by differencing and the time Structural el. P vel. [km/s] S vel. [km/s] Density [kg/m3] Suboceanic mantle 7.7 4.5 3190 Oceanic slab mantle 8.1 4.6 3340 Oceanic slab crust 6.8 3.9 2900 Subcontinental mantle 7.5 4.1 3260 Continental crust 6.4 3.6 2770 Table 4.1: Material properties of the structural elements composing the 2D subduction zone model. 84 Chapter 4. Teleseismic waveform modelling 0 25 50 75 100 D ep th [k m] 50 100 150 200 250 300 350 Horizontal distance [km] Figure 4.3: 2D subduction zone model defining the geometry of the struc- tural elements: subducting oceanic slab comprising the crust and litho- spheric mantle, the sublithospheric oceanic upper mantle, and the overriding continental crust and upper mantle. Velocities are defined in table 4.1. spacing is taken to 0.02 s with a total of 1500 increments. The pseudo- spectral method is able to accurately resolve frequencies up to 3.6 Hz. The initial wavefield is a normally incident plane P-wave Ricker pulse with a dominant wavelength of 10 km at sub-slab mantle P velocity. Although larger vertical stepsizes in the one-way method can be specified where the material properties in the vertical direction are slowly varying, we maintained the same regular spacing to facilitate direct comparison with the reference solution. The lateral model sampling is specified by the layer discontinuities on a cross-section profile of the model at a given depth. The wavefield is interpolated between adjacent model points as permitted by the natural interpolation idea. One-way extrapolation is performed using the exponential extrapolator for 129 frequencies, with a final seismogram time increment of 0.25 s. In this modelling only the first pass through the model is computed to capture the essential features of single-scattering and avoid the cluttering of the images due to a high degree of model complexity. Figure 4.4 shows the comparison of horizontal component seismograms from (a) the one-way and (b) the pseudo-spectral method. Amplitudes are scaled according to unit-normalized initial P-wavefield. Four major features are present in both simulations: (1) the direct transmitted P arrival re- 85 Chapter 4. Teleseismic waveform modelling 0 5 10 15 20 25 30 Ti m e [s] 0 50 100 150 200 250 300 350 Horizontal Distance [km] 3 31 2 4 −0.05 0.00 0.05 5 10 15 20 25 30 Ti m e [s] 50 100 150 200 250 300 350 Horizontal Distance [km] 3 31 2 4 5 −0.05 0.00 0.05 Figure 4.4: Horizontal-component synthetic wavefields calculated using the one-way algorithm (left) and the pseudo-spectral method (right). Colorbar shows the amplitude normalized to the initial incident P -wave. Only one pass is included in the modelling and the seismograms show zeroth-order scattering effects. The initial wavefield is a normally incident plane wave Ricker pulse with a dominant wavelength of 10 km. Four main features are common to both simulations: (1) positive PS conversion from the dipping Moho; (2) mixed polarities PS conversions from the various interfaces of the subducting slab; (3) diffraction patterns associated with edge discontinuities and Fourier wrap-around effects; and (4) diffracted wavefronts caused by the point scatterer at the tip of the mantle wedge. The later arrivals (5) seen in the pseudo-spectral simulation are caused by artificial reflections due to a velocity contrast between the surface and bottom velocities of the model box. The pseudo-spectral simulations further show multiple internal reflections beneath normal continental Moho. Faint arrivals at t = 0− 15 s at 100-270 km distance in the one-way synthetics are due to wrap-around effects from time-domain Fourier reconstruction. 86 Chapter 4. Teleseismic waveform modelling fracted at the dipping Moho that disappears beneath normal continental lithosphere; (2) the mixed polarity PS conversions at the various dipping in- terfaces; and (3) the bleeding tails of the wavefields due to artificial vertical discontinuities at the edges of the box and Fourier wrap-around effects. An- other noteworthy feature seen in the seismograms is the diffraction pattern (4) caused by the point scatterer at the tip of the mantle wedge, although this feature seems to have a longer duration in the pseudo-spectral syn- thetics. Moreover the one-way simulations do not reproduce the wavefronts associated with the internal reflections of the direct P wave and converted S waves (dipping positive and negative pulses east of km 250, at 14 s and onwards) bouncing off the multiple slab interfaces beneath continental litho- sphere, since reflections have been neglected in this realization. The pseudo- spectral simulation shows further back scattering effects at later times that we chose not to model with the one-way algorithm. Those are caused by artificial reflections at the free surface due to a velocity contrast between sur- face and bottom velocities of the model box. Apart from those later effects, the agreement between the two methods for the forward-scattered portion of the wavefield is remarkably good considering the level of approximation of the one-way algorithm and the complexity of the velocity model. Figure 4.5 shows wiggle plots of both horizontal and vertical compo- nent seismograms at discrete horizontal locations. Horizontal-component one-way synthetics are less noisy than the pseudo-spectral results due to 1) the absence of internal reflections in the one-way modelling, and 2) larger numerical instabilities in the pseudo-spectral algorithm caused by the dis- crete lateral grid sampling. The main features seen in horizontal-component seismograms is discussed above. The vertical-component seismograms show the horizontal moveout of the first arrival due to higher velocities beneath the subducting lithosphere, compared to lower velocities beneath the conti- nental lithosphere. From 6 s onwards the amplitudes are scaled by a factor of 10 to emphasize the finer details of the seismograms. Very low-amplitude free-surface reflections are again observed at ∼9 s on the pseudo-spectral seismograms that are absent from the one-way simulations. Reflections from the vertical boundary of the model box appear in both cases although the 87 Chapter 4. Teleseismic waveform modelling 0 5 10 15 20 25 30 0 18 36 54 72 90 108 126 144 162 180 198 216 234 252 270 288 306 324 342 Time [s] H or iz on ta l d ist an ce [k m] Horizontal component 0 5 10 15 20 25 30 0 18 36 54 72 90 108 126 144 162 180 198 216 234 252 270 288 306 324 342 Time [s] H or iz on ta l d ist an ce [k m] Vertical component Figure 4.5: Wiggle plots of horizontal (left) and vertical (right) component seismograms obtained with the one-way (black lines) and the pseudo-spectral (gray lines) methods. Horizontal-component one-way seismograms are less noisy than the pseudo-spectral results, and the main arrivals are coherent over the whole time range. Vertical-component seismograms show a move- out due to lower velocities beneath normal continental mantle. Vertical amplitudes have been scaled up 10 times starting at 6 s to emphasize the finer details of the seismograms. 88 Chapter 4. Teleseismic waveform modelling amplitudes vary between the two methods (e.g. traces 18 and 19). Ab- sorbing one-way boundary conditions could in theory be implemented to avoid artificial reflections from vertical boundaries [Clayton and Engquist , 1977]. An alternative procedure is to extend the model horizontally with constant values outside the rectangle box such that the artifacts arrive at much later times and do not contaminate the synthetic wavefields. We apply this procedure in the next section for modelling real data. Note that non-vertically incident waves can be modelled in both cases with similar levels of accuracy. Normal incidence was chosen because the waveforms are more easily interpreted in terms of scattering due to structure heterogeneity. Computation time: one-way ∼15 min using single precision Intel FC compiler on a 2.2 GHz AMD Athlon(tm) 64 processor; pseudo-spectral about the same as the one-way method on the same machine. 4.3.3 Remarks on computational efficiency In our demonstration no attempt has been made to optimize the numerical efficiency of the algorithm, i.e. we kept a regular grid and only applied the wavefield interpolation idea. The flexibility offered by using variable step- sizes in the forward direction, which would most improve upon the numerical efficiency, was not exploited to retain highest possible accuracy. Even so, the one-way algorithm performs equally fast as the pseudo-spectral technique. As intuition suggests, the computational improvements will be most signif- icant if the contrasts in material properties in the transverse and vertical directions are small [Thomson, 2005]. 4.4 Application to data from NW Canada A seismic array composed of nearly 60 broadband instruments with ∼35 km spacing was recently deployed in northwestern Canada as part of the IRIS- PASSCAL/Lithoprobe CAnada NOrthwest Experiment (CANOE) for a period of 15 to 27 months. A subset of the array stretching over the Wop- 89 Chapter 4. Teleseismic waveform modelling 30° 60° 90° 120° 220˚ 220˚ 230˚ 230˚ 240˚ 240˚ 250˚ 250˚ 260˚ 260˚ 55˚ 55˚ 60˚ 60˚ 65˚ 65˚ Edmonton Whitehorse Yellowknife Fort−Nelson NW TerritoriesYukon B. Columbia Albt. Sask. Figure 4.6: Data and station array from the CANOE seismic experiment in northwestern Canada. Left: Events used in the generation of teleseismic receiver functions. Right: Subset of stations used in the analysis composed of 20 stations with ∼ 15 km spacing. may orogen was more densely sampled with station spacing of ∼15 km (Fig- ure 4.6). This subset of 20 stations allowed teleseismic study of a Proterozoic fossil subduction zone previously characterized by Lithoprobe (SNOR- CLE) reflection/refraction surveys [Cook and Erdmer , 2005]. Figure 4.6 shows the location and events used in the study. Data were processed using standard receiver function techniques [Bostock , 1998]. Receiver functions are source-normalized seismograms obtained by deconvolving the P compo- nent from the SV (shear radial) and SH (shear transverse) components, and characterize the scattering structure of the crust and upper mantle under- lying seismic stations. Details of the methodology and interpretation are described in the paper by Mercier et al. [2008]. Individual receiver functions are ordered from left to right along the pro- file based on the geographic station locations, and, for individual stations, on the sampling point of the PS conversions projected onto the west-east directions. Radial- and transverse-component receiver function images (Fig- ure 4.7, bottom two panels) of the subsurface beneath this region shows a complex pattern of seismic discontinuities. The radial receiver function im- age (Figure 4.7-bottom left) shows a clear and nearly continuous Moho at ∼4 s, represented by a red pulse with one significant, 75 km long interruption near the suture at km 150. In the same area, the presence of crustal ma- 90 Chapter 4. Teleseismic waveform modelling terial dipping into the mantle is also inferred. The geometry of subduction is better defined, however, on the transverse component image where the subducting plate is associated with a ∼10 km thick layer exhibiting strong elastic anisotropy. This layer clearly extends from depths of 35 km (∼4 s) at the suture to at least 70 km (∼9 s), over a distance of 150 km further eastward, with some indication that it dips still more steeply thereafter to depths of >120 km. A 1-θ correction has been applied to the transverse component receiver functions to improve the visual coherence of the dipping structure [Mercier et al., 2008]. Previous reflection/refraction surveys imaged strong dipping reflectors at the suture between the two terranes down to a depth of ∼60 km that are interpreted to be a ∼10 km thick subducting crustal layer [Cook and Erdmer , 2005]. Surprizingly, but quite definitively, the anisotropic layer imaged with teleseismic receiver functions does not coincide with the crustal reflectors but, rather, parallels that feature. The top of the anisotropic layer corresponds to the base of the crustal layer, that is, the Moho of the subducting plate. The top of the crustal layer is not evident in the teleseismic image. The PS conversion from the base of the anisotropic layer exhibits a frequency dependence favouring low frequencies which suggests that it is better modelled as a first-order (vs zeroth-order) discontinuity (or linear gradient) in material properties. The argument favoring anisotropy versus dipping isotropy deserves some discussion because it provides insights into the nature of the anomalous structure. Coherent energy on the transverse component is the consequence of lateral heterogeneity, anisotropy, or both. However the presence of energy on the transverse component directly beneath the horizontal portion of the Moho at distances less than 150 km cannot be produced by an isotropic interface.The absence of strong conversion from the lower (basal) interface of the dipping layer on the radial component suggests that average veloci- ties of the layer and the underlying mantle are comparable. Moreover, the azimuthal correction applied to the transverse component image does not correspond to the signature of an isotropic layer dipping to the east. Note that we were not able to reproduce the main features on the transverse 91 Chapter 4. Teleseismic waveform modelling Radial Component Sy nt he tic s ei sm og ra m s Ti m e [s] 0 2 4 6 8 10 12 14 R ea l s ei sm og ra m s Ti m e [s] 100 150 200 250 0 2 4 6 8 10 12 14 Transverse Component 100 150 200 250 0.25 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 0.25 Distance along the profile [km] Am pl itu de re la tiv e to P ? Figure 4.7: Modelled (top) and observed (bottom) teleseismic receiver func- tions (RF) of the CANOE seismic experiment. Left panels show the radial component RF; right panels are the transverse component RF. The same set of slownesses and backazimuths was used to generate the one-way synthetic data. A butterworth filter is applied to the RF to emphasize frequency- dependent effects. The radial component shows strong signals at < 2 s due to reverberations within the sedimentary cover that are not modelled in the simulations. On the same component, a clear continuous positive (red) pulse associated with PS conversion at the Moho is visible at ∼ 4 s with one interruption near a suture between two Paleo-Proterozoic terranes, and the presence of crustal material dipping into the mantle. The geometry is better defined, however, on the transverse component where the dipping layer exhibits strong elastic anisotropy and first-order frequency effects from a gradual decrease in the percentage of anisotropy from top to bottom of the slab. The one-way simulations reproduce the main features seen in the data. 92 Chapter 4. Teleseismic waveform modelling 0 25 50 75 100 125 D ep th [k m] 50 100 150 200 250 Horizontal distance [km] Figure 4.8: Subsurface model of northwestern Canada showing the geometry of the dipping layer, represented as a series of 5 thin layers exhibiting a gradual linear decrease in anisotropy from top to bottom with an average mantle velocity (grey shades represent anisotropic layering, not velocities). The maximum amplitude of hexagonal anisotropy is 5 % at the top of the layer and the fast axis of seismic wavespeed has a plunge of 51o within the plane of the slab and a trend of 54o calculated from the north. The plunge and trend are rotated along the downgoing portion of the slab to account for the layer dip. Subslab mantle and crustal P-velocities are taken to 8.1 and 6.6 km/s, and S-velocities are those calculated for a Poisson solid (α = √ (3)β). The subcontinental wedge shows a gradual increase in velocity from 7 km/s to 8.1 km/s to the eastern edge of the model box. component using either a low- or high-velocity isotropic layer. This data set offers a convenient opportunity to test the one-way algo- rithm for various material properties including heterogeneity, anisotropy, as well as physical wave phenomena such as frequency-dependent scattering. In order to reproduce the main features of the receiver functions we model the lithospheric structure as a ∼10 km thick layer of dipping anisotropic material (average VP = 8.1 km/s) overlying a homogeneous mantle (VP = 8.1 km/s, isotropic Poisson material), overlain by a 37 km thick homogeneous crust (VP = 6.6 km/s) (Figure 4.8). We model the interruption and dipping of the Moho on the radial component by a subcrustal mantle wedge showing a linear horizontal gradient in velocity from VP =7 km/s at the suture to 8.1 km/s at the eastern boundary. We note that purely linear horizontal gra- dients need not be sampled at the same rate as dipping discontinuities, since 93 Chapter 4. Teleseismic waveform modelling only two points are necessary to model a linear gradient. This is reminiscent of the natural interpolation idea described in section 2.3.2. The first order discontinuity within the slab is modelled as a stack of thin layers showing a decrease from strong (5%) to weak anisotropy (hexagonal symmetry) from top to bottom. The fast axis of seismic wavespeed has a plunge that follows the layer dip and a trend of 220 degrees. Note that the natural interpolation idea is applied to the full elastic tensor, and not to the anisotropic parameters (symmetry axis, strength, etc.). Model discretization is 128 × 220 with spatial sampling of 2 km laterally, and 500 m vertically. Vertical grid spacing is increased to 5 km in the homogeneous crustal section. The three-component seismograms are picked according to stations location along the regular 2D profile. In order to produce results consistent with the receiver functions, we used the same events (in terms of slownesses and backazimuths) as input to the algorithm and applied the same processing as the raw data to produce similar receiver functions. As before the initial wavefields are pure plane P-wave Ricker wavelets with a dominant wavelength of 10 km. Real and synthetic receiver functions are plotted in Figure (4.7) in bottom panels. Traces are filtered between 0.2 and 1.5 Hz. The dipping anisotropic layer model clearly accounts for the essential features seen in the receiver functions, although finer details are left un- modelled. The radial component receiver function shows the discontinuous positive PS Moho conversion, and the presence of material dipping into the mantle. The dipping layer is obviously identified on the transverse compo- nent and shows the frequency-dependent signature of a velocity gradient, in this case a gradient in anisotropy. The implications of this exercise are twofold: (1) the one-way algorithm is able to reproduce all phenomena re- lated to material heterogeneity, anisotropy, and gradients; (2) the presence of dipping anisotropic material at lithospheric mantle depths suggests that shallow subduction was active in the early Proterozoic and that this process may have helped in the stabilization of subcontinental mantle root. The latter interpretation is discussed in detail in the paper by Mercier et al. [2008]. 94 Chapter 4. Teleseismic waveform modelling 4.5 Conclusion We implement a waveform modelling tool using the elastic one-way wave equation to model wave effects generated by structure heterogeneity and anisotropy in a teleseismic context. The form of wave propagation (near planar wavefronts) and scattering operators affords a number of physically motivated approximations to improve upon the computational efficiency of the algorithm. We show through synthetic modelling that the algorithm is able to reproduce waveforms from simple 1D and 2D heterogeneous models, including free surface reflections. We then proceed to model data from an array of seismic stations deployed in northwestern Canada that show a dipping layer of anisotropic material interpreted to be a fossilized Proterozoic subduction structure. We show that the algorithm reproduces the main features of the data. The elastic one-way algorithm could be used in other teleseismic appli- cations, such as SKS splitting analysis in complicated settings, simulating for example corner flow behind a subduction front [Levin et al., 2007]. Its relatively low computational cost makes it a good candidate for Monte Carlo inversion of first order scattering of teleseismic wavefields. The model imple- mentation and sampling is readily expanded in three dimensions to simulate wave propagation in more complicated media. These applications will be the focus of future studies. 4.6 Acknowledgments This paper was much improved by thorough and constructive reviews by two anonymous reviewers. 95 Bibliography Angus, D. A., A one-way wave equation for modelling seismic waveform variations due to elastic heterogeneity, Geophys. J. Int., 162, 882–898, 2005. Angus, D. A., True amplitude corrections for a narrow-angle one-way elas- tic wave equation, Geophysics, 72, T19–T26, 2007. Angus, D. A., C. J. Thomson, and R. G. Pratt, A one-way wave equation for modelling variations in seismic waveforms due to elastic anisotropy, Geophys. J. Int., 156, 595–614, 2004. Bostock, M. G., Mantle stratigraphy and evolution of the Slave province, J. Geophys. Res., 103, 21,183–21,200, 1998. Bostock, M. G., R. D. Hyndman, S. Rondenay, and S. M. Peacock, An inverted continental Moho and serpentinization of the forearc mantle, Na- ture, 417, 536–539, 2002. Chapman, C., Fundamentals of Seismic Wave Propagation, Cambridge University Press, Cambridge, UK, 2004. Claerbout, J., Imaging the Earth’s interior, Blackwell Sci. Publ., Oxford, United Kingdom (GBR), 1985. Clayton, R., and B. Engquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bull. Seism. Soc. Am., 67, 1529–1540, 1977. Cook, F. A., and P. Erdmer, An 1800 km cross section of the lithosphere through the northwestern North American plate: lessons from 4.0 billion years of Earth’s history, Can. J. Earth Sci., 42, 1295–1311, 2005. 96 Bibliography Grimbergen, J. L. T., F. J. Dessing, and K. Wapenaar, Modal expansion of one-way operators in laterally varying media, Geophysics, 63, 995–1005, 1998. Hoop, M. V. D., J. H. L. Rousseau, and R.-S. Wu, Generalization of the phase-screen approximation for the scattering of acoustic waves, Wave Mo- tion, 31, 43–70, 2000. Kennett, B. L., Seismic wave propagation in stratified media, Cambridge University Press, Cambridge, UK, 1983. Komatitsch, D., and J. Tromp, Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophys. J. Int., 139, 806– 822, 1999. Kosarev, G., R. Kind, S. V. Sobolev, X. Yuan, W. Hanka, and S. Oreshin, Seismic evidence for a detached Indian lithospheric mantle beneath Tibet, Science, 283, 1306–1309, 1999. Kosloff, D., M. Reshef, and D. loewenthal, Elastic wave calculations by the fourier method, Bull. Seis. Soc. Am., 74, 875–891, 1984. Levander, A., Fourth-order finite-difference P-SV seismograms, Geophysics, 53, 1425–1436, 1988. Levin, V., D. Okaya, and J. Park, Shear wave birefringence in wedge-shaped anisotropic regions, Geophys. J. Int., 168, 275–286, doi: 10.1111/j.1365– 246X.2006.03,224.x, 2007. McCoy, J., and L. N. Frazer, Propagation modelling based on wavefield factorization and invariant embedding, Geophys. J. R. astr. Soc., 86, 703– 717, 1986. McCoy, J., L. Fishman, and L. N. Frazer, Reflection and transmission at an interface separating transversely inhomogeneous acoustic half-spaces, Geophys. J. R. astr. Soc., 85, 543–562, 1986. 97 Bibliography McNeill, A. F., Contribution of post-critical reflections to ground motions from mega-thrust events in the Cascadia subduction zone, Master’s thesis, University of Victoria, BC, Canada, 2005. Mercier, J.-P., M. G. Bostock, P. Audet, J. B. Gaherty, E. J. Garnero, and J. Revenaugh, The teleseismic signature of fossil subduction: Northwestern Canada, J. Geophys. Res., 113, B04,308, doi:10.1029/2007JB005,127, 2008. Poppeliers, C., and G. Pavlis, Three-dimensional, prestack, plane wave mi- gration of teleseismic P -to-S converted phases: 2. Stacking multiple events, J. Geophys. Res., 108, 2267, doi:10.1029/2001JB001,583, 2003. Thomson, C. J., Notes on waves in layered media to accompany program Rmatrix, in Seismic waves in complex 3-D structures, pp. 147–162, De- partment of Geophysics, Charles University, 1996. Thomson, C. J., The ’gap’ between seismic ray theory and ’full’ wavefield extrapolation, Geophys. J. Int., 137, 364–380, 1999. Thomson, C. J., Accuracy and efficiency considerations for wide-angle wavefield extrapolators and scattering operators, Geophys. J. Int., 163, 308–323, 2005. Vinnik, L. P., R. W. E. Green, and L. O. Nicolaysen, Recent deformations of the deep continental root beneath southern Africa, Nature, 375, 50–52, 1995. Virieux, J., P-SV wave propagation in heterogeneous media; velocity-stress finite-difference method, Geophysics, 51, 889–901, 1986. Wild, A. J., and J. A. Hudson, A geometrical approach to the elastic complex screen, J. Geophys. Res., 103, 707–725, 1998. Wu, R.-S., Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, J. Geophys. Res., 99, 751–766, 1994. 98 Bibliography Yuan, X., et al., Subduction and collision processes in the central andes constrained by converted seismic phases, Nature, 408, 958–961, 2000. Zandt, G., H. Gilbert, T. J. Owens, M. Ducea, J. Saleeby, and C. H. Jones, Active foundering of a continental arc root beneath the southern Sierra Nevada in California, Nature, 431, 41–46, 2004. 99 Chapter 5 Morphology of the Explorer/Juan de Fuca slab edge in northern Cascadia: Imaging plate capture at a ridge-trench-transform triple junction 5.1 Introduction The northern end of the Cascadia subduction zone is characterized by the interaction of three oceanic plates (Pacific - PAC, Juan de Fuca - JdF, and Explorer - EXP) with North America (NA), where they meet in a set of poorly defined triple junctions [Braunmiller and Nábelek , 2002]. The EXP plate is a small oceanic fragment that detached from the subducting JdF plate around 4 Ma during creation of the left-lateral Nootka transform fault located offshore Nootka Island, south of Brooks Peninsula [Rohr and Furlong , 1995]. Since EXP inception the EXP ridges have been rotating clockwise and spreading centers have jumped northward, slowly reducing convergence with NA [Braunmiller and Nábelek , 2002]. The EXP region A version of this chapter was accepted for publication. P. Audet, M. G. Bostock, J.-P. Mercier, and J. F. Cassidy. Morphology of the Explorer/Juan de Fuca slab edge in northern Cascadia: Imaging plate capture at a ridge-trench-transform triple junction. Geology, 2008, in press. 100 Chapter 5. Morphology of slab edge in northern Cascadia Figure 5.1: Identification of major tectonic boundaries and terranes in west- ern Canada. BP Brooks Peninsula, BPfz - Brooks Peninsula fault zone, NI - Nootka Island, QCTJ - Queen Charlotte Triple Junction. Dotted lines de- lineate extinct boundaries or shear zones. Seismic stations are displayed as inverted black triangles. Station projections along LINE 1 and LINE 2 are plotted as thick white lines. White triangles represent Alert Bay volcanic field centers. Center of array locates the town of Woss. is currently undergoing major internal shear deformation as an ephemeral adjustment accommodating relative motion between the PAC and NA plates [Rohr and Furlong , 1995; Dziak , 2006]. EXP micro-plate evolution is an active example of an important process associated with ridge collision and subduction that has occured through- out geological history, typically leaving little geological records [Stock and Lee, 1994]. Given the strong dependence of plate reconstructions on oceanic data, it is critical that we understand processes associated with micro-plate formation. To gain insight into the causes of plate fragmentation, we need 101 Chapter 5. Morphology of slab edge in northern Cascadia to consider the factors controlling the balance between driving and resistive forces acting at plate boundaries. In a subduction setting the most effective force driving plate motion is the pull exerted by the subducted slab [Forsyth and Uyeda, 1975; Govers and Meijer , 2001]. Hence much could be learned about the creation and fate of the EXP plate by exploring its structure at crustal and upper mantle levels beneath northern VI and landward. Unfor- tunately, little is known about slab structure in northern Cascadia. Indirect evidence for location of the slab edge from heat flow, gravity and geochemi- cal data loosely define its surficial projection along a N-E trending corridor landward and parallel to BP [Lewis et al., 1997]. However those data lack resolution along the third (depth) dimension which is necessary to constrain the morphology of the subducting slab. Cassidy et al. [1998] used teleseismic data recorded at an array of 5 stations sparsely deployed across northern VI to provide the first direct ev- idence for a dipping low-velocity zone (LVZ) up to BP, and the resumption to normal continental-like seismic signature to the north. Here we use seis- mic data from a recently deployed POLARIS portable array together with data from permanent stations to map morphology of the subducting slab in the forearc region beneath northern VI and the interior of western British Columbia (BC) using receiver functions and P-wave traveltime tomography. 5.2 Shallow structure from receiver functions The POLARIS-BC Northern Vancouver Island (NVI) array comprises 26 seismometers along two mutually perpendicular arms (Fig. 5.1). One arm trends NW-SE in a direction parallel to strike and straddles the assumed northern end of the subduction zone (LINE 1), and the second arm trends SW-NE in the dip direction, just north of the extension of the NF beneath VI where subduction is observed (LINE 2). The array has been in opera- tion since June 2005 and to date each station has recorded an average of 50 teleseismic events with high signal-to-noise ratio. The NVI array is com- plemented with a few permanent stations run by the Geological Survey of Canada and from a previous experiment [Cassidy et al., 1998]. 102 Chapter 5. Morphology of slab edge in northern Cascadia All receiver functions (filtered between 0.05 and 0.35 Hz) are plotted in raw form along each line according to station position and, for individual stations, sorted by back-azimuth (Fig. 5.2). These data are most sensitive to structures with scale lengths between 1s and to 10s km and are dominated across both profiles by the signature of a LVZ. This signature comprises three sets of oppositely polarized pulses representing forward scattered P- to-S (Ps) (3-5 s at station VI57) conversions, and back-scattered P-to-S (Pps) (9-14 s at VI57) and S-to-S (Pss) (15-20 s at VI57) conversions af- forded by reflection of the teleseismic wavefield at the Earths free surface. These signals are interpreted as oceanic crust of the subducting JdF-EXP slab consistent with its expression in studies both further south beneath VI [Nicholson et al., 2005], Oregon [Bostock et al., 2002], and worldwide [e.g. Abers et al., 2006]. Along LINE 1 the LVZ is evident from VI57 to VI52, and disappears further north. The seismic response there is more similar to that of a single discontinuity at a typical continental crust-mantle boundary (Moho) with positive Ps and Pps arrivals at ∼4 s and ∼15 s, and a negative Pss pulse at ∼22 s. A LVZ is inferred at station PHC, although it lacks clear Ps signals and its relation to subducted oceanic crust to the south is unclear. Structural signals along LINE 2 show evidence of a well-defined LVZ dipping gently NE along the entire profile which is most easily seen in oppositely polarized reverberations (Pps, Pss) at ∼10-20 s. Oppositely polarized Ps signals are clearly imaged from stations CHM to VW01, whereafter they show polarity cross-overs as seen by the blue-red checkerboard pattern at ∼3 s from stations VI09 to VI07. This polarity reversal with respect to backazimuth of incident wavefield manifests elastic anisotropy. This change in the degree of elastic symmetry also coincides with a shallowing of the LVZ signature signaled by earlier arrivals of reverberated phases (Pps, Pss). Dipping LVZ signals resume thereafter to station VI01 where the signature disappears. The timing of scattered modes relative to P (0 time in Fig. 5.2) can be used to characterize both thickness and average VP /VS of the overlying column by assuming a dipping planar geometry of the subsurface. This is accomplished by stacking waveforms of the three scattering modes with time 103 Chapter 5. Morphology of slab edge in northern Cascadia 0 10 20 30T im e [se c] −0.3 −0.2−0.1 0.00.1 0.2 0.3 NW a) SE PH C VI3 2 VI3 1 NIM VI3 0 VI5 0 VI5 2 VI5 3 WO S VI0 8 VI0 9 VI5 4 VI5 5 VI5 6 VI5 7 0 10 20 30T im e [se c] −0.3 −0.2−0.1 0.00.1 0.2 0.3 SW b) NE CH M VI1 2 VW 03 VW 02 VW 01 VI0 9 WO S VI0 8 VI0 7 VI0 6 VI0 5 VI0 4 VI0 3 VI0 2 VI0 1 0 20 40 60D ep th [k m] NW c) SE PH C VI3 1 VI3 0 VI5 0 VI5 2 VI5 3 WO S VI0 9 VI5 4 VI5 5 VI5 6 VI5 7 0 20 40 60D ep th [k m] 0 20 40 60 80 100 120 140 Distance along line [km] SW d) NE CH M VI1 2 VW 03 VW 02 VW 01 WO S VI0 7 VI0 6 VI0 5 VI0 4 VI0 3 VI0 2 VI0 1 Figure 5.2: (a,b) Raw receiver functions along LINE 1 (a) and LINE 2 (b) sorted by station location and, for each station, by back-azimuth of incident wavefield. Distances between adjacent stations are not preserved. Inset in b) shows an enlargment of the first 6 seconds of the first Ps conversion beneath Woss to illustrate the polarity reversals that are diagnostic of anisotropy. (c,d) Best-fit estimates of depths to top (blue line) and bottom (red line) of subducted oceanic crust and, where present, continental crust-mantle discontinuity (Moho - black line) along both lines. Error bars are calculated as in Zhu and Kanamori [2000]. 104 Chapter 5. Morphology of slab edge in northern Cascadia delays that correspond to the propagation of plane waves through a range of models [Zhu and Kanamori , 2000]. Using this technique we determined the depth to top and bottom of the LVZ and mapped the morphology of the subducted oceanic crust along both lines (Fig. 5.2) and across northern Cascadia using a larger data set (Fig. 5.3). 5.3 Deep structure from P-wave tomography Deep slab structure (50 to 300 km depth) is imaged using P-wave travel- time tomography of the upper mantle in the Pacific Northwest using the method of Bostock and VanDecar [1995]. For this component of the study, data were collected over a much broader network comprising broadband and short period seismic stations from several Canadian and US portable and permanent arrays (Fig. 5.3). A total of 13,595 P-wave traveltimes were de- rived from 738 earthquakes, with good source coverage in back-azimuth and epicentral distance. The data set was inverted for upper mantle structure below Washington and western BC, extending the coverage much further north than a previous study Bostock and VanDecar [1995]. Resolution tests performed with a synthetic slab model extending along the full margin of northwestern NA indicate that deep structure is well resolved by the data. The velocity model is characterized by a quasi-planar, high-velocity layer steeply dipping beneath BC to depths of at least 300 km (Fig. 5.3). This material is inferred to represent the thermal and compositional anomaly of the subducted JdF plate. The slab signature appears to be continuous with a slight change in strike from Washington to northern VI where the signature ends abruptly. 5.4 Discussion Both shallow and deep slab results are plotted in Fig. 5.3. We will consider five main features: 1) the depth contours of oceanic crust outline a gently N-E-dipping underthrusting plate in northern VI at shallower levels than further south [see also Clowes et al., 1997]; 2) the top of oceanic crust shal- 105 Chapter 5. Morphology of slab edge in northern Cascadia 130oW 128oW 126oW 124oW 122oW 120 oW 47oN 48oN 49oN 50oN 51oN 52oN 53oN A1 A2 B1 B2 . % of P velocity −2 0 2 20 30 40 50 EXP JdF Figure 5.3: P-wave tomography image of upper mantle structure at 200 km depth. The blue high-velocity body represents the thermal and compo- sitional anomaly of the subducting JdF slab. Black dashed line indicates northern limit of subduction. White dashed line shows the approximate lo- cation of EXP detachment from JdF (Fig. 5.4). Cross-sections of the model along lines A1-A2 and B1-B2 are shown in Appendix B. Shallow oceanic crust contours from receiver functions are overlaid as black lines with depth to the top of the oceanic crust (in km) indicated. Inverted white triangles are broadband seismic stations used in receiver functions. Inverted grey tri- angles are additional short-period and broadband seismic stations used in traveltime tomography. 106 Chapter 5. Morphology of slab edge in northern Cascadia lows over a region centered on the town of Woss where anisotropy is present; 3) the contours deepen at the extrapolated location of the NF beneath VI; 4) the oceanic crustal signature disappears sharply north of BP; and 5) the edge of the deep slab is imaged at all depths, indicating continuity in structure from VI into the interior of BC. 5.4.1 Further evidence of a sharp slab edge The imaged morphology of the slab in northern VI is in agreement with all available geological and geophysical data in the region. Heat flux measure- ments show a transition from low values (∼46 mW/m2) over the subducted portion of the slab in central and southern VI to higher values (∼67 mW/m2) to the north of BP over a distance of 50 km [Lewis et al., 1997]. This shift also coincides with a transition from low to high Bouguer anomaly, indicat- ing hotter but denser material at crustal depths north of BP [Lewis et al., 1997]. Moreover, this region shows subdued topography, lower mean ele- vation, and the absence of deep seismicity, suggesting the absence of slab [Lewis et al., 1997]. Of these constraints, heat flow is the most informative regarding the time evolution of the slab margin due to the 4̃0 Ma necessary to estab- lish the surficial heat flow transition through the conductive crustal regime [Lewis et al., 1997]. The large heat flow difference in northern VI roughly coincides with the imaged slab edge, and suggests a stable location of the Queen Charlotte triple junction at BP during that time [Lewis et al., 1997]. Stability of the ridge-trench-transform triple junction at the northern end of Cascadia for ∼40 Ma and subduction of the JdF ridge likely imply toroidal mantle flow around the slab edge which can cause thermo-mechanical erosion and slab melting of the feather slab edge, generating a complex geochemical environment at the surface [Thorkelson and Breisprecher , 2005]. In northern VI the imaged edge of the slab coincides with the Alert Bay volcanic field, which possesses a geochemical signature resembling that of ocean-floor basalts and within-plate basalts and is distinct from volcanic rocks expected within the subduction zone forearc [Armstrong et al., 1985]. 107 Chapter 5. Morphology of slab edge in northern Cascadia The slab edge also coincides with a change in fault orientation northwest of the BP fault zone (Fig. 5.1), in agreement with periods of extension in the Tertiary, which, assuming triple junction stability, is presumably caused by the viscous coupling of the overriding NA crust with infilling asthenospheric material [Lewis et al., 1997]. The location of the deep slab edge beneath the interior of BC also coincides with the disappearance of normal Cascade arc volcanism. 5.4.2 Explorer tectonics: Slab stretching? The current state of EXP kinematics (EXP is converging with NA [Braun- miller and Nábelek , 2002; Mazzotti et al., 2003] or EXP is a new transform boundary between PAC and NA [Rohr and Furlong , 1995; Kreemer et al., 1998; Dziak , 2006]), is a conundrum that we cannot address directly in this study. The creation and fate of EXP, however, are governed by dynamical processes driven by the balance of driving to resistive forces acting at plate boundaries that are intimately linked with EXP structure. In a subduction zone, the ratio of negative buoyancy of the subducting slab (slab pull) to resistive forces acting on the plate drives plate motion [Forsyth and Uyeda, 1975; Govers and Meijer , 2001]. Buoyancy of the sub- ducting slab is controlled mainly by temperature and, to a lesser extent, composition. At a ridge-trench triple junction the subducted slab is young and hot, and its increased buoyancy diminishes slab pull, and hence conver- gence. However the net slab pull exerted along the much larger Cascadia margin dominates over diminished pull at the edge and the subducted ridge is forced to plunge down underneath the continent. The ensuing toroidal flow around subducted slab edge generates hot mantle upwellings and con- tributes to increased temperature near the edge, and thus further increases its buoyancy. It is then likely that the feather slab edge will remain at shal- lower levels than further south of the triple junction, consistent with the observed oceanic crustal depth contours of the EXP-JdF slab beneath VI (Fig. 5.3). These factors hold the key to the formation of EXPmicro-plate. Creation 108 Chapter 5. Morphology of slab edge in northern Cascadia of the Nootka fault implies that EXP started resisting subduction around 4 Ma due to increased slab buoyancy at shallow levels. Elevated tempera- tures most likely prevented the Nootka fault from tearing through the entire subducted JdF slab. Hence a possible cause of the creation of the EXP micro-plate is the stretching/tearing of the slab by tensile forces to accom- modate plate reconfiguration at shallower levels. Subducted slab stretching is the analogue of crustal thinning during continental rifting, with tensile forces provided by slab pull [ten Brink et al., 1999]. Assuming a linear rel- ative velocity between the EXP and JdF slabs with time from 0 to 2 cm/a, we estimate the separation to cause ∼40 km of slab stretching. Ensuing mantle upwelling results in the thermal uplift of the slab at the location of maximum extension, consistent with an increase in Bouguer anomaly and a shallowing of the subducted oceanic crust. We postulate that the shallow portion of the oceanic crust centered around Woss represents the expression of stretching of the EXP slab at subcrustal levels, and accounts for the localized strong anisotropic fabric presumably due to shearing (Fig. 5.4). Note that slab contours deepen near the extension of the NF beneath VI at station VI57, indicating a disruption in slab continuity. We interpret the subducted portion of the EXP plate to be the small underthrust segment bordered by the NF, the shallow swell near Woss where maximum stretching is inferred, and the BP fault zone. This model implies that the JdF slab remains unperturbed to the north and east of the EXP slab and still contributes to slab pull, and is consistent with the reduced convergence of the EXP plate north of the NF and EXP plate capture by NA. Such episodes of plate capture have been reported in similar ridge-trench triple junction settings, e.g. the Rivera plate north of the Cocos plate in Central America [DeMets and Traylen, 2000], and the Monterey and Arguello plates offshore Baja California [Stock and Lee, 1994; Zhang et al., 2007]. Ongoing convergence of the EXP plate with NA is probably due to a combination of viscous coupling between the EXP and JdF slab downdip and along the creeping NF beneath VI. 109 Chapter 5. Morphology of slab edge in northern Cascadia Figure 5.4: Cartoon showing a tectonic interpretation of subducted slab edge morphology in northern Cascadia. Inset shows an areal view of the region. We postulate that the shallowing of oceanic crust centered on Woss represents the expression of slab stretching. The EXP slab is detached from the JdF plate along the NF and the region of inferred stretching. 110 Chapter 5. Morphology of slab edge in northern Cascadia 5.5 Conclusion Resolving the creation and fate of the EXP micro-plate in northern Cas- cadia requires characterization of the structural elements controlling plate dynamics at a triple junction. To reconcile both models of the EXP plate, we propose a tectonic model based on seismic and geophysical data in which complexity in slab morphology is attributed to thermo-mechanical erosion of the slab edge caused by ridge subduction and toroidal mantle flow, in- creased buoyancy, reduced slab pull, slab stretching, and detachment from the JdF plate along the NF. Our model implies that JdF subduction is still active north and east of the detached EXP slab, and that the EXP plate is being captured by NA. 5.6 Acknowledgments We are grateful to Ken Dueker and George Zandt for access to BATHOLITHS data. We acknowledge thoughtful reviews by Rob Govers and two anony- mous reviewers. This work is supported by NSERC. This is Geological Survey of Canada publication 2008021. 111 Bibliography Abers, G. A., P. E. van Keken, E. A. Kneller, A. Ferris, and J. C. Stachnik, The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow, Earth Planet. Sci. Lett., 241, 387–397, doi:10.1016/j.epsl.2005.11.055, 2006. Armstrong, R. L., J. E. Muller, J. E. Harakal, and K. Muehlenbachs, The Neogene Alert Bay Volcanic Belt of northern Vancouver Island, Canada: Descending-plate-edge volcanism in the arc-trench gap, J. Vol- canol. Geotherm. Res., 26, 75–97, 1985. Bostock, M. G., and J. C. VanDecar, Upper mantle structure of the north- ern Cascadia subduction zone, Can. J. Earth Sci., 32, 1–12, 1995. Bostock, M. G., R. D. Hyndman, S. Rondenay, and S. M. Peacock, An inverted continental Moho and serpentinization of the forearc mantle, Na- ture, 417, 536–539, 2002. Braunmiller, J., and J. Nábelek, Seismotectonics of the Explorer region, J. Geophys. Res., 107, 2208, doi:10.1029/2001JB000,220, 2002. Cassidy, J. F., R. M. Ellis, C. Karavas, and G. C. Rogers, The northern limit of the subducted Juan de Fuca plate system, J. Geophys. Res., 103, 26,949–26,961, 1998. Clowes, R. M., D. J. Baird, and S. A. Dehler, Crustal structure of the Cascadia subduction zone, southwestern British Columbia, from potential field and seismic studies, Can. J. Earth Sci., 34, 317–335, 1997. 112 Bibliography DeMets, C., and S. Traylen, Motion of the Rivera plate since 10 Ma relative to the Pacific and North American plates and the mantle, Tectonophysics, 318, 119–159, 2000. Dziak, R. P., Explorer deformation zone: Evidence of a large shear zone and reorganization of the Pacific-Juan de Fuca-North American triple junction, Geology, 34, 213–216, doi:10.1130/G22,164.1, 2006. Forsyth, D., and S. Uyeda, On the relative importance of the driving forces of plate motion, Geophys. J. R. Astron. Soc., 43, 163–200, 1975. Govers, R., and P. T. Meijer, On the dynamics of the Juan de Fuca plate, Earth Planet. Sci. Lett., 189, 115–131, 2001. Kreemer, C., R. Govers, K. P. Furlong, and W. E. Holt, Plate boundary deformation between the Pacific and North America in the Explorer region, Tectonophysics, 293, 225–238, 1998. Lewis, T. J., C. Lowe, and T. S. Hamilton, Continental signature of a ridge-trench-triple junction: Northern Vancouver Island, J. Geophys. Res., 102, 7,767–7,781, 1997. Mazzotti, S., H. Dragert, J. Henton, M. Schmidt, R. Hyndman, T. James, Y. Lu, and M. Craymer, Current tectonics of northern Cas- cadia from a decade of GPS measurements, J. Geophys. Res., 108, 2554, doi:10.1029/2003JB002,653, 2003. Nicholson, T., M. G. Bostock, and J. F. Cassidy, New constraints on sub- duction zone structure in northern Cascadia, Geophys. J. Int., 161, 849– 859, 2005. Rohr, K. M. M., and K. P. Furlong, Ephemeral plate tectonics at the Queen Charlotte triple junction, Geology, 23, 1035–1038, 1995. Stock, J. M., and J. Lee, Do microplates in a subduction zones leave a geological record?, Tectonics, 13, 1471–1487, 1994. 113 Bibliography ten Brink, U., I. Shimizu, and P. Molzer, Plate deformation at depth under northern California: Slab gap or stretched slab?, Tectonics, 18, 1084–1098, 1999. Thorkelson, D., and K. Breisprecher, Partial melting of slab window mar- gins: Genesis of adakitic and non-adakitic magmas, Lithos, 79, 25–41, 2005. Zhang, X., H. Paulssen, S. Lebedev, and T. Meier, Surface wave to- mography of the Gulf of California, Geophys. Res. Lett., 34, L15,305, doi:10.1029/2007GL030,631, 2007. Zhu, L., and H. Kanamori, Moho depth variation in southern California from teleseismic receiver functions, J. Geophys. Res., 105, 2969–2980, 2000. 114 Chapter 6 Seismic evidence for overpressured subducted oceanic crust and sealing of the megathrust 6.1 Introduction In recent years, it has become increasingly apparent that significant quanti- ties of water (H2O) are stored in the mantle wedge as the serpentine mineral antigorite [Hyndman and Peacock , 2003]. Compelling evidence for the pres- ence of antigorite comes from seismic observations that reveal the presence of highly resolved, low shear-velocity anomalies in both warm (young slab - e.g. Cascadia [Bostock et al., 2002]) and cold (old slab - e.g. northeast Japan [Kawakatsu and Watada, 2007]) subduction zone environments in locations predicted by thermal and petrological models [Hyndman and Peacock , 2003]. Water in the mantle wedge most likely originates from the dehydration of hydrous minerals within subducting oceanic crust and mantle, which will occur at shallower depths (¡50 km) in warmer environments [Peacock , 1990; Kirby et al., 1996]. Hydrous minerals form within oceanic lithosphere (1-2% H2O in the crust [Peacock , 1990]) as a result of hydrothermal circulation and alteration at spreading ridges, along fracture zones, and in the trench-outer rise region. There is little constraint, however, on the in-situ abundance A version of this chapter was submitted for publication. P. Audet, M. G. Bostock, N. I. Christensen, and S. M. Peacock. Seismic evidence for overpressured subducted oceanic crust and sealing of the megathrust. Submitted to Nature, 2008. 115 Chapter 6. Overpressured subducted oceanic crust and distribution of H2O in the fluid and hydrous minerals within subduct- ing oceanic crust, even though this knowledge is likely to be important in assessing the strength of the plate interface and nature of seismogenesis. Recently, studies of episodic tremor and slip (ETS) in both Cascadia and Japan point to an origin that involves, either directly or indirectly, fluids in the general vicinity of the plate boundary, downdip of the locked zone and near the intersection with forearc mantle [Kao et al., 2005; Shelly et al., 2006; Wang et al., 2006]. Characterizing the physical state of subducted oceanic crust and distribution of fluids therein may help in understanding the temporal and spatial occurrence of ETS events. 6.2 Method We employ observations of converted teleseismic waves (or receiver func- tions) recorded over a POLARIS portable array of seismometers across Van- couver Island at approximately 10 km spacing (Fig. 6.1). All receiver func- tions from this experiment are plotted in raw form from west-south-west to east-north-east according to station position and, for individual stations, according to their sampling points (Fig. 6.2A). These data are sensitive to structures with scale lengths of 1-10 km and are dominated across the en- tire profile by the signature of an eastward dipping low-velocity zone. This signature includes 3 sets of oppositely polarized pulses comprising forward scattered P-to-S (Ps, 3-4 s at station PGC) conversions, and back-scattered P-to-S (Pps, 13-15 s at PGC) and S-to-S (Pss, 17-20 s at PGC) conversions afforded by reflection of the teleseismic wavefield at the Earth’s free surface. These signals have been previously interpreted to represent oceanic crust of the subducting Juan de Fuca plate [Nicholson et al., 2005], consistent with its expression in studies further south beneath Oregon [Bostock et al., 2002], and worldwide [Yuan et al., 2000; Abers et al., 2006; Kawakatsu and Watada, 2007]. In this study we re-examine the signature of the different scattered modes in more detail. It is widely appreciated [Zandt and Ammon, 1995; Zhu and Kanamori , 2000] that the timing (relative to the incident P-wave arrival, 116 Chapter 6. Overpressured subducted oceanic crust 233˚ 233˚ 234˚ 234˚ 235˚ 235˚ 236˚ 236˚ 237˚ 237˚ 238˚ 238˚ 48˚ 48˚ 49˚ 49˚ 50˚ 50˚ Vancouver Island Georgia Strait Vancouver Victoria Seattle Figure 6.1: Geographical map of northern Cascadia subduction zone. Loca- tion of the three-component, broadband seismic stations that recorded data used in this study are shown as inverted red triangles. Note that the array is oriented approximately perpendicular to the strike of subduction, which trends N45W. 117 Chapter 6. Overpressured subducted oceanic crust 0 time in Fig. 6.2A) of the scattered modes from a given interface can be used to constrain both depth and average P to S velocity ratio (VP /VS) of the overlying column. It can also be shown (Appendix B) that relative times between two arrivals from distinct boundaries within the lithospheric column for each set of scattered modes can be used to similarly characterize the interval properties between the two boundaries, largely independent of overlying structure, even for dipping layers. This recognition allows us to investigate the velocity structure (i.e. VP /VS) of the downgoing oceanic crust. The key constraint is the ratio of the direct conversion delay time (∆TPs, the difference in time between direct conver- sions from the top and bottom of the oceanic crust) to the corresponding back-scattered quantities (∆TPps, ∆TPss). In the two extremes as VP /VS tends to 1 and ∞, ∆TPs will tend to 0 and ∆TPps (∆TPss/2), respectively. Care must be taken in measuring these quantities, since low-pass filtering at corner periods greater than twice the true separation between pulses will bias measurements of ∆TPs, ∆TPps, ∆TPss to greater values (Appendix B). 6.3 Results Fig. 6.2B shows our estimates of oceanic crustal VP/VS beneath each sta- tion of the array. These estimates are extremely high and it becomes more practical to speak in terms of Poisson’s ratio (σ, for which there is an upper bound of 0.5 for inviscid fluids) around 0.4. These values cannot be due to compositional effects (i.e., mineralogy within the downgoing oceanic crust). The common oceanic crustal rocks basalt, diabase and gabbro have σ be- tween 0.28 and 0.29. At the prevailing temperatures and pressures (∼1 GPa, ∼500 C), these rocks will have been metamorphosed to upper amphibolite facies leading to a drop in σ to 0.26 or 0.27 [Christensen, 1996]. Antigorite, the high-pressure variety of serpentine, has a σ of 0.29 [Christensen, 2004]. Even the most extreme candidates, the low pressure serpentine minerals chrysotile and lizardite [Christensen, 1996], exhibit σ below 0.37, and their presence can be further ruled out on the basis of pressure, temperature, and compositional considerations [Rondenay et al., 2008]. The only plausible ex- 118 Chapter 6. Overpressured subducted oceanic crust 0 5 10 15 20 25 30 Ti m e (se c) a PsT PpsT PssB PsB PpsBPssT −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 PFBTWG B TWB B TSJBLZBTWK B MGC B KELBPGC SSIB SILB GOW B SNB 0.25 0.30 0.35 0.40 0.45 0.50 Po is so n’ s ra tio b 1.73 1.87 2.08 2.45 3.32 vp /v s PFB TWG B TWB B TSJB LZB TWK B MGC B KELB PGCSSIBSILBGOW BSNB −2 −1 0 1 2 (∆ v s /v s) v s* (km /s) 150 160 170 180 190 200 210 220 230 240 250 Distance from trench (km) c Figure 6.2: Radial receiver function (RF) results in northern Cascadia. (A) RFs for all data filtered between 0.05 and 0.5 Hz used in the analysis sorted by station position along the array and, for each station, by back-azimuth of incident wavefield. Amplitudes are relative to incident P-wave. RFs are measures of local discontinuities in S-wave velocity. Red and blue col- ors correspond to velocity increase and decrease with depth, respectively. Forward- and back-scattered phases show oppositely-polarized pulses that are manifest of a dipping, low-velocity layer, indicated by blue and red la- bels. (B) Estimated VP /VS and Poisson’s ratio for the low-velocity layer. Errors are standard deviations from bootstrap estimates. Elevated values (∼0.4) are interpreted as pore fluid pressures near lithostatic conditions. (C) S-velocity perturbations relative to a one-dimensional reference model (V ∗S ) using a Born approximation for radial Ps phases. Positive (red) and negative (blue) velocity contrasts from the bottom and top, respectively, of the low-velocity layer are of similar amplitude. 119 Chapter 6. Overpressured subducted oceanic crust planation is that of high pore fluid pressures. Widely documented high pore fluid pressures existing at shallower levels in accretionary prisms are due pri- marily to porosity reduction [Moore and Vrolijk , 1992] but it is more likely that our observations manifest mineral dehydration [Christensen, 1984]. Un- fortunately, a precise correspondence between σ and pore fluid pressure is rendered difficult due to the lack of laboratory measurements at confining pressures corresponding to 20-40 km depth. Nonetheless, extrapolation of measurements made at lesser pressures [Christensen, 1984] render it likely that the oceanic crust beneath Vancouver Island is overpressured (near litho- static), a result with a range of important implications. 6.4 Discussion Evidence for high pore fluid pressure in the Nankai subduction zone has been presented recently [Shelly et al., 2006; Kodaira et al., 2004] but these studies have found less extreme values of VP /VS using different forms of traveltime tomography. VP /VS measurements derived from tomography are likely to be biased toward normal values through the blurring effects of regu- larization, simultaneous earthquake relocation, and the inherent difficulties in imaging low velocity zones with traveltimes [Gerver and Markushevitch, 1966; Wielandt , 1987]. Estimates of σ obtained using scattered waves are likely to be better resolved owing to their greater sensitivity to local velocity changes. Our documentation of high fluid pressure implies that the plate interface, defined by the top of the oceanic crustal layer, likely represents a low per- meability boundary. Using Darcys law, we estimate the permeability of the interface to be ∼5×1025 to ∼5×1022 m2, assuming pore-fluid pressures equal lithostatic pressures beneath the interface, a 0.65 GPa pressure drop across the interface corresponding to hydrostatic pore-fluid pressure above the in- terface at a depth of 35 km, a slab fluid production rate of 10−4 m3/(m2 yr) [Hyndman and Peacock , 2003], an interface thickness of 1 to 1000 m, and a dynamic viscosity of 10−4 Pa s for H2O at 500 C and 1 GPa. These estimates are lower than measurements of basaltic oceanic basement permeabilities in 120 Chapter 6. Overpressured subducted oceanic crust the range 10−21 to 10−11 m2 made at length scales of 1 to 1000 m [Fisher , 1998]. The significance of the calculation is reinforced by noting that the greatest permeabilities are found within the upper few hundred meters of oceanic crust and increase with scale length [Fisher , 1998]. The low per- meability of the plate interface may reflect deformation-induced grain-size reduction [Caine et al., 1996] or the precipitation of minerals from migrating fluids [Kato et al., 2003; Meneghini and Moore, 2007]. Further evidence supporting a low permeability plate boundary at crustal levels stems from scattered wave amplitudes which are sensitive to velocity contrasts at interfaces (Fig. 6.2C). Amplitudes of arrivals for a given scat- tering mode originating from the bottom and top oceanic crustal boundaries are comparable, indicating similar S-velocity jumps of±1 km/s, respectively. The lower boundary jump is consistent with that expected for the oceanic Moho which marks the juxtaposition of lower velocity gabbro with higher velocity mantle peridotite. The lower crust of the overriding plate is thought to be formed of dominantly mafic accreted terranes [Ramachandran et al., 2005], and would thus be expected to afford little or no contrast in seis- mic velocity with oceanic crust. Prominent amplitudes of scattered waves from the upper boundary thus support a sharp contrast from low-velocity oceanic crust containing trapped fluids to higher velocity overriding crust characterized by much lower pore fluid pressure. Low permeability of the plate boundary at depths <35 km contrasts with the seismic expression of forearc mantle serpentinization that appears downdip at depths >40 km in this region [Bostock et al., 2002; Nicholson et al., 2005] and which suggests that H2O is readily transported into the wedge (Fig. 6.3). A change in the nature of the plate boundary seal can be explained by the major (>10%) volume reduction and release of H2O that accompanies eclogitization [Ahrens and Schubert , 1975], the last significant metamorphic facies transition experienced by oceanic crust. Both this reac- tion and the still greater volume expansion accompanying serpentinization of the mantle wedge [Coleman, 1971] could be expected to enhance perme- ability in the vicinity of the plate boundary through fracture generation. This interpretation is consistent with the seismic signature in the cold sub- 121 Chapter 6. Overpressured subducted oceanic crust Figure 6.3: Schematic interpretation of receiver function results. Updip limit of overpressured oceanic crust is unconstrained. Station location is approximate. duction zone setting of NE Japan, where a serpentinite layer is interpreted to form above subducting oceanic crust at onset of eclogitization near ∼90 km depth [Kawakatsu and Watada, 2007]. It is worth noting that in warm subduction zones, ETS tends to occur in the vicinity of the plate interface near but seaward of the wedge corner [Kao et al., 2005; Shelly et al., 2006]. Our documentation of high pore- fluid pressures in this region supports two groups of models that involve fluids: 1) ETS is triggered by hydrofracturing at the plate interface [Wang et al., 2006]; and 2) ETS occurs where high pore-fluid pressures extend the region of conditionally stable slip [Shelly et al., 2006; Kodaira et al., 2004]. A downdip change in permeability across the plate interface could be explained by hydrofracturing of the seal where hydraulic conductivities become sufficiently high near the onset of eclogitization. Moreover, it might also represent an important factor in altering the mode of slip along the thrust within the frictional stability transition zone. Our model suggests the possibility of linking the recurrence of ETS at regular time intervals to periodic cycles of steady pore-fluid pressure build-up from dehydration of 122 Chapter 6. Overpressured subducted oceanic crust subducted oceanic crust, fluid release from fracturing of the interface during ETS, and subsequent precipitation sealing of the plate boundary. Finally, our documentation of high-pore fluid pressure in the subducting plate has important implications for the interpretation of tomographic ve- locity images, and the position of the plate boundary interface. High pore fluid pressures result not only in elevated σ (and VP /VS) but also depressed absolute velocities (-20% and -30% for P-waves and S-waves, respectively, in the measurements of Christensen [1996]). Consequently, identification of a tomographic iso-velocity contour, that does not account for this effect, as proxy for the oceanic Moho [Ramachandran et al., 2005] will yield interface depths that are significantly overestimated [see Nicholson et al., 2005]. 6.5 Acknowledgments We thank E. Davis and M. Jellinek for discussions. Data used in this study come from the Canadian National Seismological Network and are distributed freely by the Geological Survey of Canada. 123 Bibliography Abers, G. A., P. E. van Keken, E. A. Kneller, A. Ferris, and J. C. Stachnik, The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow, Earth Planet. Sci. Lett., 241, 387–397, doi:10.1016/j.epsl.2005.11.055, 2006. Ahrens, T. J., and G. Schubert, Gabbro-eclogite reaction rate and its geo- physical significance, Rev. Geophys., 13, 383–400, 1975. Bostock, M. G., R. D. Hyndman, S. Rondenay, and S. M. Peacock, An inverted continental Moho and serpentinization of the forearc mantle, Na- ture, 417, 536–539, 2002. Caine, J. S., J. P. Evans, and C. B. Foster, Fault zone architecture and permeability structure, Geology, 24, 1025–1028, 1996. Christensen, N. I., Pore pressure and oceanic crustal seismic structure, Geophys. J. R. astr. Soc., 79, 411–423, 1984. Christensen, N. I., Poisson’s ratio and crustal seismology, J. Geophys. Res., 101, 3139–3156, 1996. Christensen, N. I., Serpentinites, Peridotites, and Seismology, Int. Geol. Rev., 46, 795–816, 2004. Coleman, R. G., Petrologic and geophysical nature of serpentinites, Geol. Soc. Am. Bull., 82, 897–918, 1971. Fisher, A. T., Permeability within basaltic oceanic crust, Rev. Geophys., 36, 143–182, 1998. 124 Bibliography Gerver, M., and V. Markushevitch, Determination of a seismic wave veloc- ity from the travel time curve, Geophys. J. Roy. astron. Soc., 11, 165–173, 1966. Hyndman, R. D., and S. M. Peacock, Serpentinization of the forearc mantle, Earth Planet. Sci. Lett., 212, 417–432, doi:10.1016/S0012–821X(03)00,263– 2, 2003. Kao, H., S.-J. Shan, H. Dragert, G. Rogers, J. F. Cassidy, and K. Ra- machandran, A wide depth distribution of seismic tremors along the north- ern Cascadia margin, Nature, 436, 841–844, 2005. Kato, A., A. Sakaguchi, S. Yoshida, and H. Mochizuki, Permeability mea- surements and precipitation sealing of basalt in an ancient exhumed fault of a subduction zone, Bull. Earthquake Res. Inst., Univ. of Tokyo, 78, 83–89, 2003. Kawakatsu, H., and S. Watada, Seismic evidence for deep-water trans- portation in the mantle, Science, 316, 1468–1471, 2007. Kirby, S. H., E. R. Engdahl, and R. Denlinger, Intermediate-depth in- traslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs (overview), in Geophys. Monogr. 96, edited by G. E. Bebout, D. W. Scholl, S. H. Kirby, and J. P. Platt, pp. 195–214, American Geophysical Union, Washington, DC, 1996. Kodaira, S., T. Iidaka, A. Kato, J.-O. Park, T. Iwasaki, and Y. Kaneda, High pore fluid pressure may cause silent slip in the Nankai Trough, Sci- ence, 304, 1295–1298, 2004. Meneghini, F., and J. C. Moore, Deformationand hydrofracturing in a sub- duction thrust at seismogenic depths: The Rodeo Cove thrust zone, Marin Headlands, California, Geol. Soc. Am. Bull., 119, 174–183, 2007. Moore, J. C., and P. Vrolijk, Fluids in accretionary prisms, Rev. Geophys., 30, 113–135, 1992. 125 Bibliography Nicholson, T., M. G. Bostock, and J. F. Cassidy, New constraints on sub- duction zone structure in northern Cascadia, Geophys. J. Int., 161, 849– 859, 2005. Peacock, S. M., Fluid processes in subduction zones, Science, 248, 329–337, 1990. Ramachandran, K., S. E. Dosso, G. D. Spence, R. D. Hyndman, and T. M. Brocher, Forearc structure beneath southwestern British Columbia: A three-dimensional tomographic velocity model, J. Geophys. Res., 110, B02,303, doi:10.1029/2004JB003,258, 2005. Rondenay, S., G. A. Abers, and P. E. van Keken, Seismic imaging of subduction zone metamorphism, Geology, 36, 275–278, doi:10.1130/G24,112A.1, 2008. Shelly, D. R., G. C. Beroza, and S. Ide, Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip, Nature, 442, 488–491, 2006. Wang, Z., D. Zhao, O. P. Mishra, and A. Yamada, Structural heterogeneity and its implications for the low frequency tremors in Southwest Japan, Earth Planet. Sci. Lett., 251, 66–78, doi:10.1016/j.epsl.2006.08.025, 2006. Wielandt, E., On the validity of the ray approximation for interpreting delay times, in Seismic tomography, edited by G. Nolet, pp. 85–98, Reidel, 1987. Yuan, X., et al., Subduction and collision processes in the central andes constrained by converted seismic phases, Nature, 408, 958–961, 2000. Zandt, G., and C. J. Ammon, Continental crust composition constrained by measurements of crustal Poisson’s ratio, Nature, 374, 152–154, 1995. Zhu, L., and H. Kanamori, Moho depth variation in southern California from teleseismic receiver functions, J. Geophys. Res., 105, 2969–2980, 2000. 126 Chapter 7 Conclusions In this Chapter we summarize key contributions of this thesis and discuss current limitations. We also outline future research opportunities to further develop ideas presented in the thesis. 7.1 Main contributions The topic of this thesis is the characterization of structure and deformation of the lithosphere at past and present convergent continental margins in order to study specific tectonic processes. The approaches used here fit into two categories: 1) Potential field analysis; and 2) Teleseismic imaging. We advocate the novel use of a wavelet transform in potential field analysis to study large-scale deformation of the lithosphere. At a smaller scale, we employ observations of scattered teleseismic waves and develop a waveform modelling tool to image lithospheric structure in subduction zone forearcs to determine the evolution of subducted plates. The reported results on synthetic data show that the methods used are well-suited for our purposes and their application to real data allows us to build an understanding of convergent margin evolution and subduction processes. Each individual contribution is detailed in the following sections. 7.1.1 Wavelet analysis of potential field data We use a wavelet transform to compute the local and azimuthal variations of the coherence between Bouguer gravity and topography. The isotropic coherence is calculated by averaging the wavelet spectra from optimally over- lapping 2D Morlet wavelets having an isotropic spectral envelope in adjacent 127 Chapter 7. Conclusions directions within 180◦, defining the so-called “fan” wavelet. The isotropic wavelet coherence spectrum is inverted to obtain local estimates of the elastic thickness (Te) of the lithosphere. We calculate the anisotropic coherence by restricting the fan wavelet over an azimuthal range of 90◦. The direction of maximum coherence is diagnostic of the direction of preferred isostatic com- pensation, or the direction where the lithosphere is weakest. The coherence is inverted using the theoretical response of a thin anisotropic plate model. We have carried out extensive tests on synthetic topography and Bouguer gravity data sets to verify that: (1) the wavelet method can recover Te for simple models with either homogeneous or spatially variable rigidity pat- terns; and that: (2) the method can determine azimuthal variations in the 2D coherence for homogeneous models with anisotropic Te. 7.1.2 Large-scale deformation of convergent margins The wavelet method is used in conjunction with the window-based mul- titaper technique to estimate Te and its anisotropy in western Canada. In addition to a nearly stepwise change in Te from the Cordillera to the Craton, we show that weak axes of Te anisotropy are parallel with most compressive directions of horizontal stress components and fast axes of upper mantle seismic anisotropy. Maxima in the magnitude of Te anisotropy are spatially correlated with the locations of most earthquakes and the locations and di- rections of maxima in electrical anisotropy and conductivity. The pattern of brittle failure and seismicity is explained as a response to plate boundary forces acting on a plate of variable rigidity and is expected to induce the observed mechanical anisotropy and enhance the flow of crustal fluids, re- sulting in locally high electrical conductivities. Our combined results thus suggest for the first time that the elastic properties of continental litho- sphere have a leading order influence on the deformation and evolution of convergent plate boundaries. 128 Chapter 7. Conclusions 7.1.3 Teleseismic waveform modelling We investigate the use of a wide-angle, one-way wave equation to model variations in teleseismic 3D waveforms due to 2D elastic heterogeneity and anisotropy. The one-way operators are derived based on a high-frequency approximation of the square-root operator and include the effects of wave propagation as well as multiple scattering. Computational cost is reduced through a number of physically motivated approximations. We present syn- thetic results from simple 1D (layer over a half-space) and 2D (subduction zone) models that are compared with reference solutions. The algorithm is then used to model data from an array of broadband seismograph stations deployed in northwestern Canada as part of the CANOE seismic experiment. In this region radial-component receiver functions show a clear continental Moho and the presence of crustal material dipping into the mantle at the su- ture of two Palaeoproterozoic terranes. The geometry of the suture is better defined on the transverse component where subduction is associated with a ∼10 km thick layer exhibiting strong elastic anisotropy. The modelling reproduces the main features of the receiver functions, including the effects of anisotropy, heterogeneity and finite-frequency scattering. 7.1.4 Slab morphology in northern Cascadia We use passive teleseismic recordings from an array of POLARIS broadband seismic stations to image crustal and upper mantle structure beneath north- ern Vancouver Island into the interior of British Columbia. A clear signature of subducted material extends N-E from Brooks Peninsula at crustal levels, beneath Georgia Strait and the mainland deep into the mantle to 300 km depth. Complexity in slab morphology results from Juan de Fuca ridge subduction and opening of a slab window, in agreement with geophysical and geological data. We propose a tectonic model for the Explorer plate in which its separation from the Juan de Fuca plate is caused by the thermo- mechanical erosion of the slab edge and slab thinning at shallow levels, which slow convergence with NA and lead eventually to plate capture. 129 Chapter 7. Conclusions 7.1.5 Subducted water and megathrust permeability We report observations of anomalously high Poisson’s ratios within the sub- ducted oceanic crust from the northern Cascadia continental margin to its intersection with forearc mantle using converted teleseismic waves. On the basis of pressure, temperature, and compositional considerations, the ele- vated Poissons ratios suggest water is pervasively present in fluid form at pore pressures near lithostatic values. Combined with our observations of a strong negative velocity contrast at the top of the oceanic crust, our results imply that the megathrust is a low-permeability boundary. The transition from a low- to high-permeability plate interface downdip into the mantle wedge is explained by hydrofracturing of the seal by volume changes across the interface caused by the onset of crustal eclogitization and mantle serpen- tinization. These results hold important implications for our understanding of seismogenesis, subduction zone structure, and the mechanism of episodic tremor and slip. 7.2 Limitations The major limitations of the work reported in this thesis reside in the inter- pretation of results which is sometimes purely conceptual (i.e. un-verified), and thus non-unique. For example, Chapters 3, 5 and 6 conclude with an intuitive (qualitative) understanding of tectonic and geodynamic pro- cesses that explain the observations, even though these are obtained using rigourous and thoroughly tested techniques. In the following sections we de- scribe alternative modelling assumptions, untested hypotheses, and incom- plete analysis of the data that could influence the interpretation of results. 7.2.1 Mechanical models of the lithosphere The mechanics of the lithosphere, and especially its rheological structure, are a matter of much debate. For example, geodynamic studies of fault slip rates and lab measurements conclude that the mantle deforms by creep- ing processes and is essentially weak, while geodetic studies of post-glacial 130 Chapter 7. Conclusions rebound and elastic support of plates suggest that the lithospheric mantle contributes significantly to the rigidity of the lithosphere [Bürgmann and Dresen, 2008]. This debate is centered on the applicability of rheological layering mod- els for the upper and lower crust, and upper mantle (the “jelly sandwich”, “crème brûlée”, or “banana split” models [Bürgmann and Dresen, 2008]). These models are based on the assumption that the yield strength envelope (YSE) of the lithosphere, which predicts the amount of stress the lithosphere can take before yielding, can be determined by a combination of brittle defor- mation due to frictional resistance and dislocation creep for each constitutive layer. Brittle failure is a linear function of depth with a slope depending on extensional or compressional regimes and varying with pore-fluid pressure. Power law creep is a non-linear function of temperature, strain rate, and activation energy (i.e. composition) of major constitutive rocks. The lesser of the sustained differential stresses governed by the two regimes defines the YSE. The YSE can be used in conjunction with flexural stress modelling to es- timate the apparent (or effective) elastic thickness of the lithosphere [Lowry and Smith, 1995; Burov and Diament , 1995]. In general the effective Te is strongly dependent on the degree of coupling between adjacent rheological layers, which is subject to considerable uncertainty and unconstrained in most cases. In Chapter 3 we argue that low values of Te in the Cordillera and large Te values in the Craton favor a model where temperature is the major factor influencing rigidity. However an alternative explanation is re- lated to the degree of coupling between the crust and upper mantle, which is also related to temperature via power law creep, but has a more subtle effect on Te. The fact that the crust has low Te (lower than crustal thickness) can be regarded as an indication of crust-mantle decoupling, with the potential result that the mantle could remain relatively strong. In any case we can hardly constrain mantle strength, or assess the validity of any one of the “gourmet” models, using information of Te alone. A second problem of our simplified model of bending-related deformation is that it does not predict finite amplitude strain, which is widely observed 131 Chapter 7. Conclusions in the Cordillera. Resolving this problem requires a thorough 3D modelling study involving the deformation of variable-thickness, visco-elastic, visco- plastic lithosphere, which is a large undertaking and beyond the scope of this work. 7.2.2 Teleseismic attributes of subducted slabs In Chapters 5 and 6 we use a simple stacking technique to characterize dip, depth and VP/VS of an overlying column. However it is apparent that some signals manifest elastic anisotropy and could contaminate the stacks in such a manner as to produce erroneous results. We tried to avoid this uncer- tainty by choosing the signals that we felt confident were not contaminated, for example by selectively stacking phases less affected by anisotropy. How- ever such signals contain useful information regarding directivity in material properties, and can supplement the list of attributes already characterized. Seismic anisotropy can be interpreted in many ways. One of the favored explanation is in terms of rock fabric acquired during finite shear deforma- tion caused by the flow of material. Another cause for anisotropy is the presence of aligned cracks. Both models predict different anisotropic signa- tures which could in theory be distinguished using seismic data. However, our results suggest that fluids, and possibly fluid motions, are important in our context, and seismic anisotropy might reflect the presence of fluid- filled cracks, or favorable conduits for water to percolate upward through the oceanic crust. This hypothesis can be tested by exploring the variablity in Poisson’s ratio observed at different back-azimuthal ranges of incident wavefields. Such measurements are difficult to make, however, and require very clear signals that can be extracted and analyzed at high frequencies (>2 Hz). If the discontinuity in structure is more gradual (e.g. described by a gradient in velocity), it will be transparent to such high frequency signals. Moreover, most stations do not record enough data to make sta- tistically significant inferences on back-azimuthal variability and have poor signal-to-noise ratios. This limitation is currently impeding us from measur- ing Poisson’s ratio along the entire Cascadia margin, even for the isotropic 132 Chapter 7. Conclusions case. It is clear, however, that characterization of the full suite of teleseis- mic attributes, which is currently limited, constitutes the ultimate goal to achieve. 7.3 Future work Here we outline 5 potential research projects that constitute extensions to the present thesis. These are ordered according to their relevance with the material presented in the thesis, and not by priority. 7.3.1 Spherical wavelet transform The directional wavelet coherence method was designed on a Cartesian grid to simplify regional analysis of lithospheric deformation [see Chapter 2; Au- det and Mareschal , 2007]. However, with the increasing availability and resolution of satellite-derived global geophysical data on the Earth (e.g. CHAMP, GRACE), and planets (e.g. Mars Global Surveyor, MESSEN- GER), localized spectral analysis of potential fields on the sphere becomes crucial. Development of the directional, spherical wavelet transform is al- ready under way [McEwen et al., 2006, 2007] and the method could be adapted to the study of planetary potential field with sufficiently resolved gravity and topography. Future work could be directed toward expand- ing the spherical, directional wavelet transform to calculate cross-spectral properties, such as the admittance and coherence. Indeed such comparative planetology is an important way to establish a thorough basic understand- ing of the general problem of the mechanisms governing the deformation of planetary lithospheres. 7.3.2 Deformation at triple junction margins Triple junctions at convergent continental margins offer a remarkable oppor- tunity to refine models of plate deformation because plate boundaries are continuously evolving from convergent to strike-slip relative motion. The western margin of North America illustrates well the complexity of plate 133 Chapter 7. Conclusions interactions at an active margin where the southern edge of the oceanic Gorda plate, which defines the southern boundary of the Cascadia sub- duction zone, is retreating northward with the Mendocino Triple Junction (MTJ), lengthening the strike-slip boundary between the Pacific and North American plates along the San Andreas fault system [Furlong and Schwartz , 2004]. Accretionary faults along the subduction margin represent potential zones of mechanical anisotropy that can evolve into strike-slip faults of the trans- form system after passage of the triple junction. It is reasonable to hypothe- size that the Te pattern and anisotropy might play a major role in continental deformation at a triple junction by concentrating strain along pre-existing zones of weakness. I propose to map Te variations and anisotropy across western United States by calculating the coherence between Bouguer gravity and topography using the directional, wavelet coherence method developed in this thesis [see Chapters 2 & 3; Audet and Mareschal , 2007; Audet et al., 2007b]. This work could also include the investigation of the effects of Te variations and anisotropy using parameterized (2D and ultimately 3D) fi- nite element models with goal of understanding how the different mechanical and geometric boundary conditions implied by these results influence plate boundary dynamics. 7.3.3 Gorda slab morphology To the south of the MTJ the fate of the subducted Gorda slab remains poorly known. The northern boundary of Cascadia is also defined by a triple junction (Queen Charlotte - QCTJ) that translates subduction of the oceanic Explorer plate into a transform boundary [see Chapter 5; Audet et al., 2008b]. As part of the Juan de Fuca subduction system the Gorda and Explorer plates are structurally linked and the fate of the subducted slab at both ends of Cascadia might play a major role in the dynamics of the subduction system (e.g. cessation of subduction and transition into a transform boundary). Results reported in Chapter 5 [see also Audet et al., 2008b] allow us to hypothesize that the seismic signature of the Gorda slab 134 Chapter 7. Conclusions can be traced beyond the triple junction to constrain the lateral and vertical (depth) extent of subduction. I propose to use seismic data from permanent and recently deployed networks and the modelling method described in this thesis [see Chapter 4; Audet et al., 2007a] to image upper mantle structure beneath North America and track the seismic signature of the Gorda slab. 7.3.4 Cascadia forearc structure and ETS Episodic Tremor and Slip (ETS) events in Cascadia and Japan occur in the general vicinity of the plate boundary, downdip of the locked zone. In devel- oping an understanding of the ETS phenomenon it is important to relate the spatial occurence of tremor to the principal structural elements within the subduction complex, in particular the subducted oceanic crust. Through- out Cascadia active and passive seismic studies image a highly-reflective, dipping, low-velocity zone (LVZ) beneath the continental forearc [e.g. see Chapters 5 & 6; Audet et al., 2008a, b; Bostock et al., 2002; Nicholson et al., 2005]. In southern Cascadia this layer is associated with the oceanic crust of the downgoing Juan de Fuca slab, in agreement with expectations from thermal and petrological models [Bostock et al., 2002]. In northern Cas- cadia, where ETS has been extensively documented [Kao et al., 2005], the low-velocity zone is associated with the E-layer, a prominent feature first identified as the subducting oceanic crust [Clowes et al., 1983; Kurtz et al., 1986], but later interpreted to represent trapped fluids in the overriding con- tinental crust [Green et al., 1986; Hyndman, 1988]. The continuity of the LVZ along the margin is currently unconstrained, however, recent evidence suggests its presence beneath northern Vancouver Island, where its signature disappears abruptly [see Chapter 5; Audet et al., 2008b]. An obvious extension to the present thesis is the mapping of the LVZ along the entire Cascadia margin. We have already assembled a large tele- seismic body-wave data set comprising stations from northern California to northern Vancouver Island. Preliminary work suggests that the LVZ is well-developed along the entire convergent margin from the coast eastward to the forearc basins (Georgia Strait, Puget Sound, Willamette Valley). Its 135 Chapter 7. Conclusions diminishing signature further landward and to greater depth is consistent in a petrologic context with an interpretation as subducting oceanic crust. The relation between slab morphology with tremor locations deserves fur- ther investiagtion but previous results [see Chapter 6; Audet et al., 2008a] suggest a link with water release at a permeability transition at the corre- sponding depth where tremors generally occur. A precise correspondence will be established in the near future. 7.3.5 Water transport and seismogenesis Dipping low-velocity zones (LVZs) have been found in several other conti- nental subduction zones, for example in South America [Yuan et al., 2000], Alaska [Ferris et al., 2003], Northeast Japan [Kawakatsu and Watada, 2007], and recently in the Mariana (oceanic) subduction zone [Tibi et al., 2008]. Such evidence suggests LVZs are generic subduction zone features that can be used to map the downgoing plate beneath arcs and forearcs at variable depths. Results of elevated Poisson’s ratio reported in Chapter 6 [Audet et al., 2008a] shows this signature carries information not only on struc- ture but also on hydrology of subducting oceanic crust via high pore-fluid pressure. This observation has numerous implications that need further in- vestigation. The link with ETS has been mentioned in the previous section. A notable effect of high pore-fluid pressure on the seismogenic behavior of faults is to lower the normal effective stress, which allows slip to occur at much lower stresses on “wet” faults than “dry” faults. Knowledge of pore- fluid pressure in subduction zones is thus of utmost importance in studies of triggered earthquakes, fault slip modelling, and earthquake hazard. Unfortunately, a relationship between Poisson’s ratio and pore-fluid pres- sure at ∼1 GPa is lacking at the present time. Although this is an impor- tant area of future research, such measurements are beyond the scope of the present thesis (and beyond the capabilities of the author). What seems reasonable to achieve as an extension of this thesis is a careful mapping of LVZs and Poissons’s ratio in various other subduction zones worldwide. This information can then be linked with the seismogenic behavior of megathrust 136 Chapter 7. Conclusions faults, and from this we can potentially make useful inferences on the seis- mic cycle of megathrust earthquakes. This is very speculative, however, but constitutes an exciting new area of research that can possibly lead to major scientific discoveries. 137 Bibliography Audet, P., and J.-C. Mareschal, Wavelet analysis of the coherence be- tween Bouguer gravity and topography: Application to the elastic thickness anisotropy in the Canadian Shield, Geophys. J. Int., 168, 287–298, 2007. Audet, P., M. G. Bostock, and J.-P. Mercier, Teleseismic waveform mod- elling with a one-way wave equation, Geophys. J. Int., 171, 1212–1225, doi:10.1111/j.1365–246X.2007.03,586.x, 2007a. Audet, P., A. M. Jellinek, and H. Uno, Mechanical controls on the defor- mation of continents at convergent margins, Earth Planet. Sci. Lett., 264, 151–166, doi:10.1016/j.epsl.2007.09.024, 2007b. Audet, P., M. G. Bostock, N. I. Christensen, and S. M. Peacock, Seis- mic evidence for overpressured subducted oceanic crust and sealing of the megathrust, Nature, p. Submitted, 2008a. Audet, P., M. G. Bostock, J.-P. Mercier, and J. F. Cassidy, Morphology of the Explorer/Juan de Fuca slab edge in northern Cascadia: Imaging plate capture at a ridge-trench-transform triple junction, Geology, p. Submitted, 2008b. Bostock, M. G., R. D. Hyndman, S. Rondenay, and S. M. Peacock, An inverted continental Moho and serpentinization of the forearc mantle, Na- ture, 417, 536–539, 2002. Bürgmann, R., and G. Dresen, Rheology of the lower crust and upper mantle: Evidence from rock mechanics, geodesy, and field observations, Annu. Rev. Earth Planet. Sci., 36, 531–567, 2008. 138 Bibliography Burov, E. B., and M. Diament, The effective elastic thickness (Te) of con- tinental lithosphere: What does it really mean?, J. Geophys. Res., 100, 3905–3927, 1995. Clowes, R. M., R. M. Ellis, Z. Hajnal, and I. F. Jones, Seismic reflections from subducting lithosphere?, Nature, 303, 668–670, 1983. Ferris, A., G. A. Abers, D. H. Christensen, and E. Veenstra, High resolution image of the subducted Pacific (?) plate beneath central Alaska, 50150 km depth, Earth Planet. Sci. Lett., 214, 575–588, 2003. Furlong, K. P., and S. Y. Schwartz, Influence of the Mendocino triple junc- tion on the tectonics of coastal California, Annu. Rev. Earth Planet. Sci., 32, 403–433, 2004. Green, A. G., C. J. Y. R. M. Clowes and, C. Spencer, E. R. Kanasewich, M. T. Brandon, and A. S. Brown, Seismic reflection imaging of the sub- ducting Juan de Fuca plate, Nature, 319, 210–213, 1986. Hyndman, R. D., Dipping seismic reflectors, electrically conductive zones, and trapped water in the crust over a subducting plate, J. Geophys. Res., 93, 13,391–13,405, 1988. Kao, H., S.-J. Shan, H. Dragert, G. Rogers, J. F. Cassidy, and K. Ra- machandran, A wide depth distribution of seismic tremors along the north- ern Cascadia margin, Nature, 436, 841–844, 2005. Kawakatsu, H., and S. Watada, Seismic evidence for deep-water trans- portation in the mantle, Science, 316, 1468–1471, 2007. Kurtz, R. D., J. M. DeLaurier, and J. C. Gupta, A magnetotelluric sound- ing across Vancouver Island detects the subducting Juan de Fuca plate, Nature, 321, 596 – 599, 1986. Lowry, A. R., and R. B. Smith, Strength and rheology of the western U. S. Cordillera, J. Geophys. Res., 100, 17,947–17,963, 1995. 139 Bibliography McEwen, J. D., M. P. Hobson, and A. N. Lasenby, A directional continuous wavelet transform on the sphere, ArXiv, p. submitted, 2006. McEwen, J. D., M. P. Hobson, D. J. Mortlock, and A. N. Lasenby, Fast directional continuous spherical wavelet transform algorithms, IEEE Trans. Sig. Proc., 55 (2), 520–529, 2007. Nicholson, T., M. G. Bostock, and J. F. Cassidy, New constraints on sub- duction zone structure in northern Cascadia, Geophys. J. Int., 161, 849– 859, 2005. Tibi, R., D. A. Wiens, and X. Yuan, Seismic evidence for widespread ser- pentinized forearc mantle along the Mariana convergence margin, Geo- phys. Res. Lett., p. in press, 2008. Yuan, X., et al., Subduction and collision processes in the central andes constrained by converted seismic phases, Nature, 408, 958–961, 2000. 140 Master Bibliography Abers, G. A., Seismic low-velocity layer at the top of subducting slabs beneath volcanic arcs: observations, predictions, and systematics, Phys. Earth Planet. Int., 149, 7–29, doi:10.1016/j.pepi.2004.10.002, 2005. Abers, G. A., P. E. van Keken, E. A. Kneller, A. Ferris, and J. C. Stachnik, The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow, Earth Planet. Sci. Lett., 241, 387–397, doi:10.1016/j.epsl.2005.11.055, 2006. Ahrens, T. J., and G. Schubert, Gabbro-eclogite reaction rate and its geo- physical significance, Rev. Geophys., 13, 383–400, 1975. Angus, D. A., A one-way wave equation for modelling seismic waveform variations due to elastic heterogeneity, Geophys. J. Int., 162, 882–898, 2005. Angus, D. A., True amplitude corrections for a narrow-angle one-way elas- tic wave equation, Geophysics, 72, T19–T26, 2007. Angus, D. A., C. J. Thomson, and R. G. Pratt, A one-way wave equation for modelling variations in seismic waveforms due to elastic anisotropy, Geophys. J. Int., 156, 595–614, 2004. Antoine, J.-P., P. Carrette, R. Murenzi, and B. Piette, Image analysis with two-dimensional continuous wavelet transform, Signal Processing, 31, 241– 272, 1993. Antoine, J.-P., R. Murenzi, and P. Vandergheynst, Two-dimensional direc- tional wavelets in image processing, J. Imag. Sys. Technol., 7, 152–165, 1996. 141 Master Bibliography Armstrong, R. L., J. E. Muller, J. E. Harakal, and K. Muehlenbachs, The Neogene Alert Bay Volcanic Belt of northern Vancouver Island, Canada: Descending-plate-edge volcanism in the arc-trench gap, J. Vol- canol. Geotherm. Res., 26, 75–97, 1985. Audet, P., and J.-C. Mareschal, Anisotropy of the flexural response of the lithosphere in the Canadian Shield, Geophys. Res. Lett., 31, L20,601, doi:10.1029/2004GL021,080, 2004a. Audet, P., and J.-C. Mareschal, Variations in elastic thickness in the Cana- dian Shield, Earth Planet. Sci. Lett., 226, 17–31, 2004b. Audet, P., and J.-C. Mareschal, Wavelet analysis of the coherence be- tween Bouguer gravity and topography: Application to the elastic thickness anisotropy in the Canadian Shield, Geophys. J. Int., 168, 287–298, 2007. Audet, P., M. G. Bostock, and J.-P. Mercier, Teleseismic waveform mod- elling with a one-way wave equation, Geophys. J. Int., 171, 1212–1225, doi:10.1111/j.1365–246X.2007.03,586.x, 2007a. Audet, P., A. M. Jellinek, and H. Uno, Mechanical controls on the defor- mation of continents at convergent margins, Earth Planet. Sci. Lett., 264, 151–166, doi:10.1016/j.epsl.2007.09.024, 2007b. Audet, P., M. G. Bostock, N. I. Christensen, and S. M. Peacock, Seis- mic evidence for overpressured subducted oceanic crust and sealing of the megathrust, Nature, p. Submitted, 2008a. Audet, P., M. G. Bostock, J.-P. Mercier, and J. F. Cassidy, Morphology of the Explorer/Juan de Fuca slab edge in northern Cascadia: Imaging plate capture at a ridge-trench-transform triple junction, Geology, accepted, 2008b. Bank, C.-G., M. G. Bostock, R. Ellis, and J. Cassidy, A reconnaissance teleseismic study of the upper mantle and transition zone beneath the Archean Slave craton in NW Canada, Tectonophysics, 319, 151–166, 2000. 142 Master Bibliography Becker, T. W., C. Facenna, R. J. O’Connell, and D. Giardini, The develop- ment of subduction in the upper mantle: insights from experimental and laboratory experiments, J. Geophys. Res., 104, 15,207, 1999. Billen, M. I., M. Gurnis, and M. Simons, Multiscale dynamic models of the Tonga-Kermadec subduction zone, Geophys. J. Int., 153, 359–388, 2003. Boerner, D. E., R. D. Kurtz, J. A. Craven, G. M. Ross, and F. W. Jones, A synthesis of electromagnetic studies in the Lithoprobe Alberta Basement Transect: constraints on Paleoproterozoic indentation tectonics, Can. J. Earth Sci., 37, 1509–1534, 1999. Bostock, M. G., Mantle stratigraphy and evolution of the Slave province, J. Geophys. Res., 103, 21,183–21,200, 1998. Bostock, M. G., Regional methods, in Seismology and structure of the Earth, Treatise on Geophysics, vol. 1, edited by G. Schubert, pp. 73–78, Elsevier, Amsterdam, The Netherlands, 2007. Bostock, M. G., and J. C. VanDecar, Upper mantle structure of the north- ern Cascadia subduction zone, Can. J. Earth Sci., 32, 1–12, 1995. Bostock, M. G., R. D. Hyndman, S. Rondenay, and S. M. Peacock, An inverted continental Moho and serpentinization of the forearc mantle, Na- ture, 417, 536–539, 2002. Braunmiller, J., and J. Nábelek, Seismotectonics of the Explorer region, J. Geophys. Res., 107, 2208, doi:10.1029/2001JB000,220, 2002. Bürgmann, R., and G. Dresen, Rheology of the lower crust and upper mantle: Evidence from rock mechanics, geodesy, and field observations, Annu. Rev. Earth Planet. Sci., 36, 531–567, 2008. Burov, E. B., and M. Diament, The effective elastic thickness (Te) of con- tinental lithosphere: What does it really mean?, J. Geophys. Res., 100, 3905–3927, 1995. 143 Master Bibliography Burov, E. B., and A. B. Watts, The long-term strength of continental lithosphere: “jelly sandwich” or “crème brûlée”?, GSA Today, 16, 4–10, 2006. Caine, J. S., J. P. Evans, and C. B. Foster, Fault zone architecture and permeability structure, Geology, 24, 1025–1028, 1996. Carr, S. D., R. Parrish, and D. L. Parkinson, Eocene extensional tectonics and geochronology of the southern Omineca Belt, British Columbia and Washington, Tectonics, 7, 181–212, 1988. Cassidy, J. F., R. M. Ellis, C. Karavas, and G. C. Rogers, The northern limit of the subducted Juan de Fuca plate system, J. Geophys. Res., 103, 26,949–26,961, 1998. Chapman, C., Fundamentals of Seismic Wave Propagation, Cambridge University Press, Cambridge, UK, 2004. Christensen, N. I., Pore pressure and oceanic crustal seismic structure, Geophys. J. R. astr. Soc., 79, 411–423, 1984. Christensen, N. I., Poisson’s ratio and crustal seismology, J. Geophys. Res., 101, 3139–3156, 1996. Christensen, N. I., Serpentinites, Peridotites, and Seismology, Int. Geol. Rev., 46, 795–816, 2004. Claerbout, J., Imaging the Earth’s interior, Blackwell Sci. Publ., Oxford, United Kingdom (GBR), 1985. Clayton, R., and B. Engquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bull. Seism. Soc. Am., 67, 1529–1540, 1977. Clowes, R. M., R. M. Ellis, Z. Hajnal, and I. F. Jones, Seismic reflections from subducting lithosphere?, Nature, 303, 668–670, 1983. Clowes, R. M., D. J. Baird, and S. A. Dehler, Crustal structure of the Cascadia subduction zone, southwestern British Columbia, from potential field and seismic studies, Can. J. Earth Sci., 34, 317–335, 1997. 144 Master Bibliography Coleman, R. G., Petrologic and geophysical nature of serpentinites, Geol. Soc. Am. Bull., 82, 897–918, 1971. Cook, F. A., and P. Erdmer, An 1800 km cross section of the lithosphere through the northwestern North American plate: lessons from 4.0 billion years of Earth’s history, Can. J. Earth Sci., 42, 1295–1311, 2005. Currie, C. A., and R. D. Hyndman, The thermal structure of subduction zone backarcs, J. Geophys. Res., 111, B08,404, doi:10.1029/2005JB004,024, 2006. Currie, C. A., J. F. Cassidy, and R. D. Hyndman, A regional study of shear wave splitting above the Cascadia subduction zone: Margin-parallel crustal stress, Geophys. Res. Lett., 28, 659–662, 2001. Currie, C. A., J. F. Cassidy, R. D. Hyndman, and M. G. Bostock, Shear wave anisotropy beneath the Cascadia subduction zone and western North American craton, Geophys. J. Int., 157, 341–353, 2004. DeMets, C., and S. Traylen, Motion of the Rivera plate since 10 Ma relative to the Pacific and North American plates and the mantle, Tectonophysics, 318, 119–159, 2000. Dziak, R. P., Explorer deformation zone: Evidence of a large shear zone and reorganization of the Pacific-Juan de Fuca-North American triple junction, Geology, 34, 213–216, doi:10.1130/G22,164.1, 2006. Eaton, D., A. Jones, and I. Ferguson, Lithospheric anisotropy structure in- ferred from collocated teleseismic and magnetotelluric observations: Great Slave Lake shear zone, northern Canada, Geophys. Res. Lett., 31, L19,614, doi:10.1029/2004GL020,939, 2004. Eaton, D. W., and A. G. Jones, Tectonic fabric of the subcontinental litho- sphere: Evidence from seismic, magnetotelluric and mechanical anisotropy, Phys. Earth Planet. Int., 158, 85–91, 2006. Efron, B., and R. Tibshirani, Statistical data analysis in the computer age, Science, 253, 390–395, 1991. 145 Master Bibliography Ferris, A., G. A. Abers, D. H. Christensen, and E. Veenstra, High resolution image of the subducted Pacific (?) plate beneath central Alaska, 50150 km depth, Earth Planet. Sci. Lett., 214, 575–588, 2003. Fisher, A. T., Permeability within basaltic oceanic crust, Rev. Geophys., 36, 143–182, 1998. Flück, P., R. D. Hyndman, and C. Lowe, Effective elastic thickness Te of the lithosphere in western Canada, J. Geophys. Res., 108, 2430, doi:10.1029/2002JB002,201, 2003. Forsyth, D., and S. Uyeda, On the relative importance of the driving forces of plate motion, Geophys. J. R. Astron. Soc., 43, 163–200, 1975. Forsyth, D. W., Subsurface loading and estimates of the flexural rigidity of continental lithosphere, J. Geophys. Res., 90, 12,623–12,632, 1985. Foufoula-Georgiou, E., and P. Kumar (Eds.),Wavelets in Geophysics, Aca- demic Press, San Diego, Calif., 1994. Frederiksen, A. W., and M. G. Bostock, Modelling teleseismic waves in dipping anisotropic structures, Geophys. J. Int., 141, 401–412, 2000. Frederiksen, A. W., M. G. Bostock, and J. F. Cassidy, S-wave velocity structure of the Canadian upper mantle, Phys. Earth Planet. Inter., 124, 175–191, 2001. Furlong, K. P., and S. Y. Schwartz, Influence of the Mendocino triple junc- tion on the tectonics of coastal California, Annu. Rev. Earth Planet. Sci., 32, 403–433, 2004. Gaspar-Escribano, J. M., J. D. van Wees, M. ter Voorde, S. Cloetingh, E. Roca, L. Cabrera, J. A. Munoz, P. A. Ziegler, and D. Garcia-Castellanos, Three-dimensional flexural modelling of the Ebro Basin (NE Iberia), Geo- phys. J. Int., 145, 349–367, 2001. 146 Master Bibliography Gerver, M., and V. Markushevitch, Determination of a seismic wave veloc- ity from the travel time curve, Geophys. J. Roy. astron. Soc., 11, 165–173, 1966. Govers, R., and P. T. Meijer, On the dynamics of the Juan de Fuca plate, Earth Planet. Sci. Lett., 189, 115–131, 2001. Green, A. G., C. J. Y. R. M. Clowes and, C. Spencer, E. R. Kanasewich, M. T. Brandon, and A. S. Brown, Seismic reflection imaging of the sub- ducting Juan de Fuca plate, Nature, 319, 210–213, 1986. Grimbergen, J. L. T., F. J. Dessing, and K. Wapenaar, Modal expansion of one-way operators in laterally varying media, Geophysics, 63, 995–1005, 1998. Holschneider, M., Wavelets: An Analysis Tool, Oxford University Press, Oxford, 1999. Hoop, M. V. D., J. H. L. Rousseau, and R.-S. Wu, Generalization of the phase-screen approximation for the scattering of acoustic waves, Wave Mo- tion, 31, 43–70, 2000. Hyndman, R. D., Dipping seismic reflectors, electrically conductive zones, and trapped water in the crust over a subducting plate, J. Geophys. Res., 93, 13,391–13,405, 1988. Hyndman, R. D., and T. J. Lewis, Geophysical consequences of the Cordillera-Craton thermal transition in southwestern Canada, Tectono- physics, 306, 397–422, 1999. Hyndman, R. D., and S. M. Peacock, Serpentinization of the forearc mantle, Earth Planet. Sci. Lett., 212, 417–432, doi:10.1016/S0012–821X(03)00,263– 2, 2003. Hyndman, R. D., P. Flueck, S. Mazzotti, T. J. Lewis, J. Ristau, and L. Leonard, Current tectonics of the northern Canadian Cordillera, Can. J. Earth Sci., 42, 1117–1136, 2005. 147 Master Bibliography Jackson, J., Strength of the continental lithosphere, GSA Today, 12, 4–10, 2002. Jones, A. G., Electromagnetic interrogation of the anisotropic Earth: Look- ing into the Earth with polarized spectacles, Phys. Earth Planet. Int., 158, 281–291, 2006. Jordan, T. H., Composition and development of the continental tecto- sphere, Nature, 274, 544–548, 1978. Kao, H., S.-J. Shan, H. Dragert, G. Rogers, J. F. Cassidy, and K. Ra- machandran, A wide depth distribution of seismic tremors along the north- ern Cascadia margin, Nature, 436, 841–844, 2005. Karner, G. D., M. S. Steckler, and J. A. Thorne, Long-term thermo- mechanical properties of the continental lithosphere, Nature, 304, 250–253, 1983. Kato, A., A. Sakaguchi, S. Yoshida, and H. Mochizuki, Permeability mea- surements and precipitation sealing of basalt in an ancient exhumed fault of a subduction zone, Bull. Earthquake Res. Inst., Univ. of Tokyo, 78, 83–89, 2003. Kawakatsu, H., and S. Watada, Seismic evidence for deep-water trans- portation in the mantle, Science, 316, 1468–1471, 2007. Kennett, B. L., Seismic wave propagation in stratified media, Cambridge University Press, Cambridge, UK, 1983. Kido, M., D. A. Yuen, and A. P. Vincent, Continuous wavelet-like filter for a spherical surface and its application to localized admittance function on Mars, Phys. Earth Planet. Inter., 135, 1–14, 2003. Kirby, J. F., Which wavelet best reproduces the Fourier power spectrum?, Comput. Geosc, 31, 846–864, 2005. 148 Master Bibliography Kirby, J. F., and C. J. Swain, Global and local isostatic coherence from the wavelet transform, Geophys. Res. Lett., 31, doi:10.1029/2004GL021,569, 2004. Kirby, J. F., and C. J. Swain, Mapping the mechanical anisotropy of the lithosphere using a 2-D wavelet coherence, and its application to Australia, Phys. Earth Planet. Int., 158, 122–138, 2006. Kirby, S. H., E. R. Engdahl, and R. Denlinger, Intermediate-depth in- traslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs (overview), in Geophys. Monogr. 96, edited by G. E. Bebout, D. W. Scholl, S. H. Kirby, and J. P. Platt, pp. 195–214, American Geophysical Union, Washington, DC, 1996. Kodaira, S., T. Iidaka, A. Kato, J.-O. Park, T. Iwasaki, and Y. Kaneda, High pore fluid pressure may cause silent slip in the Nankai Trough, Sci- ence, 304, 1295–1298, 2004. Kohlstedt, D. L., J. B. Evans, and S. J. Mackwell, Strength of the litho- sphere: Constraints imposed by laboratory experiments, J. Geophys. Res., 100, 17,587–17,602, 1995. Komatitsch, D., and J. Tromp, Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophys. J. Int., 139, 806– 822, 1999. Kosarev, G., R. Kind, S. V. Sobolev, X. Yuan, W. Hanka, and S. Oreshin, Seismic evidence for a detached Indian lithospheric mantle beneath Tibet, Science, 283, 1306–1309, 1999. Kosloff, D., M. Reshef, and D. loewenthal, Elastic wave calculations by the fourier method, Bull. Seis. Soc. Am., 74, 875–891, 1984. Kreemer, C., R. Govers, K. P. Furlong, and W. E. Holt, Plate boundary deformation between the Pacific and North America in the Explorer region, Tectonophysics, 293, 225–238, 1998. 149 Master Bibliography Kuhn, T. S., The structure of scientific revolutions, University of Chicago Press, Chicago, USA, 1962. Kumar, P., A wavelet-based methodology for scale-space anisotropic anal- ysis, Geophys. Res. Lett., 22, 2777–2780, 1995. Kurtz, R. D., J. M. DeLaurier, and J. C. Gupta, A magnetotelluric sound- ing across Vancouver Island detects the subducting Juan de Fuca plate, Nature, 321, 596 – 599, 1986. Ledo, J., and A. G. Jones, Regional electrical resistivity structure of the southern Canadian Cordillera and its physical interpretation, J. Geo- phys. Res., 106, 30,755–30,769, 2001. Levander, A., Fourth-order finite-difference P-SV seismograms, Geophysics, 53, 1425–1436, 1988. Levin, V., D. Okaya, and J. Park, Shear wave birefringence in wedge-shaped anisotropic regions, Geophys. J. Int., 168, 275–286, doi: 10.1111/j.1365– 246X.2006.03,224.x, 2007. Lewis, T. J., C. Lowe, and T. S. Hamilton, Continental signature of a ridge-trench-triple junction: Northern Vancouver Island, J. Geophys. Res., 102, 7,767–7,781, 1997. Lowry, A. R., and R. B. Smith, Flexural rigidity of the Basin and Range- Colorado Plateau-Rocky Mountain transition from coherence analysis of gravity and topography, J. Geophys. Res., 99, 20,123–20,140, 1994. Lowry, A. R., and R. B. Smith, Strength and rheology of the western U. S. Cordillera, J. Geophys. Res., 100, 17,947–17,963, 1995. Macario, A., A. Malinverno, and W. F. Haxby, On the robustness of elastic thickness estimates obtained using the coherence method, J. Geophys. Res., 100, 15,163–15,172, 1995. 150 Master Bibliography Maggi, A., J. A. Jackson, D. McKenzie, and K. Priestley, Earthquake fo- cal depths, effective elastic thickness, and the strength of the continental lithosphere, Geology, 28, 495–498, 2000. Mareschal, M., R. L. Kellett, R. D. Kurtz, J. N. Ludden, S. Ji, and R. C. Bailey, Archaean cratonic roots, mantle shear zones and deep electrical anisotropy, Nature, 375, 134–137, 1995. Mazzotti, S., H. Dragert, J. Henton, M. Schmidt, R. Hyndman, T. James, Y. Lu, and M. Craymer, Current tectonics of northern Cas- cadia from a decade of GPS measurements, J. Geophys. Res., 108, 2554, doi:10.1029/2003JB002,653, 2003. McCoy, J., and L. N. Frazer, Propagation modelling based on wavefield factorization and invariant embedding, Geophys. J. R. astr. Soc., 86, 703– 717, 1986. McCoy, J., L. Fishman, and L. N. Frazer, Reflection and transmission at an interface separating transversely inhomogeneous acoustic half-spaces, Geophys. J. R. astr. Soc., 85, 543–562, 1986. McEwen, J. D., M. P. Hobson, and A. N. Lasenby, A directional continuous wavelet transform on the sphere, ArXiv, p. submitted, 2006. McEwen, J. D., M. P. Hobson, D. J. Mortlock, and A. N. Lasenby, Fast directional continuous spherical wavelet transform algorithms, IEEE Trans. Sig. Proc., 55 (2), 520–529, 2007. McNeill, A. F., Contribution of post-critical reflections to ground motions from mega-thrust events in the Cascadia subduction zone, Master’s thesis, University of Victoria, BC, Canada, 2005. Meneghini, F., and J. C. Moore, Deformationand hydrofracturing in a sub- duction thrust at seismogenic depths: The Rodeo Cove thrust zone, Marin Headlands, California, Geol. Soc. Am. Bull., 119, 174–183, 2007. 151 Master Bibliography Mercier, J.-P., M. G. Bostock, P. Audet, J. B. Gaherty, E. J. Garnero, and J. Revenaugh, The teleseismic signature of fossil subduction: Northwestern Canada, J. Geophys. Res., 113, B04,308, doi:10.1029/2007JB005,127, 2008. Moore, J. C., and P. Vrolijk, Fluids in accretionary prisms, Rev. Geophys., 30, 113–135, 1992. Nicholson, T., M. G. Bostock, and J. F. Cassidy, New constraints on sub- duction zone structure in northern Cascadia, Geophys. J. Int., 161, 849– 859, 2005. Ojeda, G. Y., and D. Whitman, Effect of windowing on lithosphere elas- tic thickness estimates obatined via the coherence method: Results from northern South America, JGR, 107, 2275, doi:10.1029/2000JB000,114, 2002. Park, J., and V. Levin, Seismic anisotropy: Tracing plate dynamics in the mantle, Science, 296, 485–489, 2002. Peacock, S. M., Fluid processes in subduction zones, Science, 248, 329–337, 1990. Pérez-Gussinyé, M., and A. B. Watts, The long-term strength of Europe and its implications for plate- forming processes, Nature, 436, 381–384, 2005. Pérez-Gussinyé, M., A. Lowry, A. B. Watts, and I. Velicogna, On the recovery of effective elastic thickness using spectral methods: Examples from synthetic data and from the Fennoscandian Shield, J. Geophys. Res., 109, B10,409, doi:10.1029/2003JB002,788, 2004. Pérez-Gussinyé, M., A. R. Lowry, and A. B. Watts, Effective elastic thick- ness of South America and its implications for intracontinental deforma- tion, Geochem. Geophys. Geosys., 5, Q05,009, doi:10.1029/2006GC001,511, 2007. 152 Master Bibliography Perry, H. K. C., D. W. S. Eaton, and A. M. Forte, LITH5.0: a revised crustal model for Canada based on Lithoprobe results, Geophys. J. Int., 150, 285–294, 2002. Poppeliers, C., and G. Pavlis, Three-dimensional, prestack, plane wave mi- gration of teleseismic P -to-S converted phases: 2. Stacking multiple events, J. Geophys. Res., 108, 2267, doi:10.1029/2001JB001,583, 2003. Poudjom Djomani, Y. H., J. D. Fairhead, and W. L. Griffin, The flexu- ral rigidity of Fennoscandia: Reflection of the tectonothermal age of the lithospheric mantle, Earth Planet. Sci. Lett., 174, 139–154, 1999. Poudjom-Djomani, Y. H., S. Y. O’Reilly, W. L. Griffin, L. M. Nat- apov, Y. Erinchek, and J. Hronsky, Upper mantle structure beneath eastern Siberia: Evidence from gravity modelling and mantle petrology, Geochem. Geophys. Geosys., 4, 1066, doi:10.1029/2002GC000,420., 2003. Price, R. A., The southeastern Canadian Cordillera: Thrust faulting, tec- tonic wedging, and delamination of the lithosphere, J. Struct. Geol., 8, 239–254, 1986. Rajesh, R. S., and D. C. Mishra, Lithospheric thickness and mechanical strength of the indian shield, Earth Planet. Sci. Lett., 225, 319–328, 2004. Rajesh, R. S., J. Stephen, and D. C. Mishra, Isostatic response and anisotropy of the Eastsern Himalayan-Tibetan Plateau: A reap- praisal using multitaper spectral analysis, Geophys. Res. Lett., 30, 1060, doi:10.1029/2002GL016,104, 2003. Ramachandran, K., S. E. Dosso, G. D. Spence, R. D. Hyndman, and T. M. Brocher, Forearc structure beneath southwestern British Columbia: A three-dimensional tomographic velocity model, J. Geophys. Res., 110, B02,303, doi:10.1029/2004JB003,258, 2005. Ramachandran, K., R. D. Hyndman, and T. M. Brocher, Regional p wave velocity structure of the Northern Cascadia Subduction Zone, J. Geo- phys. Res., 111, B12,301, doi:10.1029/2005JB004,108, 2006. 153 Master Bibliography Ranalli, G., Rheology of the Earth, 2nd ed., Chapman and Hall, London, UK, 1995. Rohr, K. M. M., and K. P. Furlong, Ephemeral plate tectonics at the Queen Charlotte triple junction, Geology, 23, 1035–1038, 1995. Rondenay, S., G. A. Abers, and P. E. van Keken, Seismic imaging of subduction zone metamorphism, Geology, 36, 275–278, doi:10.1130/G24,112A.1, 2008. Rossi, G., G. A. Abers, S. Rondenay, and D. H. Christensen, Unusual mantle Poisson’s ratio, subduction, and crustal structure in central Alaska, J. Geophys. Res., 111, B09,311, doi:10.1029/2005JB003,956, 2006. Savage, M. K., Seismic anisotropy and mantle deformation: What have we learned from shear wave splitting?, Rev. Geophys., 37, 65–106, 1999. Schimmel, M., and H. Paulssen, Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int., 130, 497– 505, 1997. Shelly, D. R., G. C. Beroza, and S. Ide, Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip, Nature, 442, 488–491, 2006. Shragge, J., M. G. Bostock, C. G. Bank, and R. M. Ellis, Integrated tele- seismic studies of the southern Alberta upper mantle, Can. J. Earth Sci., 39, 399–411, 2002. Silver, P. G., Seismic anisotropy beneath the continents: Probing the depths of geology, Annu. Rev. Earth Planet. Sci., 24, 385–432, 1996. Simons, F. J., and R. D. van der Hilst, Age-dependent seismic thickness and mechanical strength of the Australian lithosphere, Geophys. Res. Lett., 29, 1529, doi:10.1029/2002GL014,962, 2002. Simons, F. J., and R. D. van der Hilst, Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere, 154 Master Bibliography Earth Planet. Sci. Lett., 211, 271–286, doi:10.1016/S0012–821X(03)00,198– 5, 2003. Simons, F. J., M. T. Zuber, and J. Korenaga, Isostatic response of the Aus- tralian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper spectral analysis, J. Geophys. Res., 105, 19,163–19,184, 2000. Simons, F. J., R. D. van der Hilst, and M. T. Zuber, Spatio-spectral lo- calization of isostatic coherence anisotropy in Australia and its relation to seismic anisotropy: Implications for lithospheric deformation, J. Geo- phys. Res., 108, 2250, doi:10.1029/2001JB000,704, 2003. Simpson, F., Resistance to mantle flow inferred from the electromagnetic strike of the Australian upper mantle, Nature, 412, 632–635, 2001. Soyer, W., and M. Unsworth, Deep electrical structure of the northern Cascadia (British Columbia, Canada) subduction zone: Implications for the distribution of fluids, Geology, 34, 53–56, doi:10.1130/G21,951.1, 2006. Stark, C. P., J. Stewart, and C. J. Ebinger, Wavelet transform mapping of effective elastic thickness and plate loading: Validation using synthetic data and application to the study of southern African tectonics, J. Geo- phys. Res., 108, 2258, doi:10.1029/2001JB000,609, 2003. Stegman, D. R., J. Freeman, W. P. Schellart, L. Moresi, and D. May, Influence of trench width on subduction hinge retreat rates in 3-D models of slab rollback, Geochem. Geophys. Geosys, 7, Q03,012, 2006. Stern, R. J., Subduction zones, Rev. Geophys., 40, 1012, doi:10.1029/2001RG000,108, 2002. Stewart, J., and A. B. Watts, Gravity anomalies and spatial variations of flexural rigidity at mountain ranges, J. Geophys. Res., 102, 5327–5352, 1997. Stock, J. M., and J. Lee, Do microplates in a subduction zones leave a geological record?, Tectonics, 13, 1471–1487, 1994. 155 Master Bibliography Swain, C. J., and J. F. Kirby, The effect of ”noise” on estimates of the elastic thickness of the continental lithopshere by the coherence method, Geophys. Res. Lett., 30, 1574, doi:10.1029/2003GL017,070, 2003a. Swain, C. J., and J. F. Kirby, The coherence method using a thin anisotropic plate model, Geophys. Res. Lett., 30, 2014, doi:10.1029/2003GL018,350, 2003b. Swain, C. J., and J. F. Kirby, An effective elastic thickness map of Australia from wavelet transforms of gravity and topography using Forsyth’s method, Geophys. Res. Lett., 33, L02,314, doi:10.1029/2005GL025,090, 2006. Tassara, A., C. Swain, R. Hackney, and J. Kirby, Elastic thickness structure of South America estimated using wavelets and satellite-derived gravity data, Earth Planet. Sci. Lett., 253, 17–36, 2007. ten Brink, U., I. Shimizu, and P. Molzer, Plate deformation at depth under northern California: Slab gap or stretched slab?, Tectonics, 18, 1084–1098, 1999. Thomas, W. A., Tectonic inheritance at a continental margin, GSA Today, 16, doi:10.1130/1052–5173, 2006. Thomson, C. J., Notes on waves in layered media to accompany program Rmatrix, in Seismic waves in complex 3-D structures, pp. 147–162, De- partment of Geophysics, Charles University, 1996. Thomson, C. J., The ’gap’ between seismic ray theory and ’full’ wavefield extrapolation, Geophys. J. Int., 137, 364–380, 1999. Thomson, C. J., Accuracy and efficiency considerations for wide-angle wavefield extrapolators and scattering operators, Geophys. J. Int., 163, 308–323, 2005. Thomson, D. J., Spectrum estimation and harmonic analysis, Proc. IEEE, 70, 1055–1096, 1982. 156 Master Bibliography Thorkelson, D., and K. Breisprecher, Partial melting of slab window mar- gins: Genesis of adakitic and non-adakitic magmas, Lithos, 79, 25–41, 2005. Tibi, R., D. A. Wiens, and X. Yuan, Seismic evidence for widespread ser- pentinized forearc mantle along the Mariana convergence margin, Geo- phys. Res. Lett., p. in press, 2008. Torrence, C., and G. P. Compo, A practical guide to wavelet analysis, Bull. Am. Meter. Soc., 79, 61–78, 1998. Turcotte, D. L., and G. Schubert, Geodynamics, Cambridge Univ. Press, Cambridge, UK, 2002. Vinnik, L. P., L. I. Makeyeva, A. Milev, and A. Y. Usenko, Global pat- terns of azimuthal anisotropy and deformations in the continental mantle, Geophys. J. Int., 111, 433–447, 1992. Vinnik, L. P., R. W. E. Green, and L. O. Nicolaysen, Recent deformations of the deep continental root beneath southern Africa, Nature, 375, 50–52, 1995. Virieux, J., P-SV wave propagation in heterogeneous media; velocity-stress finite-difference method, Geophysics, 51, 889–901, 1986. Wang, K., Simplified analysis of horizontal stresses in a buttressed forearc sliver at an oblique subduction zone, Geophys. Res. Lett., 23, 2021–2024, 1996. Wang, Z., D. Zhao, O. P. Mishra, and A. Yamada, Structural heterogeneity and its implications for the low frequency tremors in Southwest Japan, Earth Planet. Sci. Lett., 251, 66–78, doi:10.1016/j.epsl.2006.08.025, 2006. Watts, A. B., Isostasy and Flexure of the Lithosphere, Cambridge Univ. Press, Cambridge, UK, 2001. Watts, A. B., and E. B. Burov, Lithospheric strength and its relationship to the elastic and seismogenic layer thickness, Earth Planet. Sci. Letters, 213, 113–131, 2003. 157 Master Bibliography Wielandt, E., On the validity of the ray approximation for interpreting delay times, in Seismic tomography, edited by G. Nolet, pp. 85–98, Reidel, 1987. Wild, A. J., and J. A. Hudson, A geometrical approach to the elastic complex screen, J. Geophys. Res., 103, 707–725, 1998. Wu, R.-S., Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, J. Geophys. Res., 99, 751–766, 1994. Yuan, X., et al., Subduction and collision processes in the central andes constrained by converted seismic phases, Nature, 408, 958–961, 2000. Zandt, G., and C. J. Ammon, Continental crust composition constrained by measurements of crustal Poisson’s ratio, Nature, 374, 152–154, 1995. Zandt, G., H. Gilbert, T. J. Owens, M. Ducea, J. Saleeby, and C. H. Jones, Active foundering of a continental arc root beneath the southern Sierra Nevada in California, Nature, 431, 41–46, 2004. Zhang, X., H. Paulssen, S. Lebedev, and T. Meier, Surface wave to- mography of the Gulf of California, Geophys. Res. Lett., 34, L15,305, doi:10.1029/2007GL030,631, 2007. Zhu, L., and H. Kanamori, Moho depth variation in southern California from teleseismic receiver functions, J. Geophys. Res., 105, 2969–2980, 2000. 158 Appendix A Approximate amplitude normalizing factor In this section we derive a first order expression for the differential trans- mission operator of forward going waves, following the same approach as Thomson [2005]. The derivation is based on the propagation through a ver- tically layered medium, and the expressions are extrapolated to the case of a laterally varying medium. The anisotropic elastic wave equation at a single frequency ω in the absence of source is written ∂j(cijkl∂kul) + ω 2ρui = 0. (A.1) The elastic tensor cijkl can be written in a more compact form via (Cjk)il = cijkl. In a laterally homogeneous medium, equation (A.1) can be simplified to give C11∂ 2 1u+(Cα1+C1α)∂α∂1u+∂1C11∂1u+∂1C1α∂αu+Cαβ∂α∂βu+ω 2ρu = 0. (A.2) This is equation 2.11 in Thomson [1999] with the terms consisting of lateral derivatives of the elasticity tensor removed. By multiplying on the left by C−111 we get ∂21u+C −1 11 (Cα1 +C1α)∂α∂1u+C −1 11 Cαβ∂α∂βu+C −1 11 ∂1C11∂1u+ C−111 ∂1C1α∂αu+ ω 2ρC−111 u = 0. (A.3) 159 Appendix A. Approximate amplitude normalizing factor This equation may be exactly rewritten (∂1+iD+E)(∂1−iD)u+i(∂1D)u−D2u+iEDu+Mu = ∂21u+E∂1u+Mu = 0, where M = C−111 (ω 2ρI+ ∂1C1α∂α +Cαβ∂α∂β), E = C−111 ((Cα1 +C1α)∂α + ∂1C11), and D(x1, ∂α) is an unknown PDO. If D can be found such that ( D2 − i(∂1D)− iED ) u = Mu, (A.4) then (∂1 − iD)u = 0 is an exact one-way equation for waves propagating in the forward direction. The Fourier representation of the PDOs E and M are simply M = ω2C−111 (ρI−Cαβpαpβ) + iωC−111 ∂1C1αpα, E = iωC−111 (Cα1 +C1α) pα +C −1 11 ∂1C11. The square-root D is required in the factorization (A.4). Recalling that when inhomogeneities exist the root operator can only be found from a high-frequency approximation Thomson [1999], we examine an asymptotic expansion of the operator D in the form D = ωD0 +D1 +O( 1 ω ). (A.5) We then calculate the root operator from (A.4) by substituting the asymp- totic approximation (A.5) and matching terms with equal order in ω. We get for the highest order in ω 0 = ( ω2D20 + ω 2C−111 (Cα1 +C1α)pαD0 − ω2C−111 (ρI−Cαβpαpβ) ) u. (A.6) We introduce the Christoffel equation to help explain the derivation of the 160 Appendix A. Approximate amplitude normalizing factor factorization as (cijklpjpk − ρδil)gl = 0, (A.7) which can be written in the matrix form as [ p21I+C −1 11 (Cα1 +C1α)pαp1 +C −1 11 Cαβpαpβ − ρC−111 ] g = 0. (A.8) Recall that the Christoffel equation relates allowable polarizations g and slownesses p for plane waves of the form: u(x) = g(m) exp(iωP (m) 1 x1) exp(iωpαxα), (A.9) where m indicates the polarization type (P , S1, S2). Introducing the 3× 3 matrix G with columns given by the 3 eigen-displacements of the forward- going waves, and P1 the diagonal matrix of corresponding eigenvalues p1, equation (A.8) becomes the augmented Christoffel equation [ GP21 +C −1 11 (Cα1 +C1α)pαGP1 +C −1 11 [Cαβpαpβ − ρI]G ] = 0. (A.10) We now have solutions for u of the form u = Gv, (A.11) where v are the upgoing wave vector coefficients. Inserting this expression into equation (A.8) and comparing with the augmented Christoffel equation (A.10) we obtain D0 = GP1G −1, (A.12) as expected. This leading order term properly accounts for the phase propa- gation and the coupling between wave modes Thomson [1999]. Returning to equation (A.4), retaining the terms with coefficient ω1, and grouping terms depending on D1 on the left-hand side yields (D0 + 2A)D1 +D1D0 = iH, (A.13) 161 Appendix A. Approximate amplitude normalizing factor where the right-hand side is written H = ∂1D0 +N0D0 +Nα, (A.14) and the matrices A,N0 and Nα are given by A = 1/2C−111 (Cα1 +C1α)pα, N0 = C −1 11 ∂1C11, Nα = C −1 11 ∂1C1αpα. Equation (A.13) is equivalent to equation (B2) of Thomson [1999], al- though the defining matrices are slightly different. An intuitive understand- ing of the term D1 can be found by restricting the propagation to normal incidence. In this case, pα = 0, and equation (A.13) reduces to iC−111 ∂1 ( C11D ′ 0 ) = D′0D ′ 1 +D ′ 1D ′ 0, (A.15) where the primes denote vertical incidence. The propagator matrix D′0 is diagonal and we can write D′1 = i 2 ( C11D ′ 0 )−1 ∂1 ( C11D ′ 0 ) . (A.16) This is the energy-flux normalization term, which includes the effects of elasticity and density gradients. Assuming the vertical gradient of C11 is negligible compared to that of D0, this term can be rewritten D′1 = i 2 G′ ( P′−11 ∂1P ′ 1 ) G ′−1, (A.17) which is the differential transmission coefficient analogous to the acoustic wide-angle wave equation (14) of Thomson [2005]. This term is in fact the zeroth order term in a Taylor expansion about pα, i.e. the narrow-angle approximation. Higher order terms describe the amplitude correction asso- ciated with wave coupling from lateral gradients and wavefront curvature. 162 Appendix A. Approximate amplitude normalizing factor In our teleseismic application, we can reasonably assume that wavefront curvature is negligible. As stated by Thomson [1999], a simple explicit form for the solution of (A.13) appears to be unattainable and the inversion of a 9 × 9 matrix is required to obtain D1. This approach, although more exact, can significantly increase the computational labor of the wide-angle one-way method. An approximate solution may be sought through the narrow-angle expansion of the propagator as discussed above. This approach has been treated by Angus [2007]. We provide an alternate derivation of the energy-flux normalization term D1: D1 ≈ i 2 G ( (C11P1) −1 ∂1 (C11P1) ) G−1, (A.18) This form is a crude approximation of the exact form (A.13) and neglects possible effects generated from the coupling of waves at non-vertical inci- dence in isotropic media. However we adopt this form because it is easily implemented in the wide-angle formulation and the associated inaccuracies, observed from the comparison of numerical results with a reference solution, are small, as demonstrated by the numerical examples. 163 Bibliography Angus, D. A., True amplitude corrections for a narrow-angle one-way elas- tic wave equation, Geophysics, 72, T19–T26, 2007. Thomson, C. J., The ’gap’ between seismic ray theory and ’full’ wavefield extrapolation, Geophys. J. Int., 137, 364–380, 1999. Thomson, C. J., Accuracy and efficiency considerations for wide-angle wavefield extrapolators and scattering operators, Geophys. J. Int., 163, 308–323, 2005. 164 Appendix B Supplementary Information for Chapter 5 B.1 Receiver functions Three-component seismograms are collected for all events with magnitude > 5.5 in the epicentral distance range 30-100◦. The vertical and horizontal (radial and transverse) components of motion are decomposed into upgoing P, SV (radial), and SH (transverse) wave components using the free-surface transfer matrix [Bostock , 1998]. We retain seismograms with P-component signal-to-noise ratios > 10 dB of the first 10 seconds in the 0-2 Hz frequency band. Recorded events come mainly from the Western Pacific, Fiji-Tonga, and Central and South American regions. Backazimuthal coverage is good from 225◦ to 360◦ and from 125◦ to 175◦, and less regularly sampled other- wise. Because the strike of Cascadia subduction is approximately NW-SE, slab structure is sampled mostly from the updip and strike parallel direc- tions. Individual single-event seismograms are processed using the receiver func- tion method, which uses the P component as an estimate of the source to deconvolve the SV and SH components and recover the receiver-side S scattering structure. This procedure is performed using a Wiener spectral deconvolution [Bostock , 1998]. The resulting SV and SH receiver function time series represent forward (Pxs) and backscattered (Ppxs, Psxs) waves from discontinuities in physical properties, identified by x (Fig. B.1). The timing and amplitude of each converted phase constrain overlying velocity structure (depth and VP/VS of the overlying column) and velocity contrast, 165 Appendix B. Supplementary Information for Chapter 5 Surface Station Pps Pss P Ps Figure B.1: Ray diagram for P-to-S converted phases and free-surface mul- tiples from a single dipping interface, for example the oceanic Moho. Rays travelling as P waves are represented by solid lines, whereas rays represent- ing S travelling waves are shown as dashed lines. For a dipping, low-velocity layer bounded by parallel interfaces, this set of wave conversions is gener- ated twice with opposite polarities for waves scattered at the top and bottom layer boundaries, as in Fig. B.2. respectively. To visualize coherent structural signals beneath individual sta- tions we sort receiver functions by backazimuth of incoming wavefield and display the time series as adjacent blue-red colorcoded traces (Fig. B.2). Amplitudes are relative to incident P-wave. These data are sensitive to structures with scale lengths of 1-10 km and are dominated by the signature of an eastward dipping low-velocity zone (LVZ). This signature includes 3 sets of oppositely polarized pulses comprising forward scattered Pxs conver- sions, and backscattered Ppxs and Psxs conversions from the top (x=t) and bottom (x=b) of the LVZ. 166 Appendix B. Supplementary Information for Chapter 5 0 5 10 15 20 25T im e [se c] a −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0 5 10 15 20 25T im e [se c] b −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0 60 120 180 240 300 360 BA Z [de g] c 0.04 0.05 0.06 0.07 0.08 p x [s/ km ] d 0 5 10 15 20 25 Pxs Ppxs Psxs -0.4 -0.2 0 0.2 0.4 e b b b ttt Figure B.2: Radial(a) and transverse (b) component receiver functions at station PGC (southern Vancouver Island) sorted by backazimuth (c) of inci- dent wavefields. (d) shows distribution in horizontal slowness. (e) illustrates arrivals of converted phases (Pxs, Ppxs, Psxs) from the top [t] and bottom [b] of the LVZ within the radial-component (a) receiver functions. 167 Appendix B. Supplementary Information for Chapter 5 B.2 Phase stacking The stacking technique is based on the work by Zhu and Kanamori [2000] for horizontal layering, and later adapted by Rossi et al. [2006] for planar dipping layers. The procedure is based on stacking moveout-corrected traces for each weighted radial and transverse converted phase (Pxs, Ppxs, Psxs) in the receiver functions for a given model, and searching through a range of parameters. Because the moveout depends on the propagation path through the model for a given incident wavefield, all data at a given station can be used simultaneously. By stacking forward and backscattered phases with different moveout curves, the stacks will interfere constructively and produce a maximum at the correct model. This procedure is effective for simple layering and isotropic velocity struc- ture. If anisotropy is present, the periodic backazimuthal response of con- verted phases can in practice be corrected to account for polarity reversals, which would otherwise produce destructive interference in the stacks and thus bias the results (Fig. B.3). Alternatively, one can carefully select phases that are less affected by anisotropy, for example by applying larger relative weight to those phases, in order to avoid this effect. Strike direction along the margin is constrained by polarity reversals on transverse-component receiver functions (Fig. B.2,B.3). We use the phase- weighted stacking technique of Schimmel and Paulssen [1997] for individual phases but employ the median instead of the mean to improve stack co- herency. Our first target is the positive velocity contrast representing the Moho of the subducted oceanic plate [Abers, 2005] because it is thought to be the largest velocity contrast and thus afford the maximum amplitude in the stacks. We use a dipping, single-layer model with mantle and crustal P velocities of 7.8 and 6.5 km/s [Ramachandran et al., 2006] to compute trav- eltimes, although the results are only weakly dependent on this assumption [Zandt and Ammon, 1995]. Traveltimes are calculated from analytical solu- tions for a dipping interface [Frederiksen et al., 2001]. The results at station PGC show that, on top of recovering depth and VP /VS from the strong posi- tive maximum (bottom of the LVZ), we also recover the secondary, shallower 168 Appendix B. Supplementary Information for Chapter 5 0 5 10 15 Ti m e [se c] SV component 0 5 10 15 Ti m e [se c] 30 60 90 120 150 180 210 240 270 300 330 Back−Azimuth [deg] SH component Figure B.3: Averaged radial (SV) and transverse (SH) component receiver functions within 10◦ bins for stations VI07, VI08, VI09, VI53, VI54, and WOS. One polarity reversal in Pxs signals (26 s) on the SV-component is clearly evident at around 270◦. Polarity reversals in SH-component appear at different backazimuths due to strike of the LVZ. SV-component Pxs sig- nals from these stations thus manifest seismic anisotropy. We therefore give no weight to SV-component signals to avoid bias in the stacks. Color scale is the same as in Fig. B.2. 169 Appendix B. Supplementary Information for Chapter 5 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 Vp /V s Ps (SV) Ps (SH) Pps (SV) 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 Vp /V s 15 20 25 30 35 40 45 50 55 Depth Pps (SH) 15 20 25 30 35 40 45 50 55 Depth Pss (SV) 15 20 25 30 35 40 45 50 55 Depth Weighted average Figure B.4: Example of phase stacking for station PGC. In this case each phase is given equal weight and the moveouts are calculated for a 14◦ dipping discontinuity. The positive and negative extrema represent the solutions for the depth to bottom and top of the low-velocity zone, respectively, and the VP /VS of the overlying column. and lower VP /VS , negative extremum from a negative velocity contrast at the top of the LVZ (Fig. B.4). B.3 Tomography cross-sections We show 2 cross-sections of the tomographic model along lines A1-A2 and B1-B2 (Fig. 5.3,B.5,B.6). The positive velocity anomaly representing the thermal and compositional signature of the subducted Juan de Fuca slab extends to a depth of 300 km. The topography shown on top is exaggerated. 170 Appendix B. Supplementary Information for Chapter 5 Figure B.5: Cross-section of the tomographic model along line A1-A2. The slab extends to a depth of 300 km. Its association with deeper anomalous material (positive velocity anomaly at 300-400 km) is unclear. 171 Appendix B. Supplementary Information for Chapter 5 Figure B.6: Cross-section of the tomographic model along line B1-B2. The slab extends to a depth of 300 km. 172 Bibliography Abers, G. A., Seismic low-velocity layer at the top of subducting slabs beneath volcanic arcs: observations, predictions, and systematics, Phys. Earth Planet. Int., 149, 7–29, doi:10.1016/j.pepi.2004.10.002, 2005. Bostock, M. G., Mantle stratigraphy and evolution of the Slave province, J. Geophys. Res., 103, 21,183–21,200, 1998. Frederiksen, A. W., M. G. Bostock, and J. F. Cassidy, S-wave velocity structure of the Canadian upper mantle, Phys. Earth Planet. Inter., 124, 175–191, 2001. Ramachandran, K., R. D. Hyndman, and T. M. Brocher, Regional p wave velocity structure of the Northern Cascadia Subduction Zone, J. Geo- phys. Res., 111, B12,301, doi:10.1029/2005JB004,108, 2006. Rossi, G., G. A. Abers, S. Rondenay, and D. H. Christensen, Unusual mantle Poisson’s ratio, subduction, and crustal structure in central Alaska, J. Geophys. Res., 111, B09,311, doi:10.1029/2005JB003,956, 2006. Schimmel, M., and H. Paulssen, Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int., 130, 497– 505, 1997. Zandt, G., and C. J. Ammon, Continental crust composition constrained by measurements of crustal Poisson’s ratio, Nature, 374, 152–154, 1995. Zhu, L., and H. Kanamori, Moho depth variation in southern California from teleseismic receiver functions, J. Geophys. Res., 105, 2969–2980, 2000. 173 Appendix C Supplementary Information for Chapter 6 C.1 Receiver function analysis We use high signal-to-noise ratio (> 10) broadband seismic recordings from 369, M > 6 earthquakes at teleseismic distances that occured between 2002 and 2006. Three-component seismograms are decomposed into P and S upgoing wave modes and deconvolved using Wiener spectral deconvolution to obtain radial and transverse P receiver functions (RFs). At each station we used an average of 100 highest quality RFs which were sorted by back- azimuth of incident wavefield and filtered at corner frequencies of 0.05 to 0.5 Hz, 0.05 Hz to 1 Hz, and 0.05 Hz to 2 Hz (Fig. C.1). The resulting RFs show coherent signals representing direct P-to-S (Ps) conversions and free- surface P-to-S (Pps) and S-to-S (Pss) reverberations from velocity contrasts within the underlying column that are most easily seen at lower frequencies (Fig. C.1A). Each scattered phase displays oppositely-polarized (blue-red for Ps and Pps, red-blue for Pss) pulses that are characteristic of a prominent, dipping low-velocity layer. C.2 Phase weighted stacking The timing of scattered modes (Ps, Pps, Pss) relative to P can be used to estimate depth, and dip of major interfaces and the average VP /VS of the overlying column, [Zhu and Kanamori , 2000; Rossi et al., 2006], by stacking moveout-corrected RFs from a range of models and identifying which model produces the maximum constructive interference. Strike direction (N45W) 174 Appendix C. Supplementary Information for Chapter 6 0 5 10 15 20 25 30 Ti m e [se c] A −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 PFBTWG B TWB B TSJBLZBTWK B MGC B KELBPGC SSIB SILB GOW B SNB 0 5 10 15 20 25 30 Ti m e [se c] B −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 PFBTWG B TWB B TSJB LZBTWK B MGC BKELBPGC SSIB SILB GOW B SNB 0 5 10 15 20 25 30 Ti m e [se c] C −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 PFBTWG B TWB B TSJBLZBTWK B MGC BKELB PGC SSIB SILB GOW B SNB Figure C.1: Receiver functions for all data used in the analysis sorted by station position along the line and, for each station, by back-azimuth of incident wavefield. Receiver functions are filtered at cutoff frequencies of: A) 0.05 - 0.5 Hz; B) 0.05 - 1 Hz; C) 0.05 - 2 Hz. Only data from B) and C) are used in the analysis; A) is plotted for visual purposes. Amplitudes are relative to incident P wave. Note that distances are not preserved on this plot. 175 Appendix C. Supplementary Information for Chapter 6 is constrained by polarity reversals on transverse RFs. We adopt the phase- weighted stacking technique of Schimmel and Paulssen [1997] but employ the median instead of the mean to improve stack coherency. Our first tar- get is the positive velocity contrast representing the Moho of the subducted oceanic plate [Abers, 2005]. We use a dipping, single-layer model with man- tle and crustal P velocities of 7.8 and 6.5 km/s [Ramachandran et al., 2006] to compute travel-times, although the results are only weakly dependent on this assumption [Zandt and Ammon, 1995]. Our results suggest that the material overlying the oceanic Moho (oceanic crust + continental crust) has a vertically-averaged VP/VS of 1.82 ± 0.07 for the range of stations along the array (Fig. C.2). On the stacks, the strong negative trough corresponding to the top of low-velocity oceanic crust appears at shallower depths with an average VP /VS of overlying material (continental crust) of 1.73 ±0.07 along the array that is consistently lower by an average of 0.09 than VP /VS of material overlying the oceanic Moho. Discrepancy in the average VP /VS of the overlying column suggests that the low-velocity layer bounded by both discontinuities has significantly higher VP /VS than overlying material. The vertical thickness of the low-velocity layer is approximately constant along the array and has mean value of 4.6 ±1.1 km. Based on these findings we further investigate the structure of the low- velocity layer by studying the relative travel-times between positive and negative signals for each scattered phase. C.3 Auto-correlation stacking We adopt a similar stacking methodology to characterize the layer interval Poisson’s ratio by incorporating a few modifications. In the usual stack- ing approach [Zhu and Kanamori , 2000] travel-times are computed relative to an absolute 0 time that is conveniently provided through deconvolution. VP /VS for the low-velocity layer is sensitive to relative travel-time differ- ences from two scattered arrivals (for a given scattering mode) originating at the top and bottom of that layer. We obtain time series of differen- tial S-phase time lags by extracting oppositely-polarized pulses from filtered 176 Appendix C. Supplementary Information for Chapter 6 0 10 20 30 40 50 60 70 80 D ep th [k m] A PFB TWG B TWB B TSJB LZB TWK B MGC B KELB PGCSSIBSILB GOW BSN 1.4 1.6 1.8 2.0 2.2 V p /V s 150 160 170 180 190 200 210 220 230 240 250 Distance from trench [km] B Figure C.2: Depth (A) to top and bottom interface and Vp/Vs (B) of over- lying column obtained by stacking receiver functions filtered between 0.05 - 1 Hz. 177 Appendix C. Supplementary Information for Chapter 6 RF windows for a given scattering mode centered around the solution for the oceanic Moho found previously, and calculating the normalized auto- correlation (Fig. C.3, C.4). All positive peaks are zeroed out to eliminate interference from arrivals that do not conform with the LVZ model, and the auto-correlations are shifted and stacked using travel- time expressions for a dipping, low-velocity layer with fixed dip and P and S velocity estimates for both underlying and overlying media (Fig. C.5). The stacking approach benefits from efficient computation of travel-times through dipping layer models. In the next section we develop analytic ex- pressions for plane-wave travel-times that avoid solution of an eigenvalue problem [Frederiksen and Bostock , 2000]. C.4 Plane-wave travel-times Let us consider the propagation of a plane P-wave incident from below upon a stack of planar, dipping (and not necessarily parallel) layers. The timing of the Ps phase from the Lth and L + 1st interfaces can be written as (see Frederiksen and Bostock [2000]) tLPs = L∑ i=1 ∆zi[ξ̂ S i − ξ̂Pi ] (C.1) tL+1Ps = L+1∑ i=1 ∆zi[ξ̂ S i − ξ̂Pi ], (C.2) respectively, where ξ̂i is the absolute value of the vertical component of slowness of the appropriate (P or S) upgoing wave in the ith layer, and ∆zi is layer thickness. The differential travel-time ∆t L Ps between two Ps conversions from the top and bottom of the Lth layer is thus ∆tLPs = ∆zL[ξ̂ S L − ξ̂PL ]. (C.3) 178 Appendix C. Supplementary Information for Chapter 6 0 5 10 15 20 25T im e [se c] A −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0 1 2 3 Ti m e [se c] B −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0 1 2 3 Ti m e [se c] C −1.00 −0.75 −0.50 −0.25 0.00 0 60 120 180 240 300 360 BA Z [de g] D 0.04 0.05 0.06 0.07 0.08 p x [s/ km ] E Figure C.3: Window extraction and auto-correlation of Ps phase filtered with cutoff frequency of 2 Hz at station PGC. In A) a window (bounded by black lines) is chosen on the basis of predicted Ps arrival times for the subducted oceanic Moho (Fig. C.2). The extracted window (B) is tapered and autocorrelated (C). Each autocorrelation is normalized by positive peak amplitude at 0 time (not shown here) and all positive peaks are subsequently zeroed out. D) and E) show the back-azimuthal and horizontal slowness distribution of incident wavefields for each corresponding trace plotted in A). 179 Appendix C. Supplementary Information for Chapter 6 0 5 10 15 20 25T im e [se c] A −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0 1 2 3 4 5 Ti m e [se c] B −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0 1 2 3 4 5 Ti m e [se c] C −1.00 −0.75 −0.50 −0.25 0.00 0 60 120 180 240 300 360 BA Z [de g] D 0.04 0.05 0.06 0.07 0.08 p x [s/ km ] E Figure C.4: Same as Fig. C.3 for Pps phase and cutoff frequency of 1 Hz at station PGC. 180 Appendix C. Supplementary Information for Chapter 6 0.25 0.30 0.35 0.40 Po is so n’ s ra tio 3 4 5 6 7 8 Thickness [km] Ps 3 4 5 6 7 8 Thickness [km] Pps 3 4 5 6 7 8 Thickness [km] Average Figure C.5: Autocorrelation stacking results for station PGC. We use only Ps and Pps phases from FIg. S3 and Fig. S4 because Pss auto-correlations are not coherent at 1 Hz. Each phase is given equal weight and the minimum of the averaged stack gives the estimated Poisson’s ratio and thickness of low-velocity layer. Similar expressions hold for Pps and Pss reverberations, namely ∆tLPps = ∆zL[ξ̂ S L + ξ̌ P L ] (C.4) ∆tLPss = ∆zL[ξ̂ S L + ξ̌ S L], (C.5) where ξ̌ signifies absolute value of vertical component of slowness for a down- going wave. In a 1-D medium, these expressions simplify to ∆tLPs = ∆zL [√ R2L/V 2 P,L − p2 − √ 1/V 2P,L − p2 ] (C.6) ∆tLPps = ∆zL [√ R2L/V 2 P,L − p2 + √ 1/V 2P,L − p2 ] (C.7) ∆tLPss = ∆zL [ 2 √ R2L/V 2 P,L − p2 ] , (C.8) where p is horizontal slowness of the incident wavefield, and R = VP /VS of the layer. Note that only two equations are independent in this case since ∆tLPss = ∆t L Ps +∆t L Pps. For a planar dipping layer with parallel boundaries and unit normal 181 Appendix C. Supplementary Information for Chapter 6 n = (n1, n2, n3), these expressions generalize to ∆tLPs = n3∆zL [√ R2 V 2P,L − 1 V 2P,L+1 + (n · pI)2 − √ 1 V 2P,L − 1 V 2P,L+1 + (n · pI)2 ] (C.9) ∆tLPps = n3∆zL [√ R2 V 2P,L − 1 V 2P,L−1 + (n · pR)2 + √ 1 V 2P,L − 1 V 2P,L−1 + (n · pR)2 ] (C.10) ∆tLPss = n3∆zL [ 2 √ R2 V 2P,L − 1 V 2S,L−1 + (n · pS)2 ] , (C.11) where the slowness vectors pI , pR, pS correspond to the incident P-wave below the layer, the free-surface refelcted P-wave above the layer and the free-surface reflected S-wave above the layer, respectively. pR and pS are computed by applying Snell’s law across the interfaces. These travel-times depend on the assumed background velocity model, but, as before, only weakly. C.5 Limitating cases To characterize R (or Poisson’s ratio) within the layer the key parameter is the ratio between ∆tPs and ∆tPps (or ∆tPss). Consider a single, horizontal layer located within an otherwise homogeneous half-space, and a normally- incident plane P-wavefield (ignoring for a moment that no conversions are generated in such a circumstance). In this case the ratio of equations (6) and (7) is independent of VP and ∆z and yields X ≡ ∆tPs ∆tPps = (R− 1) (R+ 1) . (C.12) 182 Appendix C. Supplementary Information for Chapter 6 Solving for R we have R = 1 +X 1−X . (C.13) As X tends to zero (i.e. ∆tPs tends to 0 or ∆tPps ≫ ∆tPs), R tends to 1, and VP → VS . As X tends to one (i.e. ∆tPs → ∆tPps), R tends to ∞, and VS → 0. Poisson’s ratio (σ) is calculated from R using σ = 1 2 [ 1− 1 R2 − 1 ] (C.14) C.6 Frequency bias To avoid bias in the determination of R, one must be careful not to low- pass filter RFs at corner periods larger than twice the true time separation between pulses (Fig. C.6). If ∆tPs is biased to higher values (i.e. cutoff frequency is too low), both X and R will be biased upward. The opposite relation is true for the reverberations. For example, if ∆tPps is biased up- ward, both X and R will be biased downward. It is thus crucial that the time series for Ps be filtered at frequencies high enough to resolve the small- est differential travel-times between the oppositely-polarized signals. Using synthetic tests we verified that a cutoff frequency of 2 Hz is appropriate and sufficient for our purposes. Because Pps and Pss reverberations are more band-limited, we use a cutoff frequency of 1 Hz. Using slightly lower cutoff frequencies for Pps and Pss will not affect the recovery of R because X is always smaller than 1, and any bias will yield more conservative estimates of VP /VS . C.7 Bootstrapping Standard errors are calculated using the bootstrap method [Efron and Tib- shirani , 1991]. For each set of receiver functions (or auto-correlations) we randomly resample the original data set with replacement and perform phase-weighted stacking. A best fit value is determined for the resampled 183 Appendix C. Supplementary Information for Chapter 6 Figure C.6: Illustration of frequency bias due to filtering. T is corner period of filter, ∆t is minimum value of recovered time difference. 184 Appendix C. Supplementary Information for Chapter 6 data set and this procedure is repeated 100 times to calculate the standard deviation around the solution. 185 Bibliography Abers, G. A., Seismic low-velocity layer at the top of subducting slabs beneath volcanic arcs: observations, predictions, and systematics, Phys. Earth Planet. Int., 149, 7–29, doi:10.1016/j.pepi.2004.10.002, 2005. Efron, B., and R. Tibshirani, Statistical data analysis in the computer age, Science, 253, 390–395, 1991. Frederiksen, A. W., and M. G. Bostock, Modelling teleseismic waves in dipping anisotropic structures, Geophys. J. Int., 141, 401–412, 2000. Ramachandran, K., R. D. Hyndman, and T. M. Brocher, Regional p wave velocity structure of the Northern Cascadia Subduction Zone, J. Geo- phys. Res., 111, B12,301, doi:10.1029/2005JB004,108, 2006. Rossi, G., G. A. Abers, S. Rondenay, and D. H. Christensen, Unusual mantle Poisson’s ratio, subduction, and crustal structure in central Alaska, J. Geophys. Res., 111, B09,311, doi:10.1029/2005JB003,956, 2006. Schimmel, M., and H. Paulssen, Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int., 130, 497– 505, 1997. Zandt, G., and C. J. Ammon, Continental crust composition constrained by measurements of crustal Poisson’s ratio, Nature, 374, 152–154, 1995. Zhu, L., and H. Kanamori, Moho depth variation in southern California from teleseismic receiver functions, J. Geophys. Res., 105, 2969–2980, 2000. 186
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Seismic and mechanical attributes of lithospheric deformation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Seismic and mechanical attributes of lithospheric deformation and subduction in western Canada Audet, Pascal 2008
pdf
Page Metadata
Item Metadata
Title | Seismic and mechanical attributes of lithospheric deformation and subduction in western Canada |
Creator |
Audet, Pascal |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | Convergent continental margins are regions of intense deformation caused by the interaction of oceanic plates with continents. The spatial extent of deformation is broadly commensurate with the specific time scale of the causative phenomenon. For example, subduction-related short-term deformation is limited to <200 km from the margin, whereas long-term plate convergence cause deformation over ∼1000 km landward. Deformation is thus manifested in multiple ways, with attributes depending on the scale of measurement. In this thesis we investigate the use of two geophysical approaches in the study of deformation: 1) The analysis of potential-field anomalies to derive estimates of the elastic thickness (Te) of the lithosphere, and 2) The structural study of past and present subduction systems using seismic observations and modelling. Both approaches involve the development of appropriate methodologies for data analysis and modelling, and their application to the western Canadian landmass. Our findings are summarized as follows: 1) We develop a wavelet-based technique to map variations in Te and its anisotropy; 2) We show how a step-wise transition in Te and its anisotropy from the Cordillera to the Craton is a major factor influencing lithospheric deformation; 3) We implement a waveform modelling tool that includes the effects of structural heterogeneity and anisotropy for teleseismic applications, and use it to model the signature of a fossil subduction zone in a Paleoproterozoic terrane; 4) We use teleseismic recordings to map slab edge morphology in northern Cascadia and show how slab window tectonism and slab stretching led to the creation of the oceanic Explorer plate; 5) We use seismic signals from the subducting oceanic crust to calculate elevated Poisson’s ratio and infer high pore-fluid pressures and a low-permeability plate boundary within the forearc region of northern Cascadia. |
Extent | 6481632 bytes |
Subject |
Plate tectonics Deformation |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-10-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0052387 |
URI | http://hdl.handle.net/2429/2435 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_fall_audet_pascal.pdf [ 6.18MB ]
- Metadata
- JSON: 24-1.0052387.json
- JSON-LD: 24-1.0052387-ld.json
- RDF/XML (Pretty): 24-1.0052387-rdf.xml
- RDF/JSON: 24-1.0052387-rdf.json
- Turtle: 24-1.0052387-turtle.txt
- N-Triples: 24-1.0052387-rdf-ntriples.txt
- Original Record: 24-1.0052387-source.json
- Full Text
- 24-1.0052387-fulltext.txt
- Citation
- 24-1.0052387.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0052387/manifest