UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Improved teleseismic Green's functions and western Canada mantle structure and evolution Mercier, Jean-Philippe 2008

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

Item Metadata

Download

Media
24-ubc_2009_spring_mercier_jean_philippe.pdf [ 21.94MB ]
Metadata
JSON: 24-1.0052560.json
JSON-LD: 24-1.0052560-ld.json
RDF/XML (Pretty): 24-1.0052560-rdf.xml
RDF/JSON: 24-1.0052560-rdf.json
Turtle: 24-1.0052560-turtle.txt
N-Triples: 24-1.0052560-rdf-ntriples.txt
Original Record: 24-1.0052560-source.json
Full Text
24-1.0052560-fulltext.txt
Citation
24-1.0052560.ris

Full Text

Improved teleseismic Green’s functions and western Canada mantle structure and evolution by Jean-Philippe Mercier ´ B.Eng., Ecole Polytechnique de Montr´eal, 2001 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Earth and Ocean Sciences)  The University Of British Columbia (Vancouver) November, 2008 c Jean-Philippe Mercier 2008  Abstract The present thesis is divided into three distinct parts and focuses both on the improvement of deconvolution technique in a teleseismic context for crustal and upper-mantle studies and on the understanding of western Canada structure and evolution through seismic imaging. The first part presents estimates of the P-component of the teleseismic-P Green’s functions for three stations of the Canadian National Seismic Network (CNSN) obtained using a new deconvolution technique. Our results show evidence of the principal, first-order scattered Moho phases and, in particular, the P pM p. The second part presents teleseismic receiver functions from 20 broadband three-component seismometers deployed along the MacKenzie-Liard Highway in Canada’s Northwest Territories as part of the joint LithoprobeIRIS CAnadian NOrthwest Experiment (CANOE). These stations traverse a Paleoproterozoic suture and subduction zone that has been previously documented in detail to mantle depths using seismic reflection profiling. Our results reveal the response of the ∼1.8 Ga subduction zone on both the radial and transverse component. The identification of this structure and its comparison with fine-scale mantle layering below the adjacent Slave province and from a range of Precambrian terranes provides an unambiguous connection between fossil subduction and fine-scale, anisotropic mantle layering found beneath cratons. Previous documentation of similar layering below the adjacent Slave province provides strong support for the thesis that early cratonic blocks were stabilized through processes of shallow subduction. The last part presents P - and S wave velocity models for western Canada. In this part, we focus our attention on two distinct features: 1) the transition from Phanerozoic to cratonic mantle in northwestern Canada and 2) the complex tectonic environment at the northern terminus of the Cascaii  Abstract dia subduction zone where the plate boundary changes from convergent to transform. We find that the main transition from Phanerozoic to cratonic mantle in northwestern Canada occurs at the Cordilleran deformation front. In northern Cascadia, we have imaged and characterized the signature of the subducting Juan de Fuca plate and observed evidence of subduction beyond the northern edge of the slab. Our result show that the Anahim hotspot track is underlain by a -2% low-velocity zone.  iii  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  List of Tables  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Preface  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii  Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3  2 Improved Green’s functions for passive source structural studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5  2.1  Introduction  . . . . . . . . . . . . . . . . . . . . . . . . . . .  5  2.2  Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . .  8  2.2.1  Blind Deconvolution . . . . . . . . . . . . . . . . . . .  9  2.2.2  Multichannel Deconvolution  2.3  . . . . . . . . . . . . . .  12  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  16  2.3.1  Station ULM (Lac du Bonnet, Manitoba) . . . . . . .  18  2.3.2  Station WHY (Whitehorse, Yukon)  22  Results  . . . . . . . . . .  iv  Table of Contents 2.3.3  . . . . . . . . . .  22  Discussion and Concluding Remarks . . . . . . . . . . . . . .  24  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  27  2.4  Station EDM (Edmonton, Alberta)  3 The teleseismic signature of fossil subduction: Northwestern Canada  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.1  Introduction  3.2  Tectonic  3.3  Data  30  . . . . . . . . . . . . . . . . . . . . . . . . . . .  30  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  33  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  33  3.4  Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . .  34  3.5  Results  36  3.6  Modelling and Numerical Simulation  3.7  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  3.6.1  Inversion . . . . . . . . . . . . . . . . . . . . . . . . .  42  3.6.2  One-way modeling . . . . . . . . . . . . . . . . . . . .  45  Interpretation  . . . . . . . . . . . . . . . . . . . . . . . . . .  49  3.7.1  SNORCLE active source studies . . . . . . . . . . . .  49  3.7.2  Comparison of CANOE and Yellowknife array results  54  3.8  Discussion and Conclusions . . . . . . . . . . . . . . . . . . .  55  3.9  Acknowledgments  . . . . . . . . . . . . . . . . . . . . . . . .  60  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  62  4 Body-wave tomography of western Canada  . . . . . . . . .  66  . . . . . . . . . . . . . . . . . . . . . . . . . . .  66  4.1  Introduction  4.2  Data and Method  4.3  . . . . . . . . . . . . . . . . . . . . . . . .  67  . . . . . . . . . . . . . . . . . . . . . . . . . . .  67  4.2.1  Data  4.2.2  Model parametrization  . . . . . . . . . . . . . . . . .  70  4.2.3  Traveltime inversion . . . . . . . . . . . . . . . . . . .  70  4.2.4  Choice of inversion parameters . . . . . . . . . . . . .  72  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  73  Results 4.3.1  Resolution Test  . . . . . . . . . . . . . . . . . . . . .  4.3.2  Raw average station residuals  . . . . . . . . . . . . .  73 75  v  Table of Contents 4.3.3  P and S Upper Mantle Velocity structure beneath western Canada  4.4  Discussion  . . . . . . . . . . . . . . . . . . . . .  77  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  81  4.4.1  Cordillera/craton transition  . . . . . . . . . . . . . .  81  4.4.2  Subsurface geometry of Northern Cascadia . . . . . .  86  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .  95  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  96  4.5  5 Conclusion 5.1  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103  Key Contributions . . . . . . . . . . . . . . . . . . . . . . . . 103 5.1.1  Improved Green’s functions for passive source structural studies  5.1.2  . . . . . . . . . . . . . . . . . . . . . . . 104  The teleseismic signature of fossil subduction: Northwestern Canada  5.1.3 5.2  . . . . . . . . . . . . . . . . . . . . . 104  Body-wave tomography of western Canada . . . . . . 105  Limitations and Future work . . . . . . . . . . . . . . . . . . 106 5.2.1  Improved Green’s functions for passive source structural studies  5.2.2  . . . . . . . . . . . . . . . . . . . . . . . 106  The teleseismic signature of fossil subduction: Northwestern Canada  5.2.3  . . . . . . . . . . . . . . . . . . . . . 106  Body-wave tomography of western Canada . . . . . . 107  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108  Appendices A Tomography Resolution Test  . . . . . . . . . . . . . . . . . . 109  vi  List of Tables 3.1  Model parameters . . . . . . . . . . . . . . . . . . . . . . . .  46  vii  List of Figures 2.1  First order scattered phased in a layer over half-space and resulting synthetic P and SV -component of the teleseismic-P of the Green’s function . . . . . . . . . . . . . . . . . . . . . .  2.2  7  Comparison of true and estimated cross-correlogram using a white and a statistically white spectrum . . . . . . . . . . . .  11  2.3  Whitening parameter effect . . . . . . . . . . . . . . . . . . .  13  2.4  First-order crustal scattered phases as a function of slowness for a layer over half-space model . . . . . . . . . . . . . . . .  2.5  Estimate of the SV-component of the Green’s function for station ULM . . . . . . . . . . . . . . . . . . . . . . . . . . .  2.6  23  Estimate of the P and SV-component of the Green’s function for station EDM . . . . . . . . . . . . . . . . . . . . . . . . .  3.1  21  Estimate of the P and SV-component of the Green’s function for station WHY . . . . . . . . . . . . . . . . . . . . . . . . .  2.8  20  Estimate of the P-component of the Green’s function for station ULM . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  2.7  18  25  Map of western Canadian geological provinces illustrating the distribution of the broad-band three-component CANOE stations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  32  3.2  Event distribution map  . . . . . . . . . . . . . . . . . . . . .  35  3.3  Illustration of the displaying approach . . . . . . . . . . . . .  37  3.4  Receiver function results . . . . . . . . . . . . . . . . . . . . .  39  3.5  Transverse component azimuthal response . . . . . . . . . . .  41  3.6  Cartoon of the proposed subsurface model . . . . . . . . . . .  44  viii  List of Figures 3.7  Comparison between real and synthetic back-azimuthal response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.8  46  Subsurface model of the Wopmay Orogen showing the geometry of the dipping layer . . . . . . . . . . . . . . . . . . . . .  47  Comparison of synthetic and observed receiver functions . . .  48  3.10 Comparison of synthetic and real response . . . . . . . . . . .  50  3.9  3.11 Near-vertical reflection superimposed on top of the receiver function results . . . . . . . . . . . . . . . . . . . . . . . . . .  52  3.12 Comparison of H and AM response . . . . . . . . . . . . . .  55  4.1  Map of western Canada illustrating the distribution of braodband and short-period seismic stations used in this study . .  68  4.2  Event distribution map  . . . . . . . . . . . . . . . . . . . . .  69  4.3  Model grid  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  71  4.4  Ray density map . . . . . . . . . . . . . . . . . . . . . . . . .  74  4.5  Resolution test input model . . . . . . . . . . . . . . . . . . .  75  4.6  Resolution test results . . . . . . . . . . . . . . . . . . . . . .  76  4.7  Relative station residuals . . . . . . . . . . . . . . . . . . . .  78  4.8  Horizontal slices within the P velocity model for the whole survey area . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4.9  79  Horizontal slices within the S velocity model for the whole survey area . . . . . . . . . . . . . . . . . . . . . . . . . . . .  80  4.10 Cross section within the P velocity showing the location of the Cordillera/craton transition . . . . . . . . . . . . . . . . .  83  4.11 Plot of the velocity variation at 100 km depth in northwestern Canada centered on the Cordillera deformation front . . . . .  85  4.12 Depth slice at 100 km depth within the P velocity model in southwestern British Columbia with tectonics of northern Cascadia detailed . . . . . . . . . . . . . . . . . . . . . . . . .  87  4.13 Vertical cross sections within the P velocity model in the northern Cascadia region. . . . . . . . . . . . . . . . . . . . . 4.14 Vertical cross sections within the P velocity model.  . . . . .  89 93  ix  List of Figures A.1 Input P resolution test model for the depth slices presented in Figure 4.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 A.2 Output of the resolution test for the P data-set for the depth slices presented in Figure 4.8. . . . . . . . . . . . . . . . . . . 111 A.3 Input S resolution test model for the depth slices presented in Figure 4.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 A.4 Output of the resolution test for the S data-set for the depth slices presented in Figure 4.9. . . . . . . . . . . . . . . . . . . 113 A.5 Input P resolution test model for the cross sections presented in Figure 4.10. . . . . . . . . . . . . . . . . . . . . . . . . . . 114 A.6 Output of the resolution test for the S data-set for the depth slices presented in Figure 4.10. . . . . . . . . . . . . . . . . . 114 A.7 Input P resolution test model for the cross sections presented in Figure 4.13. . . . . . . . . . . . . . . . . . . . . . . . . . . 115 A.8 Output of the resolution test for the P data-set for the depth slices presented in Figure 4.13. . . . . . . . . . . . . . . . . . 116 A.9 Input P resolution test model for the cross sections presented in Figure 4.14. . . . . . . . . . . . . . . . . . . . . . . . . . . 117 A.10 Output of the resolution test for the P data-set for the depth slices presented in Figure 4.14. . . . . . . . . . . . . . . . . . 118  x  Preface This thesis is a collection of three published or submitted papers written in collaboration with numbers of coauthors from several institutions both in Canada and the United States. The content of the articles was not altered to comply with journal copyrights. All chapters were formatted to follow the regulations of the Faculty of Graduate studies. The code used to produce the work presented here were written using mostly Matlab and Fortran languages. Most of the figures were produced using Matlab, Generic Mapping Tool (GMT), Inkscape, and XFig.  xi  Acknowledgements First and foremost, I would like to acknowledge the support of Michael Bostock, my supervisor, who during the four and half years of my PhD, has, by his constant support and his invaluable input, immensely contributed this thesis. I would also like to thank the members of my supervisory committee, Ron Clowes, Elizabeth Hearn, and Felix Herrmann for their advice and technical support; my present and former colleagues in the teleseismic group, Pascal Audet, Adam Baig, Julien Chaput, Todd Nicholson and Andrew Schaeffer, whithout whom my life at UBC would not have been so rich; Special thanks to Michael Bostock, Andrew Schaeffer and Pascal Audet, for their precious advice; Andy Langlois and Issamm Al-Koubbi for their remarkable support in the field; James Ni, Eric Sandvol, and Thomas Hearn for the opportunities to see and work on the Tibetan plateau, one of the most remote and beautiful parts of this world; John Cassidy for the access to Nechako data, and George Zandt and Kenneth Dueker for the access to Batholiths data; and Noel Barstow and Jingle Jinwar for the great moments we shared together and for their friendship. J’aimerais aussi remercier mes parents, Claudette et Michel, ainsi que mon fr`ere Jean-Pascal, sans eux, cette th`ese n’aurait sans doute pas vu le jour. Leur amour, conseils et encouragements ont ´et´e d’une valeur inestimable tout au long de ma vie. Je ne peux passer sous silence la contribution de Karine, mon ´epouse, qui a partag´e au cours des quatre derni`eres ann´ees les aleas de la vie de doctorant. Je lui serai toujours reconnaissant de amour et du support dont elle a fait preuve.  xii  Dedication ` Karine, Claudette, Michel, Jean-Pascal et Lysanne. A  xiii  Statement of Co-Authorship Chapter 2 was written in collaboration with Michael Bostock and Adam Baig. Michael identified the research objectives. He proposed a new deconvolution technique. Adam and I worked on its implementation. I applied our implementation to real data from the Canadian National Seismic Network. I wrote the manuscript with input from Michael Bostock. Chapter 3 was written in collaboration with Michael Bostock, Pascal Audet from UBC, James Gaherty from Lamont Doherty Earth Observatory, Edward Garnero from Arizona State University, and Justin Revenaugh from University of Minnesota. Michael Bostock identified the research objective in collaboration with co-principal investigators on the CANOE experiment James Gaherty, Edward Garnero, and Justin Revenaugh. Pascal Audet, produced a synthetic model that was used to compare assess the validity of the model we proposed. I performed the research and I wrote the manuscript with input from Michael. Chapter 4 was written with Michael Bostock and John Cassidy from the Pacific Geoscience Center provided useful comments. Michael, in collaboration with co-principal investigators on the CANOE experiment James Gaherty, Edward Garnero, and Justin Revenaugh, and I identified the research objectives. Kenneth Dueker and George Zandt granted us access to data from the BATHOLITHS experiment and John Cassidy provided data from the POLARIS Nechako experiment. I performed the research and wrote the manuscript with input from Michael.  xiv  Chapter 1  Introduction The present thesis is a collection of three papers that were written in collaboration with researchers from several different institutions and published in or submitted to scientific journals. This thesis presents my major contributions to seismology and tackles some important issues that are both local and global in scope. It focuses on the development and application of innovative signal processing techniques in a teleseismic context and on the study of the lower crust and upper-mantle structure of the western Canada subsurface using receiver functions and travel-time tomography. The first paper featured in this document proposes an original treatment of the deconvolution problem for passive source structural studies. The deconvolution problem is regarded by many seismologists as a key issue and an important step toward the accurate study of the structure and composition of the subsurface. In the teleseismic context, deconvolution aims to unveil the earth response to a quasi-planar, impulsive wave field illuminating seismic discontinuities from below. This signal is often referred to as the teleseismic Green’s function. Nowadays, in teleseismic studies, the deconvolution operation is performed through the so-called “receiver function” analysis. This technique appeared in the mid-1960’s (Phinney, 1964) and reached maturity by the end of 1970’s (Vinnik , 1977; Langston, 1979). The receiver function analysis is now widely used to study the structure and the properties of the earth subsurface in a large variety of settings and contexts (e.g Owens et al. (1984); Ammon and Zandt (1993); Bostock (1997); Zhu and Kanamori (2000); Zandt et al. (2004); Schulte-Pelkum et al. (2005)). They suffer, however, from one major drawback, namely: they preclude the recovery of the P -component of the teleseismic Green’s function and consequently of important subsurface elastic properties. To mitigate this problem, 1  Chapter 1. Introduction several techniques, mainly relying on stacking of seismograms recorded at two or more stations, have been proposed (e.g. Bostock and Sacchi (1997); Revenaugh (1995); Li and Nabelek (1999); Langston and Hammer (2001). However, they all suffer from incorrectly mapping, in the source-time function, the P -component of teleseismic-P Green’s function that display a little move-out from station to station. The new technique proposed in this document successfully extracts information from the P -component of the Green’s function without suffering from the problems encountered by the aforementioned alternative deconvolution techniques. The second paper included in this thesis presents an important contribution to the understanding of the formation of the first fragments of continents, i.e. cratons, during the Archean. This study presents receiver functions computed at broadband seismic stations deployed along a segment of the leg A of the CAnadian NOrthwest Experiment, CANOE, array, in the Wopmay orogen region. It unveils and characterises the teleseismic signature of a Paleo-proterozoic subduction zone that had previously been documented by Cook et al. (1998) along the Lithoprobe SNORCLE transect. This study reveals striking similarities between the teleseismic signature of the Proterozoic subduction beneath the Wopmay orogen and layers infered to be Archean in age documented by Bostock (1997) found beneath the Slave Craton. This observation leads to speculation on the importance of the role played by subduction processes in the stabilization of the cratons. The last paper included in this document, presents a P and S body-wave tomography study of western Canada utilizing data from several permanent and portable seismic networks deployed across the region. This study has allowed the localization of the transition from Phanerozoic to Proterozoic basement in Northwestern Canada. Our velocity models also reveal the structure of the Northern Cascadia region showing evidence for the presence of subducted material north of Vancouver Island. We also image the subsurface beneath the Anahim hotspot which reveals the presence of a large low-velocity anomaly that likely relates to the volcanism observed at the surface.  2  Bibliography Ammon, C., and G. Zandt, Receiver structure beneath the southern Mojave Block, California, B. Seismol. Soc. Am., 83 (3), 737–755, 1993. Bostock, M., and M. Sacchi, Deconvolution of teleseismic recordings for mantle structure, Geophysical Journal International, 129 (1), 143–152, 1997. Bostock, M. G., Anisotropic upper mantle stratigraphy and architecture of the slave craton, Nature, 390, 392–395, 1997. Cook, F. A., A. J. Van der Velden, K. W. Hall, and B. J. Roberts, Tectonic delamination and subcrustal imbrication of the precambrian lithosphere in northwestern canada, mapped by lithoprobe, Geology, 26, 839–842, 1998. Langston, C. A., Structure under mount Rainier, Washington, inferred from teleseismic body waves, Journal of Geophysical Research, 84, 4749– 4762, 1979. Langston, C. A., and J. K. Hammer, The vertical component P-wave receiver function, Bulletin of the Seismological Society of America, 91, 1805– 1819, 2001. Li, X.-Q., and J. L. Nabelek, Deconvolution of teleseismic body waves for enhancing structure beneath a seismometer array, Bulletin of the Seismological Society of America, 89, 190–201, 1999. Owens, T., G. Zandt, and S. Taylor, Seismic evidence for an ancient rift beneath the cumberland plateau, Tennessee: A detailed analysis of broadband teleseismic P waveforms, Journal of Geophysical Research, 89 (B9), 1984. Phinney, R. A., Structure of the Earth’s crust from spectral behavior of long-period body waves, Journal of Geophysical Research, 69, 2997–3017, 1964.  3  Bibliography Revenaugh, J. A., Scattered-wave image of subduction beneath the Transverse Ranges, California, Science, 268, 1888–1892, 1995. Schulte-Pelkum, V., G. Monsalve, A. Sheehan, M. Pandey, S. Sapkota, R. Bilham, and F. Wu, Imaging the Indian subcontinent beneath the Himalaya, Nature, 435 (7046), 1222, 2005. Vinnik, L. P., Detection of waves converted from P to SV in the mantle, Phys. Earth Planet. In, 67, 39–45, 1977. Zandt, G., H. Gilbert, T. Owens, M. Ducea, J. Saleeby, and C. Jones, Active foundering of a continental arc root beneath the southern Sierra Nevada in California, Nature, 431 (7004), 41–46, 2004. Zhu, L., and H. Kanamori, Moho depth variation in southern California from teleseismic receiver functions, Journal of Geophysical Research, 105, 2969–2980, 2000.  4  Chapter 2  Improved Green’s functions for passive source structural studies  2.1  Introduction  For the past few decades, teleseismic receiver functions have provided valuable information on crust and lithospheric structure. In the mid-1960’s, Phinney (1964) proposed that the spectral ratio of horizontal/vertical component seismograms (that is, the frequency domain “receiver function”) could be used to gauge the effect of crustal layering on frequency-dependent polarization. More than a decade later, Burdick and Langston (1977) investigated the time-domain response of impulsive, plane-wave sources at near-vertical incidence, and demonstrated that a full suite of forward- and back-scattered phases could be exploited to better constrain crustal structure. Shortly thereafter, Langston (1979) introduced the time-domain “receiver function” as a means of recovering first-order, P-to-S scattering (note also parallel work by Vinnik (1977) on transition zone discontinuities) and thereby generalized earlier work to complex source-time functions. Inspired by this work, Owens and Zandt (1985) (and countless others after them) A version of this paper has been published in the journal Geophysics in 2006. Mercier, J.-P., M. G. Bostock, and A. M. Baig (2006), Improved Green’s functions for passive source structural studies, Geophysics, 71, doi:10.1190/1.2213951. c 2006 Society of Exploration Geophysics  5  Chapter 2. Improved Green’s functions for passive source structural studies exploited the receiver function methodology to examine the continental crust-mantle boundary by coupling receiver function observations with forward modeled synthetic seismograms for geologically realistic models. They demonstrated that receiver functions could supply information complementary to that produced using active-source seismic methods. Most early teleseismic crustal studies employing receiver functions focused on characterization of the forward-scattered connected PD s phase (where a given discontinuity is designated by “D”) and considered freesurface reflected multiples (or “back scattering”) as a source of noise (see Figure 2.1 for details). The study of the PD s phase in isolation suffers from at least two important limitations. First, there is a significant trade-off between depth to a subsurface discontinuity and the overlying velocity profile. Second, the use of this single phase with a narrow, near-vertical sampling in incidence angle provides information only on the relative S-velocity contrast across an interface, with no constraint placed on either compressional modulus or density. The timing of free-surface multiples can be used to address the first limitation to a certain degree. For example, Zandt and Ammon (1995) investigated bulk crustal composition using estimates of average Poisson’s ratio and crustal thickness. Those quantities were derived using measured times of both forward- and free-surface-reflected, back-scattered waves (that is P pD s and P sD s). This approach has been subsequently generalized and rendered more robust through waveform stacking by Zhu and Kanamori (2000). The amplitudes of the multiples P pD s and P sD s also afford some sensitivity to density contrast but remain insensitive to compressional modulus. Most recently, efforts have been directed toward exploiting receiver functions in the delineation of 2D and 3D Earth structure using migration/inversion methods like those routinely employed in hydrocarbon exploration (Revenaugh, 1995; Dueker and Sheehan, 1997; Bostock et al., 2001; Poppeliers and Pavlis, 2003). The teleseismic application holds the advantage that the full vector wavefield is generally recorded thereby allowing, at least in principle, recovery of a spatially varying suite of elastic parameters. Incorporation of both forward- and back-scattered energy from classic receiver 6  Chapter 2. Improved Green’s functions for passive source structural studies  Figure 2.1: (a) Ray paths for first order scattered phases in a layer over halfspace model. (b) and (c) show synthetic P-and SV-components for incident impulsive P-wave generated for a model that is a Poisson’s solid featuring a 36 km thick crust with P velocity αc = 6.5 km s−1 and density ρc = 2.7 g cm−3 overlying a mantle half space with P velocity αm = 8.1 km s−1 and density ρm = 3.3 g cm−3 . The wavefield is characterized by a constant horizontal slowness of 0.06 s km−1 .  7  Chapter 2. Improved Green’s functions for passive source structural studies functions within an imaging/inversion formulation permits the estimation of two isotropic elastic parameters, namely S-velocity and density, or combinations thereof. Information on subsurface P-impedance structure is also present within the recorded wavefield but its extraction therefrom will demand that we move beyond the classic “receiver function” to a more fundamental quantity, that is, the Earth’s “Green’s function”. This task, in turn, requires that we are better able to remove the effect of the earthquake source-time function. In this paper we demonstrate application of a method described by Baig et al. (2005), which relies on certain spectral properties of the teleseismic wavefield, to recover the P-component of the Green’s function. We begin by providing a motivation and brief explanation of the method including its multichannel implementation. Simple numerical simulations of Green’s functions that identify the key elements of teleseismic scattering from lithospheric discontinuities are then presented. These simulations are used to facilitate the interpretation of deconvolved data computed for 3 stations of the Canadian National Seismograph network that clearly demonstrate the retrieval of improved Green’s function estimates.  2.2  Methodology  This section outlines a new deconvolution technique (Baig et al., 2005) that allows recovery of the teleseismic-P Green’s function and, in particular, its P-component. The method is based on the recognition that the latter quantity is minimum-phase (Sherwood and Trorey, 1965; Bostock , 2004) and, consequently, that it can be reconstructed uniquely and completely from the knowledge of its amplitude spectrum alone (note here that no assumptions are made concerning the nature of the source spectrum only that of the Green’s function). Physically, the minimum phase property stems from the fact that the incident or “ballistic” contribution to the teleseismic-P Green’s function can be characterized by a single, quasi-planar wavefront that contains more energy at all frequencies than the scattered waves that follow in its wake. The validity of this assumption is reaffirmed by the 8  Chapter 2. Improved Green’s functions for passive source structural studies success of receiver function analysis which implicitly relies on the minimumphase property. Since the phase spectrum of a minimum-phase time series is uniquely prescribed by its amplitude spectrum, the source/Green’s function deconvolution problem can be simplified to the isolation of their respective amplitude spectra. Our initial approach to addressing this problem relied on the use of the cross spectrum of two seismograms sharing a common component, e.g., two components of the same three-component seismogram. This approach was motivated by the recognition that the source, which represents that convolutional element common to both seismograms, is represented within the zero-phase component of the cross spectrum. Consequently, if we were able to extract the zero-phase component of the cross-spectrum we would have an estimate of the source spectrum that is accurate to the extent that no common (i.e. zero-phase) signal exists in the underlying Green’s functions correlogram. In theory, as explained by Baig et al. (2005), there are methods for achieving this zero-phase/non-zero-phase component separation. In practice, we have found that noise, coupled with uncertainties in the length of effective support interval occupied by the Green’s function cross-correlogram, render those approaches rather less effective than an alternative blind deconvolution method that we describe below.  2.2.1  Blind Deconvolution  In the approach adopted here we shall consider a three-component teleseismicP recording as resulting from the convolution of a band-limited source with respective components of a Green’s function that possess statistically white amplitude spectra. The motivation for this breakdown stems from the absorptive properties of the mantle which tend to attenuate higher frequencies of the teleseismic wavefield, and from our expectation that the receiver-side impulse response of the crust and lithospheric mantle is full-band. Consider, for example, the cross spectrum of two components of a single threecomponent seismogram, P (ω)V ∗ (ω) = |S(ω)|2 GP (ω)GV ∗ (ω),  (2.1) 9  Chapter 2. Improved Green’s functions for passive source structural studies where P (ω), V (ω) are the P- and SV-components of the seismogram in the frequency domain, GP (ω) and GV (ω) are the corresponding components of the Green’s function, |S(ω)|2 is the power spectrum of the source, and  ∗  denotes complex conjugation. The transformation from particle displacement into estimates of the up-going P-, SV- and SH-components of motion can be accomplished by coordinate rotation (Vinnik , 1977) or by assuming propagation through a laterally homogeneous Earth model (Kennett, 1991). The accuracy of this transformation is not important except insofar as energy in the incident P-wave is effectively isolated to the P-component. This isolation improves the likelihood that the minimum-phase assumption will be honored for this component; we may restore the recovered Green’s function to standard Cartesian displacement through the corresponding inverse transformation after deconvolution has been completed. The simplest means of estimating the source spectrum might be to assume that the amplitude spectrum of the individual Green’s function components is constant in frequency and, thus, that the source power spectrum, |S(ω)|2 , is equivalent to the amplitude spectrum of the seismogram cross-correlogram, |P (ω)V ∗ (ω)|. In synthetic tests comprising simple, horizontally layered models (Baig et al., 2005), we have gauged the accuracy of this approach by examining the quality of reproduction of the recovered Green’s function cross-correlogram (i.e. the inverse Fourier transform of GP (ω)GV ∗ (ω)). This estimate does generally possess many of the features of the true cross-correlogram especially at negative lags which are dominated by interaction of the incident P-pulse with the first-order Pto-S-conversions. This observation is perhaps not surprising given that it is the quantity P ∗ (ω)V (ω) that appears in the numerator of a standard, least-squares, frequency-domain, receiver-function deconvolution. The reproduction at positive lags is poorer as the signal levels are lower because the dominant, incident-P arrival does not contribute to this portion of the cross-correlogram. Lowest-order scattering interactions are identifiable but so are a small number of spurious artificial arrivals (see Figure 2.2) that could undermine confidence in geological interpretations. A less extreme approach is to assume that the Green’s function com10  Chapter 2. Improved Green’s functions for passive source structural studies  (a)  −4  −3  −2  −1  0  1  2  3  4  −3  −2  −1  0  1  2  3  4  −3  −2  −1  0 Time [s]  1  2  3  4  (b)  −4  (c)  −4  Figure 2.2: Comparison of original cross-correlograms, (a) True F −1 {GP (ω)GV ∗ (ω)} cross-correlogram and in (b) and (c) estimated F −1 {GP (ω)GV ∗ (ω)} cross-correlograms using flat spectrum and spectral smoothing, respectively. Note the good correspondance at negative lags for both estimates and a better correspondance for spectral smoothing at positive lags.  11  Chapter 2. Improved Green’s functions for passive source structural studies ponents possess spectra that are statistically white (as opposed to strictly constant), whereas the source power spectrum is band-limited, and then to estimate the latter quantity by spectral smoothing (Claerbout, 1992). In our experiments, we have found that a simple, time-domain low-pass Gaus2 2 sian filter (e(t /tc ) , where tc is the “cut-off time”) applied to the amplitude spectrum of the cross-correlogram, works as well as any other for smoothing purposes. The choice of the cut-off time, tc , influences the quality of the deconvolution. In the limit that tc extends to the full length of the signal, the Green’s function spectrum becomes constant and then, the recovered Green’s function cross-correlogram (F −1 {|GP (ω)GV ∗ (ω)|eφ(P (ω)V F −1  denotes inverse Fourier transformation and φ(P (ω)V  ∗ (ω))  ∗ (ω))  }, where  is the phase  of the original cross-correlogram) approaches that for a constant amplitude spectrum replete with spurious arrivals as discussed above. As tc decreases, the amount of smoothing increases, and the spurious arrivals diminish in size leaving the lowest order scattering interactions more clearly identifiable above ambient noise levels. At very short tc , the source spectrum whitens and the recovered Green’s function correlogram approaches that of the raw seismograms. From experiment, we have found that the choice of optimal tc corresponds to inclusion of the main lobe of the zero-phase quantity F −1 {|GP (ω)GV ∗ (ω)|} with secondary lobes strongly attenuated. The effect of the length of tc is shown in Figure 2.3  2.2.2  Multichannel Deconvolution  In order to further improve our estimates of the source and teleseismic Green’s functions spectra, we shall exploit the availability of multiple recordings and employ the multichannel spectral deconvolution approach originally presented by Andrews (1986). We shall consider a collection of threecomponent seismograms sharing the same impulse responses (that is the same source-receiver geometry) but characterized by different source-time functions. This approach contributes to increased redundancy and an improvement in the quality of the Green’s function spectral reconstruction. For a single station and a set of M independent sources we then have  12  Chapter 2. Improved Green’s functions for passive source structural studies  Seismogram Spectrum  Optimal tc  Long tc  Green’s Function Spectrum Estimate  Source Spectrum Estimate  Short tc  Figure 2.3: Effect of the “cut off time”, or smoothing parameter, on the source and Green’s function estimate.  13  Chapter 2. Improved Green’s functions for passive source structural studies 3 × M convolution equations which can be written for amplitude spectra in the frequency domain as Pm (ω) = |Sm (ω)|G(ω),  (2.2)  where Pm (ω) = [|Pm (ω)|, |Vm (ω)|, |Hm (ω)|]T is the mth three-component seismogram, G(ω) = |GP (ω)|, |GV (ω)|, |GH (ω)| function, and Sm (ω) the  mth  T  the three-component Green’s  source.  Transformation of the system in equation 2.2 into the log-spectral domain allows us to linearize it as log{Pm (ω)} = log{Sm (ω)}1 + log{G(ω)},  (2.3)  where 1 = [1, 1, 1]T . These equations are then augmented with source estimates, S˜m , obtained by the approach of the previous section, to yield a highly overdetermined system:           I3 1 · · · 0 .. .  1 ··· .. . . . .  I3 0 · · · 0  0 ···  0          1   1 0 .. .   log{P1 (ω)}   log{G(ω)} ˜     log{ S1 (ω) }  log{|S1 (ω)|}    .. = , . ..     .    log{PM (ω)}    log{|SM (ω)|} log{ S˜M (ω) }       (2.4)  where I3 is the 3 × 3 identity matrix. The solution of equation 2.4 recovers an improved estimate of the threecomponent teleseismic Green’s functions spectrum (see Baig et al. (2005)). By invoking the minimum-phase property of the P impulse response, one can reconstruct the phase of that signal and hence completely characterize it in the time domain. The phase spectra of the SV-and SH-components of the Green’s function are recovered from the estimates of their respective amplitude spectra, |GV (ω)| and |GH (ω)|, by applying a set of sequential all-pass filters. The construction of these filters can be accomplished via Kolmogorov  14  Chapter 2. Improved Green’s functions for passive source structural studies spectral factorization (e.g. Claerbout (1992)) of the initial, unprocessed P, SV, and SH seismograms. Defining K{·} as an operator which transforms a time series to minimum-phase and exploiting the inherent minimumphase nature of P-component Green’s function, GP (ω), we may write the minimum-phase-transformed P-component seismogram as K{Pm (ω)} = K{Sm (ω)}GP (ω);  (2.5)  that is, we generate the seismogram corresponding to a minimum-phase source, K{Sm (ω)}, with the same amplitude spectrum as Sm (ω). Accordingly, we may therefore define an all-pass filter that rearranges the phase of source m to minimum phase as S  eiφm (ω) =  K{Sm (ω)} K{Pm (ω)} = . Pm (ω) Sm (ω)  (2.6)  Application of this filter to the SV and SH seismograms, V (ω), H(ω), will likewise transform the source on these components to minimum-phase. One subsequent minimum-phase transformation of the resulting time series alV  lows us to define two more all-pass filters, eiφ  (ω)  H (ω)  and eiφ  the non-minimum phase Green’s function components,  , that transform  GV (ω)  and GH (ω),  to their minimum-phase counterparts, that is: S  V  eiφ  (ω)  =  K{eiφm (ω) Vm (ω)} K{K{Sm (ω)}GV (ω)} K{GV (ω)} = = , (2.7) K{Sm (ω)}GV (ω) GV (ω) eiφSm (ω) Vm (ω)  and similarly, S  H (ω)  eiφ V  Note that eiφ  (ω)  =  K{GH (ω)} K{eiφm (ω) Hm (ω)} = . GH (ω) eiφSm (ω) Hm (ω)  H (ω)  and eiφ  (2.8)  must be independent of source m as they are  ultimately defined solely in terms of the Green’s function. It is then straightforward to use them to reconstruct the phase of the quantities |GV (ω)| and |GH (ω)|, obtained from solution of equation 2.4 by applying their inverses  15  Chapter 2. Improved Green’s functions for passive source structural studies to the minimum-phase versions of their respective amplitude spectra: V  GV (ω) = e−iφ  (ω)  K{|GV (ω)|} and,  H (ω)  GH (ω) = e−iφ  (2.9)  K{|GH (ω)|}.  (2.10)  In practice, noise will render estimates of eiφV (ω) and eiφH (ω) obtained from different seismograms inconsistent and so averaging them is likely to improve the accuracy of the filters. Since we are ultimately interested in using the reciprocal quantities to restore the proper phase to the components GV (ω) and GH (ω) in equations 2.9 and 2.10 we construct the averages as iφV (ω)  V  A (ω)e  iφH (ω)  H  A (ω)e V  1 = M  −1  M −iφV m (ω)  e  and,  (2.11)  m=1  1 = M  −1  M −iφH m (ω)  e  ,  (2.12)  m=1  H  where eiφm (ω) and eiφm (ω) represent all-pass filters derived from the mth seismogram. The all-pass nature of the filters is now lost, however, it is easily shown that this choice of average when applied to equations 2.9 and 2.10 leads to a downweighting at frequencies where estimates are least consistent (i.e. AV (ω), AH (ω) > 1), that is, most affected by noise.  2.3  Results  In this section, we demonstrate how the deconvolution method previously presented can be applied to real three-component teleseismic data to successfully recover three-component Green’s functions. We will focus on, but not limit our attention to, the continental crust-mantle boundary or “Moho” and present results from three stations of the Canadian National Seismograph Network (CNSN), namely ULM, WHY, and EDM, that clearly reveal the pure P Moho multiple P pM p. Before discussing results from real data, we shall, however summarize the properties of the principal scattered phases 16  Chapter 2. Improved Green’s functions for passive source structural studies generated by an incident plane P-wave within an idealized isotropic layerover-a-half-space crustal model. Figure 2.1a illustrates the ray-paths of the first-order scattered phases that result from interaction between the incident wavefield and the crustmantle boundary. These phases are all scattered once from the Moho and exist on either P- or SV-components of the teleseismic-P Green’s function, depending on the nature of their final leg. Figures 2.1b and 2.1c present typical P- and SV-component Green’s function seismograms with the same scattered arrivals identified for the simple crustal model whose physical properties are provided in the caption. The magnitude and the polarity of the different phases depend both on the reflection and transmission coefficients at the Moho and free surface. The polarities of P , PM s and P pM s are positive whereas those of P pM p, P sM p and P sM s are negative. The amplitudes of the different first-order scattered phases, excluding P sM p, are comparable and represent approximately 5% of the incident P amplitude for this value of horizontal slowness or ray parameter. The amplitude of P sM p is somewhat smaller at near-vertical incidence, since this multiple is converted twice en route to the receiver, unlike the other phases that experience at most one conversion. We note once more that P pM p is the only phase whose amplitude is sensitive to the contrast in compressional modulus (more specifically, the P-impedance) at the Moho; all the other phases shown in Figure 2.1 are sensitive only to combinations of S-velocity and density. Both amplitude and timing of the different phases vary with horizontal slowness (ray parameter). For the narrow range in slownesses represented by teleseismic P (0.04 − 0.08 s km−1 ), the utility of amplitude versus slowness analysis will be limited. It is generally possible, however, to trace the move-out of a given arrival over this slowness band as a useful means of identifying it as either a forward-scattered direct conversion or back-scattered multiple. Figure 2.4 shows a plot of arrival time, relative to P , of the direct PM s conversion and the first-order Moho multiples for the simple crustal model described in Figure 2.1. We note that the PM s traveltime increases with slowness while those of the other first-order back-scattered 17  Chapter 2. Improved Green’s functions for passive source structural studies 0  P-component only SV-component only P- and SV-components  2 4 PM s  Time [s]  6 8 10 P pM p 12 14 P sM p + P pM s 16 18 P s s M 20 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km]  0.07  0.075  0.08  Figure 2.4: Timing relative to incident P of the first-order crustal scattered phases as a function of slowness for the model described in Figure 2.1. phases decrease. For typical continental crust (as described in Figure 2.1), PM s arrives between 4 and 5 s after direct P followed by P pM p near 10 s on the P component. The kinematic analogues, P sM p and P pM s follow near 14 s on separate components with the final first-order scattered phase P sM s coming in close to 19 s.  2.3.1  Station ULM (Lac du Bonnet, Manitoba)  We now direct attention to the first of CNSN stations to be discussed. Station ULM is located at Lac du Bonnet, Manitoba in the western Superior province of the Canadian Shield. From previous work (Bank and Bostock , 2003), we know that the site is characterized by a relatively simple, approximately 1D crustal structure. For our analysis, we used seismograms from 344 events representing a wide range of back-azimuths. This data set was sorted into 150 equally spaced bins in horizontal slowness between 0.04 to 0.08 s km−1 based on the IASPEI reference model (Kennett and Engdahl , 1991). Thus, an average of just over two seismograms contributes to  18  Chapter 2. Improved Green’s functions for passive source structural studies each bin. All seismograms falling within a given slowness bin were then simultaneously deconvolved to produce Green’s function estimates for the corresponding slowness. Each component of the Green’s function is filtered between 0.01 and 1.0 Hz and plotted in a 2D image in Figures 2.5a and 2.6a where the principal phases have been identified according to their timing and move-out. Note that the P-component image is muted with a cosine taper between 0 and 1.0 s before filtering to remove the effect of the dominant, incident P arrival. The SV-component Green’s function presented in Figure 2.5a shows coherent signals with timing and move-out corresponding to the three expected, first-order scattered phases from the Moho (PM s, P pM s, P sM s). Figure 2.5b shows the result of stacking these time series along the PM s, P pM s and P sM s move-out curves defined in Figure 2.4. These plots further support the interpretation of these scattered Moho phases with the correct polarization and arrival time, although, their relative amplitudes diverge somewhat from expectation. P pM s and P sM s magnitudes are significantly smaller than PM s and approximately one third those predicted from synthetics obtained using an isotropic, non-dissipative model. This discrepancy may be due to energy loss through incoherent scattering and/or intrinsic absorption, or possibly to an extended Moho transition which scatters forward propagating energy more efficiently than back scattered energy (Bostock , 1999). We note that a receiver function image would appear very similar to that in Figure 2.5a since the receiver function is a first-order approximation to the SV-component Green’s function. The P-component image in Figure 2.6a also reveals a coherent signal at around 10 s following the P arrival with negative polarity relative to direct P . It appears stronger at steep angles of incidence and tends to fade as slowness increases in agreement with theoretical expressions for the reflection coefficient presented by Aki and Richards (1980). The assumption of 1D structure could also contribute to the diminished signal at larger slowness values. According to its timing and move-out it can be identified with the receiver-side P pM p multiple. This phase is also clearly evident at the expected arrival time on stacked plots of the Figure 2.6b where a strong 19  Chapter 2. Improved Green’s functions for passive source structural studies  2 4  PM s  Time [s]  6 8 10 12 14  P pM s  16 18  P sM s  20 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km]  0.07  0.075  (a)  0.1 0  Amplitude Relative to P  PM s −0.1 0.1 0  P pM s −0.1 0.1  P sM s 0 −0.1 0  5  10  15 Time [s]  20  25  30  (b)  Figure 2.5: (a) Image of the SV-component Green’s function as a function of slowness for station ULM and (b) stack along move-out curves of PM s, P pM s and P sM s.  20  Chapter 2. Improved Green’s functions for passive source structural studies  2 4  Time [s]  6 8 10  P pM p  12 14 16 18 20 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km]  0.07  0.075  (a)  0.1  Amplitude Relative to P  0.05  P pM p  0 −0.05 −0.1  0.1 0.05  P sM p?  0 −0.05 −0.1 5  10  15  20  25  30  Time [s] (b)  Figure 2.6: (a) Image of the P-component Green’s function as a function of slowness for station ULM and (b) stack along move-out curves of P pM p and P sM p.  21  Chapter 2. Improved Green’s functions for passive source structural studies negative peak is observed. In contrast, the weaker, reflected-converted P sM p phase is not so clearly defined neither on the 2-D image nor on Figure 2.6b. To the best of our knowledge, Figure 2.6a,b represents the first unambiguous identification of the P pM p phase.  2.3.2  Station WHY (Whitehorse, Yukon)  Station WHY is situated near Whitehorse, Yukon, within the Cenozoic Canadian Cordilleran orogen (Gabrielse and Yorath, 1992). The images for WHY presented in Figure 2.7a and 2.7b were produced using 415 events, again covering a large range of back azimuths, that have been sorted into 150 equally spaced slowness bins and processed similarly to the ULM data set. Unlike ULM, the Moho at this station shows evidence of a more complex crust-mantle boundary. On the S-component image, the PM s phase is found between 4 s and 5 s with a well defined move-out and polarization. The main peak is followed immediately by a second, negative peak. A comparable signal with opposite move-out and smaller amplitude, interpreted to be the P pM s phase, is observed at 14 s showing a similar waveform. P sM s, which is expected at ∼ 18 s, is not readily apparent on the image. On the P-component image, P pM p is clearly defined at 10 s, and like PM s and P pM s, is composed of two distinct pulses of opposite polarity. This bipolar PM s signal was previously observed by Lowe and Cassidy (1995) using a reduced data set. We note that it suggests the presence of a high P- and S-velocity layer immediately underlying the Moho.  2.3.3  Station EDM (Edmonton, Alberta)  Station EDM is located near the town of Leduc roughly 25 km to the southwest of Edmonton, Alberta in the Western Canada Sedimentary Basin. The sedimentary platform in this region is known to be approximately 2.5 km thick (Cassidy, 1995). Our analysis for this station is based on 307 seismograms distributed among 115 equally spaced slowness bins. The data set is again processed using the parameters described for station ULM. 22  Chapter 2. Improved Green’s functions for passive source structural studies  Time [s]  5  PM s  10  15  P pM s  20  P sM s  25 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km]  0.07  0.075  0.07  0.075  (a)  Time [s]  5  10  P pM p  15  20  25 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km] (b)  Figure 2.7: Image of the (a) SV-component, and (b) P-component Green’s function as a function of slowness for station WHY.  23  Chapter 2. Improved Green’s functions for passive source structural studies The P- and SV-component images shown in Figure 2.8a and 2.8b display the usual first-order scattered Moho phases described in preceding sections, albeit the back-scattered multiples are rather weak. In addition, there are several signals that clearly originate from shallower crustal levels. As discussed in previous receiver function studies (Cassidy, 1995; Eaton and Cassidy, 1996), these signals are due to the basin sediments overlying a crystalline crust that exhibits a pronounced low-velocity zone between 12.517 km depth. The sedimentary package (denoted by “S”) is evident as coherent energy immediately following direct P on both components. On the SV-component, PS s energy associated with the sediments dominates the Green’s function for 1.5 s, while its duration on the P-component as P pS p extends to 3 s as a result of the longer path length for the multiple. Over this frequency band it is difficult to discern details of basin structure but both signals are likely dominated by the scattering from the sediment-basement contact. The low-velocity zone (denoted by “L”) appears as the 2 or 3-lobed signal P pL p of opposite polarity on the P image near 4 s and PL s on the SV images at approximately 2 s. Eaton and Cassidy (1996) combined receiver function observations with Lithoprobe crustal-scale reflection data to interpret the layering as serpentinized former oceanic lithosphere thrust up into mid-crustal levels during an episode of Proterozoic subduction. Our observation of a coherent P pL p thus bridges a gap in their observations as pure P reflections observed at the frequencies of conventional receiver functions.  2.4  Discussion and Concluding Remarks  In this paper we have described a new source deconvolution technique which exploits the minimum-phase property of the teleseismic-P Green’s function to enable recovery not only of S components but also the P component of this quantity. We have demonstrated that the approach can be applied to real teleseismic data to extract more accurate estimates of the impulse response than possible using the conventional receiver function approach. In 24  Chapter 2. Improved Green’s functions for passive source structural studies  PS s PL s PM s  Time [s]  5  10  15 P pM s  20  P sM s  25 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km]  0.07  0.075  0.07  0.075  (a)  P pS p P pL p  Time [s]  5  10 P pM p  15  20  25 0.04  0.045  0.05  0.055 0.06 0.065 Slowness [s/km] (b)  Figure 2.8: Image of the (a) SV-component, and (b) P-component Green’s function as a function of slowness for station EDM. Note the presence of the intra-crustal phases (i.e. PS s, PL s, P pS p and P pL p).  25  Chapter 2. Improved Green’s functions for passive source structural studies particular, we have for the first time, to the best of our knowledge, unambiguously identified and documented the receiver-side P pM p phase from teleseismic observations. The amplitude of this phase provides a measure of the compressional impedance contrast across the crust-mantle boundary which is not afforded by other signals in the teleseismic P wavetrain. We note that previous work on extracting P-component information from teleseismic P (Bostock and VanDecar , 1994; Revenaugh, 1995; Li and Nabelek , 1999; Langston and Hammer , 2001) has generally involved stacking seismograms recorded at many stations to produce an approximate source-time function under the assumption that the P-coda is incoherent across stations. The main problem with this latter approach is that scattered signals recorded across an array that exhibit little temporal moveout (e.g. the P pM p phase in many instances) are incorrectly mapped to the source and are thus absent from the resulting Green’s function. As we have demonstrated here, our approach does not suffer from this shortcoming. Results have been presented for three stations of the Canadian National Seismic Network, each representing a different geologic environment. Each station, namely ULM (Lac du Bonnet, Manitoba), WHY (Whitehorse, Yukon) and EDM (Edmonton, Alberta), shows evidence of the principal, first-order scattered Moho phases and, in particular, the P pM p multiple. Station ULM exhibits a very simple crustal structure that is reasonably well modelled as a layer over a half-space. The crust-mantle boundary at station WHY, as portrayed through PM s and P pM p, is revealed to be rather more complex, representing a thin, high-velocity lid. Intra-crustal scattering dominates the Green’s functions recovered at station EDM. The sedimentarybasement contact produces prominent arrivals at early times followed by signals from a pronounced low-velocity layer that has been previously interpreted from seismic reflection and receiver functions to represent a serpentinized sliver upthrust during Proterozoic subduction.  26  Bibliography Aki, K., and P. G. Richards, Quantitative seismology, Theory and Methods: Volume 1, W. H. Freeman and Company, 1980. Andrews, D. J., Objective determination of source parameters and similarity of earthquakes of different sizes, 37 ed., 259-268 pp., American Geophysical Union, 1986. Baig, A. M., M. G. Bostock, and J.-P. Mercier, Spectral reconstruction of teleseismic-P Green’s function, J. Geophys. Res., 110, doi:doi:10.1029/ 2005JB003625, 2005. Bank, C.-G., and M. G. Bostock, Linearized inverse scattering of teleseismic waves for anisotropic crust and mantle structure: 2. Numerical examples and application to data from Canadian stations, Journal of Geophysical Research, 108, 10.1029/2002JB001,950, 2003. Bostock, M. G., Seismic waves converted from velocity gradient anomalies in the Earth’s upper mantle, Geophysical Journal International, 138, 747– 756, 1999. Bostock, M. G., Green’s functions, source signatures, and the normalization of teleseismic wave fields, Journal of Geophysical Research, 109, 10.1029/2003JB002,783, 2004. Bostock, M. G., and J. C. VanDecar, The influence of crust and upper mantle heterogeneity on short period waveform distortion, Physics of the Earth and Planetary Interiors, 83, 225–247, 1994. Bostock, M. G., S. Rondenay, and J. Shragge, Multiparameter twodimensional inversion of scattered teleseismic body waves: 1. Theory for oblique incidence, Journal of Geophysical Research, 106, 30,771–30,782, 2001. Burdick, L. J., and C. A. Langston, Modeling crustal structure through the use of converted phases in teleseismic body-wave forms, Bulletin of the Seismological Society of America, 67, 677–691, 1977. 27  Bibliography Cassidy, J. F., A comparison of the receiver structure beneath stations of the Canadian National Seismograph Network, Canadian Journal of Earth Sciences, 32, 938–951, 1995. Claerbout, J., Earth Sounding Analysis: Processing Versus Inversion, Blackwell, Oxford, 1992. Dueker, K. G., and A. F. Sheehan, Mantle discontinuity structure from midpoint stacks of converted P and S waves across the Yellowstone hotspot track, Journal of Geophysical Research, 102, 8313–8326, 1997. Eaton, W. E., and J. F. Cassidy, A relic Proterozoic subduction zone in western Canada: New evidence from seismic reflection and receiver function data, Geophysical Research letters, 23, 3191–3794, 1996. Gabrielse, H., and C. Yorath (Eds.), Geology of the Cordilleran Orogen in Canada, vol. 4, 847 pp., Geological survey of Canada, 1992. Kennett, B. L. N., The removal of free surface interactions from threecomponent seismograms, Geophysical Journal International, 104, 153–163, 1991. Kennett, B. L. N., and E. R. Engdahl, Traveltimes for global earthquake location and phase identification, Geophysical Journal International, 105, 429–465, 1991. Langston, C. A., Structure under mount Rainier, Washington, inferred from teleseismic body waves, Journal of Geophysical Research, 84, 4749– 4762, 1979. Langston, C. A., and J. K. Hammer, The vertical component P-wave receiver function, Bulletin of the Seismological Society of America, 91, 1805– 1819, 2001. Li, X.-Q., and J. L. Nabelek, Deconvolution of teleseismic body waves for enhancing structure beneath a seismometer array, Bulletin of the Seismological Society of America, 89, 190–201, 1999. Lowe, C., and J. F. Cassidy, Geophysical evidence for crustal thickness variations between the Denali and Tintina fault systems in west-central Yukon, Tectonics, 14, 909–917, 1995. Owens, T. J., and G. Zandt, The response of the continental crust-mantle boundary observed on broadband teleseismic receiver functions, Geophysical Research Letters, 12, 705–708, 1985. 28  Bibliography Phinney, R. A., Structure of the Earth’s crust from spectral behavior of long-period body waves, Journal of Geophysical Research, 69, 2997–3017, 1964. Poppeliers, C., and G. L. Pavlis, Three-dimensional, prestack, plane wave migration of teleseismic P-to-S converted phases, Journal of Geophysical Research, 108, 10.1029/2001JB000,216, 2003. Revenaugh, J. A., Scattered-wave image of subduction beneath the Transverse Ranges, California, Science, 268, 1888–1892, 1995. Sherwood, J. W. C., and A. W. Trorey, Minimum-phase and related properties of a horizontally stratified absorptive Earth to plane acoustic waves, Geophysics, 30, 191–197, 1965. Vinnik, L. P., Detection of waves converted from P to SV in the mantle, Physics of the Earth and Planetary Interiors, 67, 39–45, 1977. Zandt, G., and J. A. 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, Journal of Geophysical Research, 105, 2969–2980, 2000.  29  Chapter 3  The teleseismic signature of fossil subduction: Northwestern Canada  3.1  Introduction  Northwestern Canada is composed of a variety of geological terranes that formed and accreted over the past 4.0 Ga. This unique collage spans over 5000 km from the Pacific coast in the west to the Slave craton in the east, and represents the most nearly complete and contiguous sampling of geological time on the surface of the Earth. Consequently, the region is particularly well suited to address fundamental questions concerning the nature of continental evolution and variations in structure and geometry of the subsurface from early Archean (3.5-4.0 Ga) to present. In the 1990’s, the Lithoprobe Slave-Northern Cordillera Lithospheric Evolution transect, hereafter referred to as SNORCLE, was undertaken over this region (Cook and Erdmer , 2005). It comprised a series of geophysical and geological experiments that led to several important discoveries. In particular, an active-source, seismic reflection experiment in the Wopmay orogen revealed the presence of several clearly defined and laterally coherA version of this chapter has been published in the Journal of Geophysical Research in 2008. Mercier, J.-P., M. G. Bostock, P. Audet, J. B. Gaherty, E. J. Garnero, and J. Revenaugh (2008), The teleseismic signature of fossil subduction: Northwestern Canada, J. Geophys. Res., 113, doi:10.1029/2007JB005127. c 2008 American Geophysical Union.  30  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada ent mantle reflectors (Cook et al., 1998, 1999). The most prominent among them is an eastward dipping layer in the lithospheric mantle merging with the crust at the suture of Fort Simpson and Hottah terranes, and extending eastward beneath the Great Bear magmatic arc where it reaches a depth of ∼100 km (see Figure 3.1b for geographic location). Geological and geophysical evidence indicates that this feature is a remnant of Proterozoic subduction. Although several other similar reflectors were observed further east toward and below the Slave province, their relation with the dipping structure described above has remained unclear. Also beneath the Slave province, a teleseismic-P receiver function study by Bostock (1998) using data recorded at the Yellowknife (YKA) mediumaperture seismic array shows a well developed mantle stratigraphy. Three major layered structures, defined by discontinuous boundaries between azimuthally anisotropic domains, were identified, and labeled H, X and L at depths of 75, 135, and 170 km, respectively. Labels H and L are used generically to represent the so-called Hales and Lehman discontinuities, see also, e.g., Revenaugh and Jordan (1991) and Gaherty and Jordan (1995). The structures beneath the Slave craton were hypothesized to represent stranded, former oceanic crustal material emplaced during a series of subduction episodes. This discovery prompted speculation on a possible continuation of the dipping layer, defined on the SNORCLE reflection profile, eastward beneath the Slave province, and on the importance of subduction processes in the formation and stabilization of the underlying cratonic root. Our study presents teleseismic receiver-function images of the mantle below the Wopmay orogen from 20 three-component, broadband seismic stations installed in the Northwest Territories as part of the joint Incorporated Research Institutions for Seismology (IRIS)-Lithoprobe CAnadian NOrthwest Experiment (hereafter referred to as CANOE) which complements the earlier SNORCLE near-vertical reflection (Cook et al., 1998) and refraction/wide-angle reflection (Fern´ andez-Viejo and Clowes, 2003) profiles. In what follows, we provide further constraints on the nature and geometry of the subsurface structure mentioned above which yield insight into the importance of subduction processes in the stabilization of the Slave 31  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  230˚  Wopm ay  240˚  65˚  ˚ 260 70˚  250˚  Thelon  220˚  Slave  70˚  (a)  65˚  Yellowknife  Leg A Leg B  Fort Nelson  Proterozoic  Leg C Edmonton  Phanerozoic  50˚  55˚  Hudso  Archean  n  Cordil lera  55˚  60˚  Trans−  Whitehorse  Hearne  60˚  50˚ 260  220˚  250˚  230˚  (b) 237˚ 62˚30'  238˚  239˚  0  240˚  240˚  241˚  242˚  ˚  km 500  243˚  244˚ 62˚30'  Great bear magmatic arc A17  Hottah Terrane  Fort Simpson Terrane  62˚00'  62˚00'  Fort Simpson  U  61˚00'  C  09  UBC10 A15  B  Nahanni Terrane  61˚30' Fort Providence  U  61˚00'  A16  A U 10 B C 01 A BC 11 02 U B C 03 A U 12 B C 0 U 4 B C U 05 B C A 06 1 U 3 B U C0 B 7 C A 0 14 8  A  08  A  09  61˚30'  km 60˚30' 237˚  238˚  239˚  240˚  241˚  242˚  0 243˚  50  60˚30' 244˚  Figure 3.1: (a) Map of western Canadian geological provinces illustrating the distribution of the broad-band three-component CANOE stations along legs A, B and C. The three legs radiate away from Fort Nelson to Yellowknife (leg A), Whitehorse (leg B) and Edmonton (leg C). The study area is identified by the inset portion of leg A. (b) Inset of study region where 20 broad-band three-component seismic stations were deployed at ∼35 km interval at the ends and ∼15 km in the center. The dashed lines delineate different geological terranes within the Wopmay Orogen.  32  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada craton.  3.2  Tectonic  The area of interest in this study spans, from west to east, the entire Wopmay orogen complex and the western portion of the Slave province, known as the Anton terrane (see Figure 3.1b). The Slave province is a relatively small craton that forms, together with the Nain, Superior and Rae provinces, one of the nuclei of the Canadian Shield, the core of the North American continent. It contains the oldest dated rocks on Earth, the Acasta gneiss, estimated at more than 4.03 Ga in age (Bowring and Williams, 1999). Its assembly was complete by 2.6 Ga, upon the accretion of volcanic arcs and micro-continents (Isachsen and Bowring, 1994). The formation of Wopmay orogen, a Paleoproterozoic collage of domains, occurred between 2.1 and 1.84 Ga (Hildebrand et al., 1987). It is composed of three distinct elements: the Hottah, Fort Simpson and Nahanni terranes. The Hottah terrane developed between 1.92 and 1.91 Ga as a magmatic arc. It collided with the Slave province during the Calderian Orogeny (1.901.88 Ga), shortening and displacing the sedimentary rocks of the Coronation margin. Subduction of an oceanic plate beneath the Hottah terrane between 1.88 and 1.84 Ga resulted in formation of the Great Bear magmatic arc located at the eastern edge of the Hottah terrane (Hoffman and Bowring, 1984). The Fort Simpson terrane accreted to the western margin of Hottah between 1.845 and 1.745 Ga (Villeneuve et al., 1991). The Nahanni terrane is located west of the Fort Simpson terrane and is the underlying basement to the Northern Cordillera. It reveals no surface exposure and remains an enigmatic geological feature.  3.3  Data  The CANOE project involved the deployment of an expansive seismic array comprising nearly 60 three-component broadband stations (see Figure 3.1a). Instrumentation was a mix of Guralp CMG-3T, CMG-ESP, and CMG-40T 33  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada seismometers and RefTek recorders, provided by various institutions (see Acknowledgements). The array consisted of 3 legs which radiate outward from the array center at Fort Nelson, to Yellowknife (leg A), Whitehorse (leg B) and Edmonton (leg C) spanning more than 3500 km in aperture and traversing an extensive suite of geological domains. The instruments were deployed in June 2003 and 2004 during two separate campaigns and were decommissioned in September 2005 after 27 and 15 months, respectively. The station spacing is on average 35 km over most of the array with the exception of a stretch of Mackenzie Highway between Fort Simpson and Fort Providence on leg A where spacing was reduced to ∼15 km. In the present study, we employ a subset of the data from 20 CANOE stations (A08-A17 and UBC01-UBC10) (Figure 3.1b) traversing the entire Wopmay orogen complex and the western edge of the Slave craton. Over the duration of the experiment teleseismic-P waves from ∼250 earthquakes larger than magnitude 5.5 were recorded in the epicentral distance range 30◦ -100◦ , mainly from the Western Pacific, Fiji-Tonga, Central and South American regions. Consequently, the back-azimuthal coverage is good from 225◦ to 40◦ and from 125◦ to 175◦ and less regularly sampled otherwise (see Figure 3.2). Seismogram selection based on signal-to-noise ratio (SNR) (threshold set to 15 dB for the P -component seismogram in the 0-2 Hz band) followed by a visual inspection, resulted in a data set that comprises 1428 fair-to-good quality, broadband, three-component recordings with an average of 71 seismograms per station.  3.4  Methodology  The isolation of scattered S -waves in the teleseismic-P coda produced through interaction of the near-vertically propagating incident P -wave with subhorizontal discontinuities requires accurate removal of the earthquake source time function. This operation is conventionally performed through a spectral equalization procedure referred to as receiver function analysis (Vinnik , 1977; Langston, 1979). It consists of accurate decomposition of the incident wavefield into upgoing P and S modes (e.g. Kennett (1991)) which, when 34  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  30° 60° 90° 120°  Figure 3.2: Equidistant azimuthal projection centered on station UBC06 showing the distribution of events recorded over the period 2003-2005 with sufficiently high signal-to-noise ratio at one or more stations of the array to be used in this study. Note the relatively uniform azimuthal sampling between back-azimuth of 225◦ and 40◦ (events mainly from Pacific ring of fire, Asia, middle east and Northern Africa), and between 125◦ to 175◦ (events mainly from the Americas).  35  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada combined with a Wiener spectral deconvolution (Wiener , 1949) of the S components by the P -component, allows the recovery of an estimate of the S -components (radial and transverse) of the receiver-side Earth’s Green’s function. The timing of direct P -to-S conversions in the resulting time series depends on the velocity structure and is proportional to the depth of discontinuity whereas their amplitudes scale with the magnitude and sign of the contrasts in elastic moduli and density. Although a formal 2-D or 3-D inversion for anisotropic parameters is, in principle, possible (Burridge et al., 1998), teleseismic datasets are generally too poorly sampled in slowness and back-azimuth to afford a well-posed inverse problem. We have opted for a simpler imaging scheme that involves the projection of receiver functions to a 2-D profile wherein the amplitudes are color-coded/shaded and displayed as a function of time and horizontal distance. In this way the entire receiver function data set is directly represented and laterally coherent signals, even if weak, are readily apparent. Receiver functions are ordered from left to right along the resulting profile based on the geographic position of the stations and, for individual stations, on the sampling point of the P -S conversions projected onto the westeast direction (Figure 3.3a). The station spacing along leg A varies from ∼45 km (A08-A10 A15-A17) at the ends to ∼15 km in the centre (stations A10 to A15). Differences in station spacing are accounted for by a variable horizontal stretching proportional to the station separation (Figure 3.3b). This approach allows comparison of the profiles with time-migrated reflection data of Cook et al. (1999) through a simple scaling between the direct P -to-S conversion time and the 2-way vertical P reflection time as described in more detail in section 3.7.  3.5  Results  Figures 3.4a,b present the raw radial (SV) and transverse (SH) component receiver-function images obtained from the CANOE data. Figure 3.4c displays a version of the tranverse component wherein polarities have been corrected to account for the periodicity of the signal identified as AM. The 36  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  (a)  Surface west  east  Transverse Component Time [s]  (b)  7  6  A1  U  A1  BC U 09 BC 10 A1 5  08 4 A1  07 BC  BC U  U  06 3 A1  05 BC  BC U  U  A1 02 1 U BC 03 U BC A1 04 2  BC U  0 BC U  A1  9 A0  A0  8  01  Discontinuity at depth  0  5  10  15  0  50  100  150  200  250  300  350  Distance along the east/west direction [km]  Figure 3.3: (a) Illustration of the approach adopted to display receiver functions along the array. The receiver functions are ordered from west to east based on i) the geographic position of the stations, and ii) the azimuthal sampling of events for individual stations. (b) Display format adopted for Figures 3.4 and 3.11. Black/white shading corresponds to positive/negative amplitudes, respectively. Seismic stations location are represented by reverse triangles, and projected to their location on the 61st parallel. To produce a continuous, west-east profile that accounts for variable station spacing, receiver functions for individual stations are horizontally scaled to cover that segment of the image proportional to the local interstation distance (ie. half the sum of the distances separating a given station from its two nearest neighbours). Black vertical lines mark the plotting limits of the receiver functions for an individual station.  37  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada nature of this correction is described in more detail later in this section. For the sake of comparison with the SNORCLE reflection profile, we have plotted both time and approximate depth axes where, for the latter, we √ have assumed a Poisson solid (V p/V s = 3) with P-velocities of 6 km/s and 8 km/s for crust and mantle, respectively. Note that the assumption of constant Poisson’s ratio may not hold along the entire profile (see Fern´ andezViejo et al. (2005)), but resulting image degradation should be minor given the low frequency content of teleseismic data and the rudimentary imaging scheme. The conversion between delay time ∆tps and depth interval ∆z within each layer is thus given by: ∆z =  Vp ∆tps 3 − Vp2 p2 −  ,  (3.1)  1 − Vp2 p2  where p is horizontal slowness and Vp is the P-wave velocity. In the following sections we direct our attention to five prominent and laterally coherent signals labelled S, TM, AM, A1 and A2. Signal S is defined on the radial component by a succession of strong positive (red) and negative (blue) pulses within the first 2 s along most of the array, with the exception of stations A16 and A17. Since the occurrence of signal S coincides with the presence of a shallow ∼1 km thick layer of Phanerozoic sediments that extends from the western limit of our array to the vicinity of Fort Providence (see Figure 3.1b), it likely represents free-surface reverberations from this structure. The relatively large pulse amplitudes and duration preclude structural investigation of the crust to a depth of ∼15 km.  38  Radial Component Time [s] 15  10  5  0  15  10  5  0  15  10  5  0  0  50  TM  9 A0  100  150  AM  AM  U  05  BC U  BC  TM2  A1  07 08 3 BC BC 14 A1 U U A  200  06  Distance along the east/west direction [km]  TM1  S  p. 01 02 03 04 m Si 10 BC BC 11 BC BC 2 t F A U U A U U A1  250  A2  . 09 10 rov BC BC t P 15 U U F A  300  A1  6  140  120  100  80  60  40  20  0  140  120  100  80  60  40  20  0  140  120  100  80  60  40  20  0  7  A1  Estimated Depth [km] Estimated Depth [km]  39  Figure 3.4: Radial component in (a), transverse component in (b), and polarity corrected (stations A10 through A15 only) transverse component in (c). Red/blue colors correspond to positive/negative amplitudes, respectively. In (c), the correction is applied to account for the polarity cross-over in signal AM documented in section 3.6.3.5. Responses are filtered between 0.2 and 1.5 Hz.  Corrected Transverse Component Time [s]  Estimated Depth [km]  (c)  Uncorrected Transverse Component Time [s]  (b)  (a)  8 A0  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada On the same component, a clear positive (red) pulse labeled TM (Teleseismic Moho) is visible at around 4 s (∼30 km) beneath station A08 at the western end of the profile, and corresponds to the P -to-S conversion from the crust-mantle boundary (Moho). It shifts to later times beneath station A09 to ∼4.5 s, and is relatively constant in time along the rest of the profile with the exception of a 75 km long disruption beginning near kilometer 125. At this point, the signal appears to dip eastward into the mantle and can be traced to kilometer 175 where it reaches a depth of ∼50 km (TM1 ). A shallower, intermittent Moho signal (TM2 ) becomes evident a few 10’s of km to the east before becoming better defined and more coherent by kilometer 200. The transverse component (Figure 3.4b) is dominated by a high amplitude signal comprising at least two parallel, oppositely polarized pulses collectively labelled AM (Anisotropic Mantle). The signal appears first at 4 s beneath station A10, is flat for 75 km and begins to dip thereafter. Its geometry is well defined until kilometer 250 whereupon it becomes unclear. The dominant early pulse, whose timing corresponds to the direct conversion from the Moho on the radial component (i.e. signals TM and TM1 ), is sharp and rich in high frequencies whereas the second pulse is more diffuse. This latter signature is characteristic of an interface defined by a more gradual variation (e.g. first-order versus zeroth-order discontinuity) in material properties (Bostock , 1999). The polarity of AM varies as a function of back-azimuth, alternating from positive to negative (red to blue). Where the layer dips and its geometry is well defined (i.e. roughly between kilometers 150 and 220), a dominantly 1-θ (i.e. 360◦ ) periodicity is observed. This behaviour is illustrated in Figure 3.5 where data from stations UBC03, UBC04 and A12 have been aligned and sorted by back-azimuth. The polarity of the first arrival is clearly negative between 250◦ and 60◦ , weakly positive between 230 − 240◦ , and likely positive between 130◦ and 180◦ . Polarity reversals must evidently occur where no data are recorded between 60◦ and 130◦ , and between 240◦ and 250◦ . The precise back-azimuths at which reversals occur cannot be determined but, if the signal is strictly 1-θ (i.e. polarity reversals occurring 40  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada 0 2 4 Time [s]  AM 6 8 10  −  −  12  + 14 50  100  150 200 250 Back Azimuth [°]  300  350  Figure 3.5: Transverse component azimuthal response for stations UBC03, UBC04, and A12 in the frequency band 0.2-1.5 Hz. Receiver functions are aligned along signal AM and stacked into 100 back-azimuth bins. Transparent mask was applied to the signal of interest between 6.5 and 8.5 s. Vertical lines mark estimated back-azimuth of polarity cross-over. Expected 1-θ amplitude variation is shown in solid line. periodically every 180◦ ), they will be restricted to the range between 60◦ and 65◦ and between 240◦ and 245◦ . Visual coherence of the signal AM can be greatly improved if polarity corrections that account for the natural periodicity of the signal, are applied to the transverse component between stations A10 and A15 (see Figure 3.4c). Note that the correction derived from the polarity of the signal AM along the dipping portion is also somewhat successful in improving its coherence over the flat portion. On the corrected transverse component, additional energy that was not readily apparent reveals further signals A1 and A2. They appear in the mantle, above the layer defined by AM and between kilometers 200 and 300. Their combined signature resembles AM in that it exhibits a similar pulse sequence and azimuthal periodicity.  41  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  3.6  Modelling and Numerical Simulation  Estimation of elastic properties of the subsurface is an important step toward achieving an accurate interpretation of lithospheric structures from coherent signals on the receiver function images. In this section, we construct a model that succeeds in explaining the essential spatio-temporal distribution and frequency response of energy on the radial and transverse components. We draw our attention, in particular, to the layering defined by the signal AM on the transverse component. The appearance of coherent energy on the transverse-component receiver functions results from a rotation of the particle motion out of the sagittal plane (i.e. the vertical plane containing earthquake source and receiver). This rotation is the consequence of lateral heterogeneity, anisotropy or some combination thereof. In the case of AM, the signal clearly defines a dipping layer, so lateral heterogeneity must obviously contribute. There is a reason to believe, however, that anisotropy also influences the signal, in part, through the association of the continental Moho with the top of the layering AM. Note first that the absence of strong conversion from the lower (basal) interface of AM on the radial component suggests that average velocities of the layer and the underlying mantle are comparable. Moreover, the strength of AM on the transverse component does not diminish significantly as the layer dip shallows, as would be expected for a strictly isotropic layer.  3.6.1  Inversion  To develop a model that successfully reproduces the characteristic signature of AM, most clearly expressed at stations UBC03, UBC04 and A12, we generated a suite of synthetic seismograms for plane incident P -waves at a fixed slowness of 0.06 s/km and for back-azimuths varying from 0◦ to 360◦ . We used the high frequency asymptotic approach of Frederiksen and Bostock (2000) applied to a simple three-layer model composed of i) a 40 km homogeneous crust (Poisson solid, Vp = 6.6 km/s, ρ = 2900 kg/m3 ), ii) a homogeneous half-space mantle (Poisson solid, Vp = 8.1 km/s, ρ = 3300 kg/m3 ), and iii) a dipping 10 km thick layer exhibiting hexagonal anisotropy with average 42  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada mantle velocities. Hexagonal anisotropy is parametrized by the variation in P - and S -wave velocity (∆Vp and ∆Vs ), and the plunge and trend of the symmetry axis (e.g. Figure 2 of Sherrington et al. (2004)). Although we have specified Poisson solids for layers i) and ii), we note that S -velocity exerts the main control on amplitude of response. The dip and strike of the layer also influence the azimuthal response. Based on the timing of AM, we fixed the dip of the layer to 20◦ , ∆Vp and ∆Vs to ±5% and constrained the strike to the interval of −20◦ to 20◦ . The plunge and trend of the symmetry axis were allowed to vary from 0◦ to 90◦ , and 0◦ to 180◦ , respectively.  43  Figure 3.6: Cartoons showing (a) velocity ellipsoid corresponding to modelled hexagonal anisotropy with fast symmetry axis, (b) definition of trend and plunge of symmetry axis, and (c) orientation with respect to suture model.  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  44  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada An exploration of this reduced model space, employing a simple grid search with a sampling interval of 10◦ in each dimension, allowed us to obtain a model that reproduces the characteristic azimuthal signature of AM reasonably well. A solution was obtained with a model involving a strike of 0◦ , a plunge of 70◦ , and a trend of 40◦ as defined in Figure 3.6. Careful investigation revealed that no model comprising a fast symmetry axis in the plane of the dipping layer could successfully reproduce the observations. Figure 3.7 compares the azimuthal responses for the data and preferred model. It shows a good correspondence between observed and predicted polarities and amplitudes on both the transverse and radial components. Note, in particular, that the synthetic back-azimuthal response reproduces the dominantly 1-θ pattern with polarity cross-overs occurring at back-azimuths of approximately 60◦ and 240◦ .  3.6.2  One-way modeling  To further validate this model and accommodate laterally varying structural dip along the array we employed a modeling code based on the wide-angle, one-way wave equation (Thomson, 1999) as described by Audet et al. (2007). We adapted the slab geometry from Cook et al. (1998) and represented the structural elements as a 15 km thick layer of dipping anisotropic material, a 37 km crust and a half-space mantle with properties derived from the grid search model and identified in Table 1 (see Figure 3.8). The frequency dependent behavior of AM, documented in section 5, is modeled as a stack of 5 thin layers showing a decrease from strong (±5%) to weak anisotropy from top to bottom. The plunge of the anisotropy is rotated as the dip shallows to ensure that the fast axis maintains the same angle with respect to the plane of the dipping layer. We modelled signals TM1 and TM2 from the teleseismic Moho near the suture by a mantle wedge that displays a linear horizontal gradient in velocity from Vp = 7 km/s to 8.1 km/s over a distance of ∼150 km. These velocities are broadly consistent with the lowvelocity mantle wedge reported by Fern´ andez-Viejo and Clowes (2003) and Fern´ andez-Viejo et al. (2005).  45  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  Radial Component  Transverse Component 0.2  Synthetic Response  6 8  0.15  10 Amplitude Relative to P  0.1  Synthetic Time [s]  6 8 10  Observations  6  0.05 0 −0.05 −0.1  8 −0.15  10  0  100  200  300 0 100 Back Azimuth [°]  200  300  −0.2  Figure 3.7: Radial (left) and transverse (right) back-azimuthal responses. Top row shows synthetic back-azimuthal response for incident waves with a fixed slowness of 0.06 s/km. Middle row shows synthetic impulse responses individually stacked in 3.6◦ back-azimuth bins for the set of slownesses and back-azimuths corresponding to the observed distribution. Bottom row shows observed receiver functions from stations UBC03, UBC04 and A12 aligned along AM and stacked in 3.6◦ back-azimuth bins.  Vp (km/s) 8.1  Vs (km/s) 4.5  ρ (kg/m3 ) 3300  Anisotropy (±%) 5  trend (◦ ) 40  plunge (◦ ) 70  strike (◦ ) 0  dip (◦ ) 20  Table 3.1: Model parameters  46  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  Depth [km]  0 25 50 75 100 125 50  100  150  200  250  Horizontal distance [km] Figure 3.8: Subsurface model of the Wopmay Orogen 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. The maximum amplitude of hexagonal anisotropy is 5% at the top of the layer. 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. The model domain encompasses an area of 256 km×110 km in 128×220 grid-points for a horizontal sampling of 2 km and a vertical sampling of 500 m. Although this simplified model does not account for all the features within the observed receiver functions it provides a good geometrical template of the suture suitable for modelling the main arrivals.  47  Synthetic Seismograms Time [s] 10  5  0  10  5  100  150  200  200  Transverse Component  250 100 150 Distance along the profile [km]  Radial Component  250  Amplitude Relative to P  −0.25  −0.2  −0.15  −0.1  −0.05  0  0.05  0.1  0.15  0.2  0.25  Figure 3.9: Synthetic (top) and observed(bottom) receiver functions for stations A10 to A15. Left panels show the radial receiver function; right panels show the transverse receiver function. The same set of slownesses and back-azimuths was used to generate the one-way synthetic data. The traces filtered between 0.2 and 1.5 Hz. Surface sediments reverberations seen on the real radial component are not modeled in the simulation. Note that the synthetics reproduce the main features of seen in the real radial and transverse component.  Real Seismograms Time [s]  0  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  48  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada To obtain results consistent with the receiver functions, we produced synthetic seismograms for a subset of 50 events from the real data-set and ordered them in like manner by geographical location, bandwidth and P wave slowness. The simulation results are presented in Figure 3.9. The synthetic radial component captures the dominant characteristics of the Moho, namely, the variable dip and the disappearance of signature at depth below the suture. The polarity is consistent with observation, being generally positive along the profile, except where the signal dips between back-azimuths of 220◦ and 310◦ . On the synthetic transverse component, the dipping layer is clearly identified and exhibits a similar frequency dependence to AM (see Figure 3.10).  3.7  Interpretation  In this section, we place our observations in the context of several other seismic studies previously undertaken in the study area. We compare our results with the SNORCLE near-vertical reflection (Cook et al., 1998; Cook and Vasudevan, 2003) and wide-angle reflection/refraction studies (Fern´ andez Viejo et al., 1999; Fern´ andez-Viejo and Clowes, 2003), and teleseismic receiver function analysis from the Yellowknife array (Bostock , 1998).  3.7.1  SNORCLE active source studies  To better compare the SNORCLE near-vertical reflection profile and the leg A CANOE data, we superimpose the time migrated, coherency filtered, reflection section (dominant frequency >15 Hz) upon the radial- and transverse-component receiver functions (dominant frequency < 1 Hz) in Figure 3.11. This is accomplished through a simple scaling between the direct P -to-S conversion time, tps , and the 2-way P travel-time, t2p in the form 2tps t2p =  1 − Vp2 p2  3 − Vp2 p2 −  ,  (3.2)  1 − Vp2 p2  49  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  Synthetic Time [s]  4 6 8 10 12  CANOE Time [s]  4 6 8 10 12 50  100  150 200 250 Back Azimuth [°]  300  350  Normalized Amplitude  Figure 3.10: Comparison of synthetic (top) and real (bottom) response for stations UBC03, UBC04 and A12 for the same suite of back-azimuth and slowness. The receiver functions are aligned along AM and sorted by back azimuth. Left panels compare the azimuthal response. Right panels compare the amplitude of the stacked impulse responses.  50  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada where p is the horizontal slowness, and a Poisson solid is again assumed. For Vp = 6 km/s and p = 0.06 s/km, the scaling relation is t2p ≈ 2.6 tps . The signature of the topmost (1 km thick) sedimentary sequence S is manifest in different ways on the two profiles. The base of the sediments is barely visible in the reflection data as a near-continuous reflector at t2p < 1 s between stations A09 and A16. Its expression on the radial component receiver functions extends to significantly later times (t2p ∼ 5 s) because i) free-surface reverberations dominate the response and are mislocated with depth as forward conversions using (3.2); and ii) the teleseismic response exhibits lower frequency bandwidth, resulting in greater pulse widths. On the radial component, we observe a good correspondence between the teleseismic Moho, TM, and the sudden change in reflectivity with depth on the reflection profile, interpreted by Cook et al. (1998) as being the base of the crust. Over much of the profile, this feature remains relatively flat. At the Fort Simpson/Hottah terrane suture, however, TM follows the base of a dipping layer on the reflection profile over a distance of approximately 50 km. This layer is interpreted to be remnant crustal material subducted beneath the Hottah terrane during the Proterozoic. Note that the dipping signal AM visible on the transverse receiver function, does not coincide with the subducted crust, as inferred from the reflection profile, but rather parallels it, sharing a common interface that is the teleseismic (and reflection) Moho. An important implication of this observation is that the anisotropic layering identified on the teleseismic response represents a lid of shallowmost mantle material.  51  40  35  30  25  20  15  10  5  0  40  35  30  25  20  15  10  5  0  0  50  A  TM  09  100  150  AM  TM1  200  250  A1  TM2  05 06 07 08 BC BC 13 BC BC 14 U U A U U A  A2  300  09 10 BC BC 15 U U A  Distance Along the Profile [km]  p. im 01 02 03 04 S 0 . BC BC 11 BC BC 2 Ft A1 U U A U U A1  350  . Ft  .  ov Pr 6  400  TM  A1  450  A1  7  150  110  70  30  0  150  110  70  30  0  Estimated Depth [km] Estimated Depth [km]  Figure 3.11: (a) Radial and (b) polarity-corrected transverse receiver functions upon which time migrated, coherency-filtered, near-vertical reflection data are superimposed with the position of the stations along the profile shown on top. The display format is slightly different from the one adopted in Figure 3.4 in the sense that receiver function are projected along the profile instead of the west-east direction. To agree with two-way P travel time of the reflection profile a simple time scaling, described in (3.2) was applied. An estimate of the depth based on a simplistic two layers featuring a crust and mantle with velocity of 6 and 8 km/s, respectively is provided on the right axis. The teleseismic responses are filtered between 0.2 and 1.5 Hz. The position of the stations along the profile is shown on top.  Transverse Component Two Way P−Travel Time [s]  b)  Radial Component Two Way P−Travel Time [s]  a)  A  08  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  52  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada East of the Fort Simpson/Hottah terrane suture (∼ kilometer 150-200 along the profile), we observe a temporary disappearance of the Moho in the depth range of 30-50 km on both the teleseismic and the near-vertical reflection profiles suggesting a reduced velocity contrast. Wide-angle reflection/refraction data indeed show high crustal and low mantle velocities along this stretch of profile (Fern´ andez-Viejo and Clowes, 2003). Further study by Fern´ andez-Viejo et al. (2005) that combined both P - and S -wave observations yielded additional insights into the nature of the wedge. They demonstrated that the wedge region is characterized by low Vp /Vs or equivalently low Poisson’s ratios, and attributed this observation to the presence of pyroxenite. Slightly to the east (kilometers 300-400 along the profile), a re-analysis of near-vertical reflection data by Cook and Vasudevan (2003) permitted the identification of synclinal mantle reflections that are discordant with the overlying Moho and which they must therefore predate. The authors suggested that these structures may have developed during the 1.8 Ga subduction event, but the lithologies responsible are difficult to constrain. It is interesting to note that the teleseismic results present evidence for structural heterogeneity within the intermediate wedge region (signals A1 and A2 at kilometers 200-250 in Figure 3.11) that appears to possess anticlinal geometry although the true geometry is difficult to ascertain given the simple imaging approach adopted here. The 1-θ polarity correction applied to improve lateral coherence of AM on the transverse component also appears to enhance the continuity of the A1 and A2 structures and may indicate that these features are all dynamically related. East of Fort Providence, both CANOE leg A and SNORCLE reflection line 1 proceed north-east toward Yellowknife. Along this section of the Mackenzie Highway the teleseismic station distribution is sparser and the quality of the data allows only few coherent signals to be identified. Station A16 does appear to show near-horizontal layering at depths near 70 km, in particular on the radial component, that coincides with mantle reflectivity on SNORCLE line 1. As for the reflection profile, however, it is difficult to determine the relation of this structure with the main suture to the west 53  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada near kilometer 150.  3.7.2  Comparison of CANOE and Yellowknife array results  The Yellowknife array at the eastern terminus of leg A affords over 15 years of high quality broadband three-component data well sampled in back-azimuth and epicentral distance. As mentioned in the introduction, Bostock (1998) identified a series of anisotropic lithospheric mantle reflectors beneath Yellowknife: H at 70-80 km, X at 120-150 km and L at 170-230 km depth. He speculated that L could represent the continuation of the Proterozoic subducted plate observed to the west. Moreover, he interpreted the layering, and in particular H, as evidence for stranded oceanic crust that had developed anisotropy through a structural preferred orientation of garnet and clinopyroxene mineralogies during the process of eclogitization. The signature of H documented by Bostock (1998) and that of AM observed here, bear a close similarity. Figure 3.12 compares the two sets of aligned signals as AM at stations UBC03, UBC04 and A12 and H at YKW1, YKW2, YKW3 and YKW4. Both seismogram sections and their summary traces exhibit a sharp earlier arrival followed by a more diffuse later pulse (note, however, that AM displays a lower amplitude, sharp, negative precursory pulse that is not present for H ). Like the signal AM, H also has a counterpart on the near-vertical reflection profile in the form of a single reflection that is evident beneath the western Slave province (the Anton Terrane) including the vicinity of the Yellowknife array, at a depth of 75 km [Cook et al. (1998)]. It appears to be paired with a second reflection 5 km deeper over a short 10 km stretch to the east which is roughly 5 km deeper. The mantle reflection (like H ) appears to be dominantly horizontal beneath the Slave province but shallows to the west beneath the Great Bear magmatic arc. Anisotropy is clearly manifest on H in polarity reversals at ∼180◦ , 270◦ back azimuth on the transverse component and at ∼320◦ on the radial component, implying a dominantly 2-θ symmetry. In comparing this response with that of AM, we corrected for the effect of slab dip in the latter feature. The corrected orientation does produce a polarity reversal near  54  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada  Yellowknife Time [s]  4 6 8 10 12  CANOE Time [s]  4 6 8 10 12 0.045  0.05  0.055 0.06 0.065 Slowness [s/km]  0.07  0.075  Normalized Amplitude  Figure 3.12: Comparison of receiver functions for Yellowknife array (signal H, top) and CANOE (signal AM, bottom). Impulse responses are filtered between 0.2 and 1.5 Hz and plotted between 0 and 15 s and the maximum of each trace is normalized. The left panel shows events in the back-azimuth range of 275◦ and 312◦ , sorted by horizontal slowness and moveout corrected. The right panel shows stacked receiver functions in black for Yellowknife (top) and CANOE (bottom). The gray dashed line shows the stacked receiver function for the other array, i.e for CANOE (top) and Yellowknife (bottom). 240◦ back azimuth on the transverse component. The isotropic velocity contrast is, however, insufficient to produce any polarity reversal on the radial component, and the corrected symmetry axis remains at a high angle to the horizontal resulting in a complex combination of 1-θ, 2-θ, and 3-θ contributions.  3.8  Discussion and Conclusions  The observations documented in previous sections bear important implications for our understanding of the structure and origin of continents. Figures 3.11a,b demonstrate that the high-frequency seismic reflection and low-frequency teleseismic signatures of fossil subduction are distinct and 55  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada manifest different elastic properties. The reflection response illuminates finer-scale elasticity contrasts, whereas the teleseismic signature identifies the presence of an anisotropic mantle lid within the subducted lithosphere. In this section we shall examine candidate mechanisms by which this localized (10-15 km thick) zone of mantle material, AM, may acquire its anisotropy, under the assumption that it represents lattice preferred orientation of olivine (Babuska and Cara, 1991). We proceed then to examine its relation to the structure H observed beneath the Slave province and implications thereof for models of continental evolution. There are two obvious candidate mechanisms for the generation of AM, namely: 1) anisotropy has developed in situ as an immediate consequence of Paleoproterozoic subduction/suture, or, 2) it represents an original fabric developed during genesis of the (later subducted) lithosphere at the parent mid-ocean ridge. It is also conceivable that the signature of AM comprises contributions from both mechanisms; in any event, stability of mantle fabric over close to 2 Ga of Earth history is implied. Studies of upper mantle peridotites from volcanic xenoliths and orogenic massifs reveal that anisotropic fabrics generally develop at high to very high temperature (> 1100◦ C) conditions (Nicolas, 1976). In modern-day subduction zone environments, such temperatures are unlikely to exist within the structural configuration represented by AM, i.e. mantle lid of subducting plate at shallow depths. However, it has been observed that lower temperature deformation of olivine can take place within 2-10 km thick shear domains that have been ascribed to transform faults (Prinzhofer and Nicolas, 1980). It would appear reasonable, therefore, that low-temperature, localized deformation might also occur in subduction environments, especially in the near vicinity of the plate interface. The location of AM immediately below the former Moho and the reduction in strength of anisotropy with depth both suggest the development of simple shear strain created by stresses transmitted from the overlying plate boundary. It is necessary then to identify the factor controlling the scale of layer thickness (10-15 km). One candidate mechanism involves the eclogitization of the overlying oceanic crust. In this scenario, the ∼10% volume change accompanying 56  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada eclogitization within the crust creates a dimensional mismatch with underlying mantle and produces compressive stress (and consequent deformation) over a depth interval below the crust-mantle boundary that is roughly commensurate with our observations (see Kirby et al. (1996), their figure 5b). A difficulty with this interpretation is that it is inconsistent with definition of the Moho as presented on the receiver function profiles in Figure 3.11. Because eclogite and peridotite are characterized by comparable velocities (e.g. Christensen (1996)), we expect the transformation of crustal lithologies to eclogite to be marked by a disappearance of the Moho. The Moho is visible to a depth of ∼50 km beneath the suture (kilometer 175), setting the minimum depth of eclogitization, whereas the anisotropic mantle is well defined on the transverse component to significantly shallower levels corresponding to normal Moho depths (i.e. ∼ 30 km depth at kilometer 100). An alternative possibility is that bending strains at the suture play a role in development of anisotropy and layer-thickness scale, although this mechanism is difficult to evaluate without more detailed modelling. Oceanic sampling and ophiolite studies suggest that dunites and olivinerich harzburgites predominate within the oceanic lithosphere from the Moho to depths of ∼ 25km (Nicolas et al., 1980). These rocks are plastically deformed at the high temperatures prevailing near the ridge, with particularly large strains developing close to the Moho. Orientations of olivine and pyroxene crystals in the Bay of Islands ophiolite appear to show a systematic relation wherein orthopyroxene c crystallographic axes roughly correspond to the a axes of olivine grains (Christensen and Lundquist, 1982). The net effect is a reduction in elastic anisotropy as the percentage of orthopyroxene within peridotite increases. These observations present an alternative interpretation for AM, that is, that it represents a gradient zone of increasing orthopyroxene content (with corresponding reduction in anisotropy) and decreasing deformation between the Moho and 10-15 km depth developed at a Paleoproterozoic ridge. A combination of compositional and strain gradients may therefore define the thickness of AM. The steeply dipping plunge of the symmetry axis documented above may imply a slow spreading rate whereby asthenospheric flow at the ridge axis is frozen in with an 57  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada intrusive, dominantly sub-vertical attitude, since the flow geometry beneath faster spreading centers is generally subhorizontal (Nicolas and Violette, 1982; Boudier and Nicolas, 1985). A potential problem with both subduction and mid-ocean ridge mechanisms for the origin of AM concerns the fact that, to the best of our knowledge, there is no comparable documentation of depth-localized anisotropy within the mantle lid of subducting slabs at shallow depths in the modern context, despite a growing body of research. Notwithstanding the absence of a conclusive mechanism of formation for AM, its identification with paleosubduction is important in a range of contexts. Let us first consider the Slave Province and the remarkable resemblance between the structure AM and H. If AM and H represent the same subduction episode and form part of one contiguous structure, the implication is that the bulk of mantle lithosphere beneath the western Slave province is Paleoproterozoic (ca. 1.8 Ga) in age. However, the interpretation of these two structures as one and the same, although conceivable, would likely require shallowing toward the eastern end of the profile and stands at odds with recent evolutionary models for the Slave province based on other geophysical, geological and geochemical observations that favor development of a root by 2.6 Ga (Davis et al., 2003). We prefer therefore to interpret the similarity as due to a common, generic signature of fossil subduction as recorded in two distinct episodes of plate convergence, with H presumably emplaced in the late Archean and AM in the paleo-Proterozoic. It is unfortunate that neither the seismic reflection nor teleseismic profiles are able to definitively establish the geometric (continuous or discontinuous) relation between AM and H, X or L due to the highly variable and intermittent characteristics of the signals between the two locations. We should note that the original interpretation of Bostock (1998) that the anisotropic layering beneath YKA is the signature of subducted lithosphere is now corroborated in general, although it errs in detail because the anisotropic layering responsible for AM (and by inference H ) can now clearly be identified with mantle material as opposed to eclogitized, former oceanic crust (although eclogitized crustal material presumably immediately overlies H as it does AM ). This observation also renders moot another dif58  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada ficulty with the former interpretation, namely, that eclogite xenoliths from the Slave Province fail, for the most part, to produce a significant S-wave anisotropy (to which forward scattered conversions are most sensitive) as measured in the laboratory (Kopylova et al., 2004). Confirmation that ancient subducted lithosphere is characterized by finescale (∼10 km thick) anisotropic layering lends further credibility to the thesis that Archean cratons were stabilized through shallow subduction processes (e.g., Vlaar (1986); Helmstaedt and Schulze (1989); Abbott (1991)). Another few hundreds of km to the northeast near the center of the Slave province where the major diamond production areas are located, Snyder (2007) has identified similar anisotropic layering at depths near 120 km exhibiting a south-easterly dip. It would appear likely that this layering is related to one of the structures H, X, or L beneath YKA, and that it therefore signifies cratonic-scale lithospheric underthrusting. There are at least 3 published examples of comparable structures defined by teleseismic receiver functions on other cratons. On the Indian shield, Saul et al. (2000) observe a direct P-to-S conversion from the top of a similar layered structure at 90 km depth below the Dharwar craton (station Hyderabad) on both radial and transverse components. Although, the back-azimuthal distribution of events for stations in central India is limited, the strength of the signal and its moveout argue strongly for an origin in near horizontal, depth localized anisotropy. Levin and Park (2000) have modelled forward conversions from an anisotropic layer structure at 70-90 km depth beneath station RAYN on the Arabian Shield which, like AM and H, is characterized by a lower boundary that exhibits a discontinuity in material property gradient, as opposed to discontinuous material properties. Their interpretation associates the anisotropic layer with a shear zone developed in peridotitic mantle material, although they attributed its development to continental collision rather than to the underthrusting of a subducting plate. A study of discontinuity structure beneath the Kaapvaal craton using stacked data from the portable South African experiment (Gao et al., 2002) failed to find evidence for significant lithospheric layering beneath the Moho. Their analysis was based, however, on generally coarse 59  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada spatial sampling and the assumption that lithospheric mantle structures are dominantly isotropic and horizontal such that possible back azimuthal variations are ignored. A separate examination of data from the subsidiary Kimberly array (Snyder et al., 2003) does strongly suggest the presence of thin, anisotropic layering exhibiting variations over short distances consistent with the geographical patterns of anisotropy expressed through SKS splitting and PKP delay times (Fouch et al., 2004). We would suggest that similar structures may be found beneath other portions of Kaapvaal craton but that lateral variability in both structural geometry and anisotropic fabric preclude their detection with coarse station sampling and limited duration deployments. In summary, a growing body of evidence is pointing to the importance of shallow subduction in the stabilization of early continents. The results presented here from the Wopmay orogen in Canada’s Northwest Territories provide unambiguous evidence for the signature of fossil subduction in the form of highly anisotropic, shallow-most mantle characterizing the underthrust plate. This anisotropic layering is comparable in dimension and teleseismic response to structures observed both beneath the Slave province and cratons worldwide, and which we, accordingly, infer to share a common origin. The identification of this fossil signature may provide an important diagnostic for future studies of deep continental structure and evolution.  3.9  Acknowledgments  We gratefully acknowledge the contribution of instruments from IRIS-PASSCAL, POLARIS, Arizona State University, Georgia Institute of Technology, and the University of British Columbia. Andy Langlois provided invaluable assistance in coordination of the field effort. Discussions with Jounada Oueity and Ron Clowes helped to focus our initial analysis. We are grateful to Associate Editor Frederik Simons, Sebastian Rost and an anonymous reviewer for comments that helped improve manuscript presentation, and to Nikolas Christensen for discussion on the nature of localized mantle anisotropy. This work was financially supported by the Natural Sciences and Engineering 60  Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada Research Council of Canada (NSERC-Lithoprobe supporting geoscience grant CSP0006963 to MB) and National Science Foundation (grants NSF EAR-0453747 to JG, NSF EAR-0711401 to EG and NSF EAR-0003745-004 to JR).  61  Bibliography Abbott, D., The case for accretion of the tectosphere by buoyant subduction, Geophys. Res. Lett., 18, 585–588, 1991. Audet, P., M. G. Bostock, and J.-P. Mercier, Teleseismic waveform modelling with a one-way wave equation, Geophys. J. Int., doi:10.1111/j. 1365-246X.2007.03586.x, 2007. Babuska, V., and M. Cara, Seismic Anisotropy in the Earth, Kluwer, Dordrecht, 1991. Bostock, M. G., Mantle stratigraphy and evolution of the Slave province, J. Geophys. Res., 103, 183–200, 1998. Bostock, M. G., Seismic waves converted from velocity gradient anomalies in the Earth’s upper mantle, Geophys. J. Int., 138, 747–756, 1999. Boudier, F., and A. Nicolas, The harzburgite and the lherzolite subtypes in ophiolitic and oceanic environments, Earth Planet. Sci. Lett., 76, 84–92, 1985. Bowring, S. A., and I. S. Williams, Priscoan (4.00-4.03 ga) orthogneisses from northwestern Canada, Contrib. Mineral Petr., 134, 3–16, 1999. Burridge, R., M. V. D. Hoop, D. Miller, and C. Spencer, Multiparameter inversion in anisotropic elastic media, Geophys. J. Int., 134, 757–777, 1998. Christensen, N. I., Poisson’s ratio and crustal seismology, J. Geophys. Res., 101, 3139–3156, 1996. Christensen, N. I., and S. M. Lundquist, Pyroxene orientation within the upper mantle, Geol. Soc. Am. Bull., 93, 279–288, 1982. Cook, F. A., and P. Erdmer, Introduction to special issue on the lithoprobe Slave-northern cordillera lithospheric evolution (SNORCLE) transect., Can. J. Earth Sci., 42, 865–867, doi:10.1139/E05-067, 2005.  62  Bibliography Cook, F. A., and K. Vasudevan, Are there relict crustal fragments beneath the moho?, Tectonics, 22, doi:10.1029/2001TC001341, 2003. Cook, F. A., A. J. Van der Velden, K. W. Hall, and B. J. Roberts, Tectonic delamination and subcrustal imbrication of the precambrian lithosphere in northwestern canada, mapped by lithoprobe, Geology, 26, 839–842, 1998. Cook, F. A., A. J. Van der Velden, K. W. Hall, and B. J. Roberts, Frozen subduction in Canada’s Northwest Territories: Lithoprobe deep lithospheric reflection profiling of the western canadian shield, Tectonics, 18, 1–24, 1999. Davis, W. J., A. G. Jones, W. Bleeker, and H. Grutter, Lithosphere development in the Slave craton: a linked crustal and mantle perspective, Lithos, 71, 575–589, 2003. Fern´ andez-Viejo, G., and R. M. Clowes, Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay Orogen, nortwestern Canada, from a lithoprobe refraction/wide-angle reflection survey, Geophys. J. Int., 153, 1–19, 2003. Fern´ andez Viejo, G., R. Clowes, and J. Amor, Imaging the lithospheric mantle in northwestern Canada with seismic wide-angle reflections, Geophys. Res. Lett., 26, 2809–2812, 1999. Fern´ andez-Viejo, G., R. M. Clowes, and J. K. Welford, Constraints on the composition of the crust and uppermost mantle in northwestern Canada: Vp/Vs variations along lithoprobe’s SNorCLE transect, Can. J. Earth Sci., 42, 1205–1222, 2005. Fouch, M. J., P. G. Silver, D. R. Bell, and J. N. Lee, Small-scale variations in seismic anisotropy near Kimberley, South Africa, Geophys. J. Int., pp. 764–774,doi: 10.1111/j.1365–246X.2004.02,234.x, 2004. Frederiksen, A. W., and M. G. Bostock, Modelling teleseismic waves in dipping anisotropic structures, Geophys. J. Int., 141, 401–412, 2000. Gao, S. S., P. G. Silver, and K. H. Liu, Mantle discontinuities beneath southern Africa, Geophys. Res. Lett., 29, 1491, doi: 10.1029/2001GL013834, 2002. Helmstaedt, H., and D. J. Schulze, Southern African kimberlites and their mantle sample: Implications for Archean tectonics and lithosphere evolu63  Bibliography tion, in it Kimberlites and Related Rocks, 358-368 pp., Blackwell, Cambridge, Mass., 1989. Hildebrand, R. S., P. F. Hoffman, and S. A. Bowring, Tectono-magmatic evolution of the 1.9 Ga Great Bear Magmatic zone, Wopmay Orogen, NW Canada, J. Volc. Geotherm. Res., 32, 99–118, 1987. Hoffman, P. F., and S. A. Bowring, Short lived 1.9 Ga continental margin and its destruction, Wopmay orogen, Geology, 12, 68–72, 1984. Isachsen, C. E., and S. A. Bowring, Evolution of the Slave craton, Geology, 22, 917–920, 1994. Kennett, B. L. N., The removal of free surface interaction from threecomponent seismograms, Geophys. J. Int., 104, 153–163, 1991. Kirby, S., E. R. Engdahl, and R. Denlinger, Intermediate-depth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs, Geophysical monograph, 96, 195–214, 1996. Kopylova, M. G., J. Lo, and N. I. Christensen, Petrological constraints on seismic properties of the Slave mantle (Northern Canada), Tectonophysics, 77, 493–510, 2004. Langston, C. A., Structure under mount Rainier, Washington, inferred from teleseismic body waves, J. Geophys. Res., 84, 4749–4762, 1979. Levin, V., and J. Park, Shear zones in the Proterozoic lithosphere of the Arabian Shield and the nature of the Hales discontinuity, Tectonophysics, 323, 131–148, 2000. Nicolas, A., Flow in upper-mantle rocks: some geophysical and geodynamic consequences, Tectonophysics, 32, 93–106, 1976. Nicolas, A., and J. F. Violette, Mantle flow at oceanic spreading centers: models derived from ophiolites, Tectonophysics, 81, 319–339, 1982. Nicolas, A., F. Boudier, and J. L. Bouchez, Interpretation of peridotite structures from ophiolitic and oceanic environments, Am. J. Sci., 280, 192– 210, 1980. Prinzhofer, A., and A. Nicolas, The Bogota Peninsula, New Caledonia: a possible oceanic transform fault, J. Geology, 88, 387–398, 1980. 64  Bibliography Saul, J., M. R. Kumar, and D. Sarkar, Lithospheric and upper mantle structure of the Indian Shield, from teleseismic receiver functions, Geophys. Res. Lett., 27, 2357–2360, 2000. Sherrington, H. F., G. Zandt, and A. W. Frederiksen, Crustal fabric in the Tibetan Plateau based on waveform inversion for seismic anisotropy parameters, J. Geophys. Res., 109, B02,312, 2004. Snyder, D. B., Stacked uppermost mantle layers within the Slave craton of NW Canada as defined by anisotropic seismic discontinuities, Submitted to Tectonics, 2007. Snyder, D. B., S. Rondenay, M. G. Bostock, and G. D. Lockhart, Mapping the mantle lithosphere for diamond potential using teleseismic methods, Lithos, 77, 859–872, 2003. Thomson, C. J., The ’gap’ between seismic ray theory and ’full’ wavefield extrapolation, Geophys. J. Int., 137, 364–380, 1999. Villeneuve, M. E., R. J. Theriault, and G. M. Ross, U-Pb ages and Sm-Nd signature of two subsurface granites from the Fort Simpson magnetic high, NW Canada, Can. J. Earth Sci., 28, 1003–1008, 1991. Vinnik, L. P., Detection of waves converted from P to SV in the mantle, Phys. Earth Planet. In, 67, 39–45, 1977. Vlaar, N. J., Archean global dynamics, Geology, 65, 91–101, 1986. Wiener, N., Extrapolation, interpolation and smoothing of stationary time series, New York: Wiley, 1949.  65  Chapter 4  Body-wave tomography of western Canada  4.1  Introduction  Western Canada is composed of a wide variety of geologic terranes that range from stable craton to active orogen and which encompass the last 4.0 Ga of the planet’s history. Accordingly, western Canada offers an ideal setting to study subsurface variations in upper-mantle structure, geometry, and physical properties from Archean to the present day. The study of the western margin also provides an opportunity to explore the complex transition from convergent to transform boundary. In this paper we undertake a continental-scale study of the western Canada upper mantle using teleseismic body-wave tomography. Studies of western Canada upper-mantle velocity structure began more than 3 decades ago when Buchbinder and Poupinet (1977), Wickens (1977), and Wickens and Buchbinder (1980) used short-period P-waves, surfacewave, and long-period S-wave traveltime residuals, respectively, to document a low-velocity region beneath British Columbia. More recently, global and continental-scale studies that include western Canada, using both bodywaves (e.g. Grand (1994); Grand et al. (1997)), and surface-waves (e.g. A version of this chapter has been submitted to the Journal Tectonophysics in September 2008. Mercier, J.-P., M. G. Bostock, J. F. Cassidy, K. Dueker, J. B. Gaherty, E. J. Garnero, J. Revenaugh, and G. Zandt (2008), Body-wave tomography of western Canada , submitted to Tectonophysics c 2008 Elsevier.  66  Chapter 4. Body-wave tomography of western Canada Frederiksen et al. (2001); Van Der Lee and Frederiksen (2005)) have imaged a relatively sharp mantle transition from low to high velocity across the Cordilleran deformation front. At finer scales Bostock and VanDecar (1995) utilized body-wave traveltimes to constrain the geometry of the subducting Juan de Fuca plate below southwestern British Columbia whereas Frederiksen et al. (1998) and Shragge et al. (2002) focused on details of upper mantle velocity structure beneath the southern Yukon and south-central Alberta, respectively. These previous studies have provided useful insights into western Canadian lithospheric and upper-mantle structure but are either global in coverage and suffer from low resolution, or local and lack a broader perspective. Until recently, the station density in most of western Canada precluded the undertaking of a comprehensive body-wave study over the entire region as performed for example, in the western United States (Humphreys and Dueker , 1994). However, the recent deployment of stations from several portable experiments (i.e. CANOE, BATHOLITHS, POLARIS BC and POLARIS-Nechako) have contributed to filling some major gaps in station coverage and afford opportunity for a broader scale study. In this paper, we exploit both newly available and previously analysed broadband and short-period data to examine the large-scale, upper-mantle P - and S velocity structure beneath western Canada and hence provide constraints on the evolution and physical properties of this region. The results presented here also complement the most recent high resolution velocity models of the western U.S. calculated using the newly available data from the USArray (Burdick et al., 2008; Sigloch et al., 2008) by extending coverage northward.  4.2 4.2.1  Data and Method Data  The traveltime data employed in this study to construct P - and S -wave velocity models were recorded at stations of the following permanent and temporary networks (see Figure 4.1): Advanced National Seismic System  67  Chapter 4. Body-wave tomography of western Canada  65 oN  60 oN  55 oN  50 oN BATHOLITHS CANOE CNSN POLARIS  45 oN  OTHERS  144 o W 136 oW  128oW  o 120 W  o 112 W  Figure 4.1: Map of western Canada illustrating the distribution of braodband and short-period seismic stations used in this study  (ANSS), Alaska Tsunami Warning System (ATWS), Alaska Regional Network (ARN), BATHOLITHS, CAnadian NOrthwest Experiment (CANOE), Canadian National Seismograph Network (CNSN), and several experiments of the Portable Observatories for Lithospheric Analysis and Research Investigating Seismicity (POLARIS). The P -wave data-set consists of 23, 420 teleseismic P delay times from 1609 earthquakes recorded at 234 broadband and short-period stations at epicentral distances from 30◦ to 100◦ . The S -wave data-set includes 15, 805 teleseismic S delay times from 884 events recorded at 194 broadband stations. Azimuthal coverage is good for both data-sets with only one unsampled sector between 200◦ and 250◦ (Figure 4.2). To facilitate picking, P and S data were divided into independent subsets based on station location. For each subset, P and S traveltimes were obtained from visual picks of the vertical and transverse component, respectively, which were subsequently refined through multichannel crosscorrelation (VanDecar and Crosson, 1990). Frequency bands of 0.4 Hz to 1.6 Hz and 0.05 Hz to 0.4 Hz were used for P and S phases, respectively.  68  Chapter 4. Body-wave tomography of western Canada  30° 60° 90° 120°  Figure 4.2: Equidistant azimuthal projection centered at 55◦ N and 118◦ W illustrating the distribution of events with sufficiently high signal-to-noise ratio recorded at one or more stations of the array. Note the relatively uniform azimuthal coverage with only one region between 200◦ and 250◦ which is unsampled  69  Chapter 4. Body-wave tomography of western Canada The timing uncertainty is estimated using the standard deviation of the residual associated with each trace and is on average ∼30 ms for P and ∼120 ms for S.  4.2.2  Model parametrization  The region of interest in this study extends over more than 23◦ of latitude from 43.5◦ N to 67.27◦ N, over 41◦ of longitude from 105◦ W to 146◦ W and from the surface to 700 km depth. For both the P- and the S-wave inversions the velocity models are parametrized in splines under tension constrained by a series of regularly spaced knots (Figure 4.3). The grid is composed of 175,770 knots, 81 in longitude, 70 in latitude and 31 in depth. The dimensions of the smallest elements of the model are 0.33◦ longitude by 0.5◦ latitude by 20 km depth, whereas the dimension of the largest elements are 1◦ longitude by 1◦ latitude by 25 km. The absolute dimension of the element spacing varies considerably due to the large latitudinal range. The 1D radial earth model IASP91 (Kennett and Engdahl , 1991) was chosen as a reference.  4.2.3  Traveltime inversion  The traveltime inversion procedure adopted in this study to recover slowness perturbations at every point of our model grid is discussed in detail in VanDecar (1991). This technique produces a solution that represents the minimum variation in seismic velocities required to fit the data by imposing conditions on the first and second spatial derivatives of the model which control the flatness and smoothness, respectively. The inversion solves the following problem utilizing a combination of conjugate gradient iterations and robust Huber down-weighting (see Huber (1981)): min Φd + λΦm  (4.1)  x  where λ is a free parameter that controls the level of smoothness and flatness def  of the model, Φd =  Ax − b  2,  def  Φm =  ∆˙s  2  +µ  ∆¨s  2.  Here, ∆˙s  and ∆¨s are vectors containing the first and second derivatives of the slow-  70  Chapter 4. Body-wave tomography of western Canada  Figure 4.3: Three-dimensional view of the grid that constrains the velocity perturbation model with the coast line on top. The velocity model extends from 43.5◦ N to 67.27◦ N and 105◦ W to 146◦ W, and from the surface to 700 km depth. Knots are positioned at the intersection of the thin black lines. Thick black lines indicate the region to which velocity maps and cross sections are confined to. Inverse triangles indicate the position of seismic stations  71  Chapter 4. Body-wave tomography of western Canada ness perturbations vector at each model point, respectively. The remaining quantities are, A=  P C E   ∆s  ,  (4.2)      x= c  e  (4.3)  where the matrix P contains the Frechet derivatives which, at knot j and event i, is given by Pij = δti /δsj . Matrices C and E identify stations and events pertaining to a particular traveltime perturbation, respectively; ∆s, c, and e are vectors containing the slowness perturbation of every point of the model, the station correction term, and the event correction term respectively; and b is the vector of traveltime residuals. Finally, µ is a free parameter that controls the ratio of smoothness to flatness of the model. It should be noted that the procedure adopted in this study provides good estimates of lateral velocity gradients across a region but does not provide a reliable measure of absolute seismic velocity at a given depth (VanDecar , 1991).  4.2.4  Choice of inversion parameters  The choice of the regularization parameter λ is a subjective compromise between fitting the data and obtaining a model that is smooth and flat. In this study, we used λ = 350 and µ = 60 for both P and S inversions. These values were chosen, based on visual inspection of a variety of models, as they represent conservative estimates of the velocity structure favoring smooth and flat models at the expense of fitting of the data. The number of iterations was determined to ensure convergence to a stable solution. A thousand conjugate gradient iterations within 15 Huber down-weighting iterations were used to produce the P model, whereas, due to lower signalto-noise values of the S data-set, 1000 conjugate gradient iterations within 40 Huber down-weighting iterations were required for the S inversion.  72  Chapter 4. Body-wave tomography of western Canada  4.3 4.3.1  Results Resolution Test  The ability of the inversion to recover the true P - and S- velocity models depends both on the ray distribution and the ray geometry (azimuth and incidence angle). We assess the resolution afforded by our P and S data sets through a standard checkerboard resolution test composed of a succession of positive and negative ellipsoidal anomalies (see Figure 4.5). In the test model utilized here, the anomalies are separated by 2◦ in latitude, 3.5◦ in longitude and 200 km in depth. The peak velocity anomaly was set to 2% of the background velocity for both P and S phases. Synthetic traveltimes were calculated for these models using the eventstation ray geometry of the actual data. Gaussian noise with standard deviations of 30 ms and 200 ms was subsequently added to the P and S datasets, respectively. The inversion of the synthetic data-sets was performed using the same procedure and trade-off parameters as employed for the inversion of the real data. Figures 4.6a,b present the synthetic resolution test results for the P and S data sets, respectively. The quality of the recovery is variable and depends, as expected, on the local station density. The resolution is best along the CANOE and BATHOLITHS arrays, and in southwestern British Columbia including Vancouver Island, for depths from 100 to 600 km where station coverage is good and ray density is highest for both P and S models. Elsewhere, the resolution varies from fair to poor depending on the ray density and station geometry. Overall, where the ray density is greater than 0.1, the correlation coefficients between the original and the recovered models are 0.8 and 0.7 for the P and S data sets, respectively. Note that the sub-vertical geometry of the teleseismic body-wave raypaths introduces a vertical smearing that reduces resolution in the vertical direction relative to that in the horizontal plane. The recovered velocity anomalies are systematically underestimated and their magnitudes are typically reduced by about 30% compared to the input anomalies. The underestimated amplitudes are a predictable consequence of introducing the station and event 73  Chapter 4. Body-wave tomography of western Canada  Figure 4.4: Ray density map of the region of interest at 200 km depth for the P data-set. Black represents the highest density whereas white represents the lowest. Note the logarithmic scale that extends from 0 to 1000 rays. The density is the highest along the three legs of CANOE array and in soutwestern British Columbia  74  Chapter 4. Body-wave tomography of western Canada  Figure 4.5: Input model for the checkerboard or eggbox resolution test. This model presents alternating positive and negative ellipsoidale shape velocity perturbations. The blue and red isosurfaces represent a velocity perturbation to the IASP91 reference model of +2% and -2%, respectively.  correction terms and of the regularization which favors smooth, minimum energy models. Note that Figures A.1 to A.10 in the appendix present the input and output of the resolution test for results presented in Figures 4.8 to 4.14.  4.3.2  Raw average station residuals  Raw residuals provide a good preliminary indication of the velocity structure beneath the study region. Positive residuals indicate the presence of lowvelocity structures whereas negative residuals manifest high-velocity anomalies. Figures 4.7a, b combine the average station residuals for the four azimuthal quadrants (northeast, southeast, southwest and northwest) for P and S phases, respectively. We observe a general agreement between the P and S values with only few exceptions along the northern segment of the 75  Chapter 4. Body-wave tomography of western Canada  Figure 4.6: Result of the resolution test given the input model in 4.5 for 76 the P (top) and S (bottom) data-sets. The blue and red isosurfaces represent velocity perturbations of the IASP91 reference model of +1% and -1%, respectively.  Chapter 4. Body-wave tomography of western Canada north-south oriented portion of the CANOE array where residual polarities are opposite. Residuals vary considerably within western Canada suggesting fine-scale lateral heterogeneity. Positive residuals are generally observed west of the Cordilleran deformation front whereas negative residuals are observed east of that boundary. Southwestern British Columbia, where most of the mainland stations east of Vancouver Island display negative residuals, remains a notable exception. The few outliers that do not obey these trends are found, however, throughout western Canada, implying small scale complexity in mantle velocity structure.  4.3.3  P and S Upper Mantle Velocity structure beneath western Canada  Figures 4.8 and 4.9 show 4 horizontal slices through the final P- and Svelocity models, at depths of 100, 200, 300, and 400 km. Note the good correlation between the P and S images; the similarity between the two models, computed from independent data sets, suggests that the inversions are robust.  77  Chapter 4. Body-wave tomography of western Canada  65 oN  60 oN  55 oN  50 oN  +1 sec 45 oN  −1 sec  144 o W 136 oW  128oW  o 120 W  128oW  o 120 W  o 112 W  65 oN  60 oN  55 oN  50 oN  +2.5 sec 45 oN  −2.5 sec  144 o W 136 oW  o 112 W  Figure 4.7: P (top) and S wave (bottom) relative station residuals. For each individual station, four residuals are plotted for the north-east, south-east, south-west, and north-west quadrants. Blue and red symbols corresponds to early and late arrivals, respectively, with symbol size proportional to the delay. 78  −2  −1  % of P velocity 136 oW  tin au lt 128 W  aF o  nt  an ro ler n f rdil atio Co form de  120 W  o  o 112 W  144 o W  45 oN  50 oN  N  65 o  −2  −2  −1  0  1  2  60 oN  55 oN  144 o W 136 oW  % of P velocity 136 oW  400 km  128oW  tin au  aF lt 128oW  o 120 W  nt o 120 W  an ro ler n f rdil atio Co form de o 112 W  o 112 W  Figure 4.8: Depth slices covering the whole survey area at 100 km (top left), 200 km (top right), 300 km (bottom left), and 400 km (bottom right) through the P velocity model. This model was obtained from 1000 conjugate gradient iterations within 15 Huber downweightings using λ = 350 and µ = 60. Blue and red colors represent positive and negative velocity perturbations, respectively.  144 o W  45 oN  0  1  Tin  2  300 km  o 120 W  nt an ro ler n f rdil atio Co form de  60 oN  lt au 128oW  o 112 W  −1  0  % of P velocity  Tin  50 oN  N  % of P velocity 136 oW  45 oN  1  2  lt au  65 o  −2  −1  50 oN  200 km  nt an ro ler n f rdil atio Co form de  55 oN  144 o W  45 oN  0  aF tin Tin  1  2  55 oN  60 oN  65 oN  aF tin Tin  50 oN  55 oN  60 oN  65 oN  100 km  Chapter 4. Body-wave tomography of western Canada  79  −2  −1  % of S velocity 136 oW  tin au lt 128 W  aF o  nt  an ro ler n f rdil atio Co form de  120 W  o  o 112 W  144 o W  45 oN  50 oN  N  65 o  −2  −2  −1  0  1  2  60 oN  55 oN  144 o W 136 oW  % of S velocity 136 oW  400 km  128oW  tin au  aF lt 128oW  o 120 W  nt o 120 W  an ro ler n f rdil atio Co form de o 112 W  o 112 W  Figure 4.9: Depth slices covering the whole survey area at 100 km (top left), 200 km (top right), 300 km (bottom left), and 400 km (bottom right) through the S velocity model. This model was obtained from 1000 conjugate gradient iteration withins 40 Huber downweightings using a λ = 350 and µ = 60. Blue and red colors represent positive and negative velocity perturbations, respectively.  144 o W  45 oN  0  1  Tin  2  300 km  o 120 W  nt an ro ler n f rdil atio Co form de  60 oN  lt au 128oW  o 112 W  −1  0  % of S velocity  Tin  50 oN  N  % of S velocity 136 oW  45 oN  1  2  lt au  65 o  −2  −1  50 oN  200 km  nt an ro ler n f rdil atio Co form de  55 oN  144 o W  45 oN  0  aF tin Tin  1  2  55 oN  60 oN  65 oN  aF tin Tin  50 oN  55 oN  60 oN  65 oN  100 km  Chapter 4. Body-wave tomography of western Canada  80  Chapter 4. Body-wave tomography of western Canada Two important observations can be taken from the P and S models. The first is that most of the low-velocity anomalies are observed in the western portion of the study area whereas the east is mostly fast. The transition from low to high velocity is sharp and well defined at shallow depths, 100 and 200 km, in the south and occurs roughly at the Cordilleran deformation front. North of 55◦ N the transition becomes less clear due to the presence of a high-velocity anomaly west of the Tintina fault but, nonetheless it appears to correspond to the deformation front. In this area, note that the boundary between low and high-velocity regions tends to dip steeply eastward. The second observation concerns a high-velocity anomaly in southwestern British Columbia, east of Vancouver Island. This feature, also eastward dipping, corresponds to the previously documented thermal signature of the Juan de Fuca subducting plate (e.g. VanDecar (1991); Bostock and VanDecar (1995)). We observe an abrupt termination of the slab signature north of Vancouver Island, near 52◦ N. The next section will focus specifically on these two regions.  4.4 4.4.1  Discussion Cordillera/craton transition  Precambrian shields are typically underlain by lithospheric roots that are thicker, colder and more depleted in Fe-rich constituents than the young mantle beneath more recent geological formations (Jordan, 1978). The transition from Phanerozoic mantle to cratonic lithosphere is often characterized by a marked change in physical properties that is typically manifest in tomographic models by a sharp increase in seismic velocities. In recent years, several regional and global tomographic studies have provided images of such transitions in various locations around the globe, where young orogens are juxtaposed with cratons (e.g Grand (1994); Zielhuis and Nolet (1994); Polet and Anderson (1995); Ekstr¨ om et al. (1997); Grand et al. (1997); Simons et al. (1999); Fouch et al. (2004)). In western Canada, Frederiksen et al. (2001) and Van Der Lee and Frederiksen (2005) have documented the  81  Chapter 4. Body-wave tomography of western Canada significant increase in seismic velocity that roughly follows the Cordilleran deformation in the south and Tintina fault in the north. There is general consensus that the shift in seismic velocities documented by passive source imaging (Frederiksen et al., 2001; Van Der Lee and Frederiksen, 2005) defines the transition from Phanerozoic mantle to cratonic lithosphere in the south. In the northern Cordillera, however, the situation is more complex. Despite the fact that Cook et al. (2004) show good evidence for prolongation of a Proterozoic crustal wedge west of the Tintina fault, the extension of the western edge of ancestral North American mantle westward well into the northern Cordillera, as suggested by Frederiksen et al. (1998), Frederiksen et al. (2001), and Van Der Lee and Frederiksen (2005), stands at odds with other geophysical observations. Magnetotelluric (MT) (Jones, 1999; Ledo et al., 2002; Jones et al., 2005) and heat flow data (Hyndman and Lewis, 1999; Lewis et al., 2003), as well as elastic thickness studies (Fl¨ uck et al., 2003; Audet et al., 2007), suggest instead that the transition occurs roughly at the Cordilleran deformation front. The inferred presence of strong thick and cold cratonic roots beneath the northern Cordillera west of the Tintina fault also poses problems with regard to some basic isostatic considerations (Lachenbruch and Morgan, 1990) and in explaining the origin of the surface deformation, which requires a lithosphere that is hot and weak (Hyndman et al., 2005). The stations deployed along the northern portion of the CANOE array afford an excellent lateral resolution of structure across the northern Cordillera. They are consequently well positioned to locate and characterize the transition from allocthonous to autocthonous mantle in northwestern Canada, in the region where the western extent of the craton is not clearly identified. The northern portion of the CANOE array comprises roughly 40 three-component broadband seismic stations separated at variable intervals between 15 km to 60 km. Some 30 of these stations are deployed in the east/west direction, perpendicular to the dominant strike of the major sutures whereas the remaining dozen stations are deployed in the north/south direction sub-parallel to the Cordilleran deformation front. As mentioned in the previous section, the depth slices in Figures 4.8 and 82  Chapter 4. Body-wave tomography of western Canada  Tintina Fault  0  (C)  100  300 400  64 oN 62 oN (A) 60 o N A1 58 oN 56 oN 54 oN 144 o W 13 o 6 W  A2  o W o 128oW 120 W 112  3 500  1.5  600 700 A1 (60.5N, 138.7W)  0 A2 (60.5N, 108.7W)  −1.5  % of P velocity  Depth [km  ]  200  (B)  Cordilleran front  −3  Figure 4.10: Cross section within the P velocity model along the corridor A1A2 (see inset map at top right) in the northern Cordillera from the surface to a depth of 600 km. The topography plotted on top of the cross section is exaggerated 50 times and smoothed. Blue and red represent positive and negative velocity perturbations, respectively.  4.9 indicate that the velocity structure in the northern Cordillera is more complex than in the south. We do image a transition from low to high velocity west of the Tintina fault, but our results suggest that the major transition in seismic velocity occurs farther east at the Cordilleran deformation front. To better characterize this transition we provide in Figure 4.10 a vertical cross section along 60.5◦ N within the P -velocity model, centered on the deformation front, that extends over 15◦ of longitude. On this figure, three important features are labelled A, B and C. The structure A is a large high-velocity anomaly that extends to depths near 300 km beneath the Wopmay orogen and deeper still beneath the Slave craton. It is characterized by a strong velocity anomaly that is consistently 2% to 4% faster than the background value. It is juxtaposed with a low-velocity anomaly of ∼-2% labelled B that extends from beneath A to  83  Chapter 4. Body-wave tomography of western Canada the surface. The transition from A to B is sharp, occurs within less than 50 km and is roughly centered on the Cordilleran deformation front (see Figure 4.11). This transition coincides with an abrupt variation of surface topography from high to low and with a major change in mantle temperature documented by Hyndman and Lewis (1999) and Lewis et al. (2003) from heat flow measurements, and to a lesser extent by Fl¨ uck et al. (2003) and Audet et al. (2007) from the study of elastic thickness. This combination of observations strongly suggests that the transition from B to A corresponds to the boundary between Phanerozoic mantle and Proterozoic lithosphere, thereby shifting the western extent of the craton by hundreds of kilometers eastward from what was proposed by Frederiksen et al. (1998), and Frederiksen et al. (2001) and Van Der Lee and Frederiksen (2005) using body-wave or surface-wave tomography, respectively. The discrepancy in location of the Phanerozoic/craton boundary in northwestern Canada with some previous studies can be attributed to the presence of a mostly high-velocity anomaly identified by C that extends over 400 km west of Tintina fault and displays localized variations in velocities under ±2%. Although we are unable to offer a specific explanation for C, the Tintina Fault is recognized as a major strike-slip boundary (Price and Carmichael , 1986) along which seismically distinct Phanerozoic mantle lithospheres may have been juxtaposed. Some previous studies, lacking either the resolution or broader geographical context of the present work may thus have misinterpreted this feature as the craton boundary. The results presented here position the boundary between Phanerozoic and cratonic mantle at the Cordilleran deformation front in northwestern Canada and are consistent with the aforementioned heat flow and gravity observations. This is not surprising since coincident velocity transitions from high to low topography and from hot Proterozoic to cold Precambrian mantle have been observed in many other instances across the globe (e.g. Australia (Simons et al., 1999), United States (Dueker et al., 2001), and Northern Europe (Zielhuis and Nolet, 1994)).  84  Chapter 4. Body-wave tomography of western Canada  10 P velocity S velocity Frederiksen 2001 Grand 1997  8  Tintina Fault  6  % of perturbation  4  2  0  −2  −4  −6  −8 −1000  −800  −600  −400  −200  0  200  400  600  800  1000  Distance from the cordilleran deformation front [km]  Figure 4.11: Plot of the velocity variation at 100 km depth along A1-A2 (see Figure 4.10) centered on the Cordillera deformation front. Plain line represents the P velocity perturbations, dashed line the S perturbations, dashed/dotted line the Frederiksen et al. (2001) model, and dotted line Grand et al. (1997).  85  Chapter 4. Body-wave tomography of western Canada  4.4.2  Subsurface geometry of Northern Cascadia  The northern Cascadia subduction zone is a complex tectonic environment characterized by the transition from convergent to transform boundary that involves interaction between the North American continent and three oceanic plates, namely the Juan de Fuca (JdF), Explorer (EXP), and Pacific plates (PAC) (see tectonic setting in Figure 4.12). Although, the structure of Cascadia has been extensively studied in its southern and central portions, its northernmost section remains relatively unexplored due to a historical shortage of geophysical observations in general and seismic data in particular. The geometry of its subsurface is only loosely constrained by scattered heat flow, gravity and geochemical data (Lewis et al., 1997). The recent deployment of broadband three-component seismic stations on northern Vancouver Island, the Nechako area of central British Columbia and along the two arms of the BATHOLITHS array, now permits the study of subsurface structures at the northern edge of the subducting Juan de Fuca plate system and mantle to the north. Juan de Fuca slab and convergence north of Cascadia Subduction Zone The structure of the northern Cascadia region from 48.5◦ N up to 50◦ N has been investigated previously by Bostock and VanDecar (1995) using traveltime tomography. Their velocity model features a quasi-planar, highvelocity body inferred to represent the thermal and compositional anomaly of the subducted Juan de Fuca plate. This body exhibits velocity deviations of up to 3% from the background reference model and extends to depths of at least 400-500 km. Cassidy et al. (1998), using five seismic stations deployed in Northern Vancouver Island and on the adjacent mainland, provided constraints on the shallow S wave velocity structure in Northern Vancouver Island. They observed a pronounced change in the S velocity structure that was interpreted as the northern limit of the subducted oceanic plate beneath Vancouver Island. Using considerably more seismic stations Audet et al. (2008) focused on the structure of the subducted Explorer plate 86  Chapter 4. Body-wave tomography of western Canada  58 oN  56 oN  M ilb an  (G)  R a Ilginb Itc ac ow ha hu z N az ko  ke  54 oN  52 oN  (F) (D) WB  50 oN  PAC 48 oN  (E)  EXP  P velocity perturbation [%] 2 0 2  JdF o  132 W  128oW  o  124 W  o 120 W  o 116 W  Figure 4.12: Depth slice at 100 km depth within the P velocity model in southwestern British Columbia with tectonics of northern Cascadia detailed. Again Blue and red represent positive and negative velocity perturbations, respectively.  87  Chapter 4. Body-wave tomography of western Canada beneath Vancouver Island. Based on receiver function and tomography results, these authors proposed a tectonic model in which the separation of the Explorer plate from the Juan de Fuca plate is described as the consequence of thermo-mechanical erosion of the slab edge and slab thinning at shallow levels and where the slow convergence of the Explorer plate with North America will eventually lead to its capture by the Pacific plate. Bostock and VanDecar (1995), Cassidy et al. (1998), and Audet et al. (2008) did not, however, explore the region beyond the northern tip of Vancouver Island. Figure 4.12 presents an expanded view of our velocity model along the British Columbia coast for the P -velocity model shown in Figure 4.8 at a depth of 100 km. On this figure, four features are identified as D, E, F, and G. The structure D represents a high-velocity anomaly that lies south-east of the Queen Charlotte Islands. This structure may represent the subducted downdip portion of oceanic lithosphere formerly part the Pacific plate that detached from the slab during the formation, at ∼1 Ma, of the Winona block (WB) (Davis and Riddihough, 1982). The high velocity anomaly E is easier to interpret and represents the previously documented thermal and compositional signature of the Juan de Fuca slab (Bostock and VanDecar , 1995). The signature of the slab represents a velocity anomaly of ∼3% in the south and less than 2% in the north. Inland, its northwestern terminus coincides roughly with the projection of the Juan de Fuca displacement vector of the interpreted slab edge beneath Vancouver island documented by Lewis et al. (1997), Cassidy et al. (1998), and Audet et al. (2008) (see dashed line on Figure 4.12). The low-velocity anomaly identified by F is visible beneath part of the northern Cordillera and may be the expression of the slab gap that formed at 45 Ma as a result of the cessation of the subduction of the Kula/Resurrection plate along the northern Canadian margin (presently Queen Charlotte fault). G is a relatively small high-velocity anomaly located east of the Queen Charlotte Basin that may represent the signature of the Pacific plate.  88  Depth [km]  (G)  D1 132 oW o 124 W  (E)  B2 (52.3N, 120W)  B1  B2  0  A2  600 D1 (52.3N, 132W)  (G)  600 C1 (51.5N, 129W)  300  200  100  0  A2 (51.8N, 117W)  500  C2 (54N, 124.5W)  600 A1 (48N, 124W)  Depth [km]  500  )  400  o 116 W  )  400  300  200  100  o 120 W  500  400  300  200  100  0  (E  (E )  89  D2 (55N, 127W)  Figure 4.13: Vertical cross sections within the P velocity model in the northern Cascadia region.  128oW  % of P velocity −2 0 2  D2 (D)  (F)  Depth [km]  (D  600 B1 (49N, 126W)  500  400  300  200  100  0  48 oN  50 oN  52 oN  54 oN  C2 A1  56 oN  C1  Location of cross−section  Depth [km]  58 oN  Chapter 4. Body-wave tomography of western Canada  Chapter 4. Body-wave tomography of western Canada The slices A1–A2 and B1–B2 in Figure 4.13, present two cross sections perpendicular to the strike of Juan de Fuca subduction. On both images, we observe a north-eastward dipping high-velocity anomaly up to 3% that represents the signature of the subducting Juan de Fuca slab. In the south, cross section (A1–A2), the geometry of the slab is well defined and continuous. It dips at a ∼60◦ angle and is traceable to a depth of at least 400-500 km although some smearing may be present. This result is in agreement with the findings of Bostock and VanDecar (1995). In the north, the anomaly is not as clearly delineated and is only defined to depths of ∼300 km. It dips, however, at a similar angle. The degradation of the slab signature northward may be attributed to the lower station density over this region, or alternatively, could result from the action of the erosion mechanisms at the northern edge of the subducting plate. From the cross section C1–C2 of Figure 4.13, we observe the signature of a dipping structure that may reflect plate convergence north of the Cascadia subduction zone. This north-eastward dipping anomaly represents a 2% perturbation that is evident to a depth ∼200 km. As discussed previously, it may represent the signature of formerly subducting material that detached from the Pacific plate, 1 Ma ago, during the formation of the Winona block. It is less clear however, whether there is convergence farther north along the Queen Charlotte Fault. The Queen Charlotte Fault is known to be principally a transform boundary that accommodates ∼50 mm/a of strike-slip motion. It also contains, however, a significant component (∼15 mm/a) of convergence across the margin (e.g. DeMets and Dixon (1999); Kreemer et al. (2003); Mazzotti et al. (2003)). Two models have been proposed to accommodate the ∼15 mm/a of convergence. The first model, involves deformation and shortening within the Queen Charlotte Basin and the Pacific and North American plates (e.g. Rohr et al. (2000)). The second model, involves underthrusting of a Pacific slab beneath the North America Plate (e.g. (Yorath and Hyndman, 1983; Hyndman and Hamilton, 1993; Bustin et al., 2007)). The presence of a high-velocity anomaly centered at 53.5N, 130W, labeledG, within our velocity model may suggest that convergence exists. However, on the cross section D1–D2 of Figure 4.13, this anomaly 90  Chapter 4. Body-wave tomography of western Canada does not share the dipping geometry of the subducting slabs evident in the south (i.e. feature D and F ). Anahim volcanic belt and the Chilcotin group We now turn our attention to the southern edge of the low-velocity zone F, just north of the subducting plate identified by E. In this region, upper Miocene-Quaternary alkalic and peralkalic volcanism is observed along the Anahim volcanic belt and within the Chilcotin plateau basalts, alternatively referred to as the Chilcotin group, contrasting with the subduction driven calc-alkalic volcanism found in the nearby Pemberton and Garibaldi volcanic belt (Bevier , 1989). The Anahim volcanic belt is a 600 km long westeast trending chain of volcanoes that began to erupt at ∼14.5 Ma near the coast and progressed inland to the east, at a rate that varies from 20 to 33 mm/a (Bevier et al., 1979). This belt is composed of five major volcanic centers: The Milbanke Sound group (52.5◦ N 128.7◦ W, ∼14.5 Ma), the Rainbow Range (52.7◦ N 125.8◦ W, ∼8 Ma), the Ilgachuz Range (52.8◦ N 125.3◦ W, ∼5 Ma), the Itcha Range (52.7◦ N 124.9◦ W, ∼2.5 Ma) and the Nazko Range (52.9◦ N, 123.7◦ W, 0.01-0.34 Ma) which is the site of more recent volcanic (Bevier et al., 1979) and seismic activity. Magmatism in the Chilcotin Plateau Basalts, spanned about 16 million years, and occurred in the same area during three main episodes: 15-13 Ma, 9-6 Ma, and 31 Ma (Bevier , 1983b; Mathews, 1989). Although possibly fortuitous, the appearance of volcanic activity in the Anahim belt and Chilcotin group at ∼15 Ma also coincides with an acceleration of uplift rate in the central Coast Mountains (52◦ -54◦ ), from 0.1–0.2 km/Ma to 0.4 km/Ma (Parrish, 1983). Several mechanisms have been proposed to explain the origin of magmatism in this region. Rohr and Currie (1997) suggested that simple lithosphericscale shear could account simultaneously for the subsidence of the Queen Charlotte Basin, the uplift of the central Coast mountains, and the Neogene basalt flow and plutonism in the Anahim belt and Chilcotin group. Asthenospheric flow along the northern edge of the subducting slab was also hypothesized as a plausible cause for the magmatism (Stacey, 1974; Souther  91  Chapter 4. Body-wave tomography of western Canada et al., 1987). The presence of a mantle plume has also been invoked as a potential source of magmatism (Bevier , 1983a; Souther et al., 1987). This idea is supported by the observation that the 20-30 mm/a eastward migration of volcanism in the Anahim belt compares well with the 27 mm/a of the Yellowstone (Wyoming) and Raton hotspots (New Mexico), and is consistent with the predicted 25.1 mm/a rate of movement of the North American plate at the center of the volcanic belt (52◦ N, 125◦ W) (Minster et al., 1974). The peralkalic composition of the Anahim and Chilcotin basalts are also consistent with the plume origin (Bevier , 1989). Figure 4.14 presents cross sections perpendicular and parallel to the inferred hotspot track. On the north-south cross section F1-F2, we observe that the low-velocity anomaly is centered near Nazko Cone and extends to a depth of at least 400 km. The ∼-2% anomaly is flanked on both sides by high-velocity anomalies of variable amplitude. In the south, we recognise the signature of the Juan de Fuca subducting slab complex. High velocities in the north may reflect the remnant of batholithic roots that formed as a result of continuous subduction along the northern margin from ∼150 Ma to ∼50 Ma. On the west-east cross section E1-E2, the low-velocity anomaly extends from the current location of the Anahim hotspot at Nazko cone seaward along the hotspot track. East of Nazko cone, we observe a small high-velocity anomaly marking the eastern extent of the hotspot track. In the west, we recognize the signature of the high-velocity anomaly labeled D in Figure 4.12. From the images in top and middle panel of Figures 4.14 we observe that the low-velocity anomaly is visible to a depth of ∼400 km directly beneath Nazko Cone. As suggested by our model, the low-velocity anomaly may, however, extend deeper southward beneath the Juan de Fuca plate. Low-velocity anomalies extending beneath hotspots, or hotspot tracks have been imaged by regional traveltime tomography in several other locations around the globe (e.g.  Hawaii (Ellsworth and Koyanagi , 1977;  Tilmann, 1999), Yellowstone (Iyer et al., 1981; Evans, 1982; Saltzer and Humphreys, 1997), the Massif Central (Granet et al., 1995), Iceland (Wolfe et al., 1997), and in Germany (Ritter et al., 2001)). The general consensus is 92  Chapter 4. Body-wave tomography of western Canada  F2  Location of cross−section  56 oN  54 oN  nk e R a Ilginb Itc ac o h hw N a uz az ko  E2  E1  52 oN  M  ilb a  (G)  (F) (D)  50 oN  (E)  % of P velocity −2  0  2 F1  48 oN 136 oW  132 oW  128oW  o 124 W  o 120 W  o 116 W  Nazko Cone  0  (D)  (F)  100  Depth [km]  200  300  400  500  600  700  E1 (52.3N, 131.1W)  E2 (53.5N, 116.4W)  Nazko Cone  0  (F)  100  (E)  Depth [km]  200  300  400  500  600  700  F1 (48.4N, 123.4W)  F2 (57.4N, 124.1W)  Figure 4.14: Vertical cross sections within the P velocity model.  93  Chapter 4. Body-wave tomography of western Canada that these anomalies are related to or are the source of volcanism observed at the surface. By analogy, the presence in our tomographic model of a lowvelocity anomaly present beneath Nazko Cone that extends seaward along the Anahim volcanic track suggests a direct relation to the magmatism observed in this region. Our results, consequently, support models of volcanism that involve mantle-scale rather that lithosphere-scale processes (i.e. simple lithospheric-scale shear Rohr and Currie (1997)), but cannot distinguish between asthenospheric flow along the slab edge (Stacey, 1974; Souther et al., 1987) and mantle plume models (Bevier , 1983a; Souther et al., 1987). Upwelling of asthenosphere along the slab edge is postulated to exist in slab window environments where a pair of diverging plates is subducted beneath a continent (e.g. (Dickinson and Snyder , 1979; Johnson and O’Neil , 1984; Forsythe and Nelson, 1985; Thorkelson and Taylor , 1989; Thorkelson, 1996)). This concept could perhaps be applied to the current tectonic setting where the gap left north of the Juan de Fuca plate, following the end of subduction along what is now the Queen Charlotte fault, may have focussed upwelling asthenosphere toward the Anahim volcanic belt region. Shear-wave splitting observations in the southern Cascadia region are suggestive of such complex flow along the edge of subducting plates (Zandt and Humphreys, 2008). The sudden appearance of magmatism at ∼15 Ma and the steady spatio-temporal eastward progression of the hotspot along an axis coherent with the displacement of the North American plate are, however, less readily accounted for by this model. The obliquity of the Juan de Fuca plate motion to the hotspot track also poses a geometrical problem. The presence of a mantle plume either confined to the upper-mantle, or extending beneath the subducting plate below the transition zone offers a simpler explanation. Mantle plumes are known to produce stable convection during an extended period of time (e.g. Nataf (2000)) and are recognized to be the source of hotspot volcanism and perhaps also the source of flood basalts in several similar locations (Richards et al., 1989; Duncan and Richards, 1991).  94  Chapter 4. Body-wave tomography of western Canada  4.5  Conclusion  In this paper, we have presented the first regional-scale, high resolution P and S velocity models across western Canada, utilizing data collected at stations from a variety of temporary and permanent seismic networks. Our models indicate that the western portion of the study area displays a generally lower velocity, whereas the eastern region is mostly fast. The only significant exception is the presence of a broad high-velocity anomaly in the Cascadia subduction zone. In this study we focused our attention on two areas where our station distribution affords the best resolution, namely, the corridor along legs A and B of the CANOE array in Northwestern Canada, and the northern margin of Cascadia where the plate boundary changes from convergent to transform. In northwestern Canada, we have characterized the transition from Cordilleran to cratonic mantle. We argue that the main change occurs at the Cordilleran deformation front and represents a sharp increase in mantle seismic velocity greater than 4% over a distance of ∼50 km. At the northern margin of Cascadia, we have imaged three significant velocity anomalies: a high-velocity anomaly centered at 51◦ N, 128◦ W, dipping eastward, possibly representing the remnant of the slowly converging Winona block; a high-velocity anomaly representing the thermal signature of the Juan de Fuca subducting plate; and a low-velocity anomaly north of the slab. Immediately north of the northern edge of the slab, we have paid special attention to mantle structure along and perpendicular to the Anahim hotspot track. In this region, our model reveals a low-velocity zone north of the Juan de Fuca plate that extends well into mantle depth beneath Nazko Cone and may continue southward beneath the Juan de Fuca plate. Our observations suggest that a deep-seated low-velocity anomaly is the source of magmatism in the Anahim volcanic belt and Chilcotin Basalt group, and, consequently, that a mantle-scale process rather than lithospheric scale process controls surface volcanism. Coupled with the well documented temporal progression of surface volcanism, our model favours an origin for the hotspot track in the form of a mantle plume over slab edge flow.  95  Bibliography Audet, P., A. M. Jellinek, and H. Uno, Mechanical controls on the deformation of continents at convergent margins, Earth Planet. Sc. Lett., 264 (1-2), 151–166, doi:{10.1016/j.epsl.2007.09.024}, 2007. Audet, P., M. G. Bostock, and J.-P. Mercier, Morphology of the Explorer/Juan de Fuca slab edge in northern Cascadia: Imaging plate capture at a ridge-trench-transform triple junction, Submitted to Geology, 2008. Bevier, M., Implications of chemical and isotopic composition for petrogenesis of Chilcotin Group basalts, British Columbia, J. Petrol, 24, 207–226, 1983a. Bevier, M., Regional stratigraphy and age of Chilcotin Group basalts, south-central British Columbia, Can. J. Earth Sci., 20 (4), 515–524, 1983b. Bevier, M., R. Armstrong, and J. Souther, Miocene peralkaline volcanism in West-central British Columbia; its temporal and plate-tectonics setting, Geology, 7 (8), 389–392, 1979. Bevier, M. L., A lead and strontium isotopic study of the Anahim volcanic belt, British Columbia; additional evidence for widespread suboceanic mantle beneath western North America, Geol. Soc. Am. Bull., 101 (7), 973–981, 1989. Bostock, M. G., and J. C. VanDecar, Upper-mantle structure of the northern cascadia subduction zone, Can. J. Earth Sci., 32 (1), 1–12, 1995. Buchbinder, G., and G. Poupinet, P-wave residuals in Canada, Can. J. Earth Sci., 14 (6), 1292–1304, 1977. Burdick, S., et al., Upper Mantle Heterogeneity beneath North America from Travel Time Tomography with Global and USArray Transportable Array Data, Seismological Research Letters, 79 (3), 384, 2008.  96  Bibliography Bustin, A., R. Hyndman, H. Kao, and J. Cassidy, Evidence for underthrusting beneath the Queen Charlotte Margin, British Columbia, from teleseismic receiver function analysis, Geophys. J. Int., 171 (3), 1198–1211, 2007. Cassidy, J., R. Ellis, C. Karavas, and G. Rogers, The northern limit of the subducted Juan de Fuca plate system, J. Geophys. Res., 103 (B11), 26,949–26,961, 1998. Cook, F., R. Clowes, D. Snyder, A. van der Velden, K. Hall, P. Erdmer, and C. Evenchick, Precambrian crust beneath the Mesozoic northern Canadian Cordillera discovered by Lithoprobe seismic reflection profiling, Tectonics, 23 (2), 2004. Davis, E., and R. Riddihough, The Winona Basin: structure and tectonics, Can. J. Earth Sci., 19 (4), 767–788, 1982. DeMets, C., and T. Dixon, New kinematic models for Pacific-North America motion from 3 Ma to present, I: Evidence for steady motion and biases in the NUVEL-1A model, Geophys. Res. Lett, 26 (13), 1921–1924, 1999. Dickinson, W., and W. Snyder, Geometry of subducted slabs related to San Andreas transform, Journal of Geology, 87 (6), 609–627, 1979. Dueker, K., H. Yuan, and B. Zurek, Thick-Structured Proterozoic Lithosphere of the Rocky Mountain Region, GSA Today, 11 (12), 4–9, 2001. Duncan, R., and M. Richards, Hotspots, mantle plumes, flood basalts, and true polar wander, Reviews of Geophysics, 29 (1), 31–50, 1991. Ekstr¨ om, G., J. Tromp, and E. Larson, Measurements and global models of surface wave propagation, J. Geophys. Res., 102 (B4), 8137–8158, doi: 10.1029/96JB03729, 1997. Ellsworth, W., and R. Koyanagi, Three-dimensional crust and mantle structure of Kilauea Volcano, Hawaii, J. Geophys. Res, 82 (33), 5379–5394, 1977. Evans, J., Compressional wave velocity structure of the upper 350 km under the eastern Snake River Plain near Rexburg, Idaho, J. Geophys. Res., 87 (B4), 2654–2670, 1982. Fl¨ uck, P., R. Hyndman, and C. Lowe, Effective elastic thickness T e of the lithosphere in western Canada, J. Geophys. Res., 108 (B9), 2430, doi: 10.1029/2002JB002201, 2003. 97  Bibliography Forsythe, R., and E. Nelson, Geological manifestations of ridge collision: evidence from the Golfo de Penas-Taitao Basin, southern Chile, Tectonics, 4 (5), 477–495, 1985. Fouch, M., D. James, J. VanDecar, and S. van der Lee, Mantle seismic structure beneath the Kaapvaal and Zimbabwe Cratons, S. Afr. J. Geol., 107 (1-2), 33–44, doi:10.2113/107.1-2.33, 2004. Frederiksen, A. W., M. G. Bostock, J. C. VanDecar, and J. F. Cassidy, Seismic structure of the upper mantle beneath the northern Canadian Cordillera from teleseismic travel-time inversion, Tectonophysics, 294 (12), 43–55, 1998. Frederiksen, A. W., M. G. Bostock, and J. F. Cassidy, S-wave velocity structure of the Canadian upper mantle, Phys. Earth Planet. In., 124 (34), 175–191, 2001. Grand, S., V. der Hilst R.D., and W. S., High resolution global tomography: a snapshot of convection in the Earth, GSA Today, 7, 1–7, 1997. Grand, S. P., Mantle shear structure beneath the america and surrounding oceans, J. Geophys. Res., 99 (B6), 11,591–11,621, 1994. Granet, M., M. Wilson, and U. Achauer, Imaging a mantle plume beneath the French Massif Central, Earth Planet. Sc. Lett., 136 (3-4), 281–296, 1995. Huber, P., Robust Statistics, 1981. Humphreys, E. D., and K. G. Dueker, Western united-states upper-mantle structure, J. Geophys. Res., 99 (B5), 9615–9634, 1994. Hyndman, R., and T. Hamilton, Queen Charlotte area Cenozoic tectonics and volcanism and their association with relative plate motions along the northeastern Pacific margin, J. Geophys. Res., 98 (B8), 14–257, 1993. Hyndman, R., and T. Lewis, Geophysical consequences of the Cordillera– Craton thermal transition in southwestern Canada, Tectonophysics, 306 (34), 397–422, doi:10.1016/S0040-1951(99)00068-2, 1999. Hyndman, R., C. Currie, and S. Mazzotti, Subduction zone backarcs, mobile belts, and orogenic heat, GSA Today, 15 (2), 4–10, 2005. Iyer, H. M., J. R. Evans, G. Zandt, R. M. Stewart, J. M. Coakley, and J. N. Roloff, A deep low-velocity body under the Yellowstone caldera, Wyoming: 98  Bibliography Delineation using teleseismic P-wave residuals and tectonic interpretation: Summary, B. Seismol. Soc. Am., 92 (11), 792–798, 1981. Johnson, C., and J. O’Neil, Triple junction magmatism: a geochemical study of Neogene volcanic rocks in western California, Earth Planet. Sc. Lett., 71 (2), 241–262, 1984. Jones, A., Imaging the continental upper mantle using electromagnetic methods, Lithos, 48 (1-4), 57–80, doi:10.1016/S0024-4937(99)00022-5, 1999. Jones, A., et al., The electrical resistivity structure of Archean to Tertiary lithosphere along 3200 km of SNORCLE profiles, northwestern Canada, Can. J. Earth Sci, 42 (6), 1257–1275, 2005. Jordan, T., Composition and development of the continental tectosphere, Nature, 274 (5671), 544–548, doi:10.1038/274544a0, 1978. Kennett, B. L. N., and E. R. Engdahl, Traveltimes for global earthquake location and phase identification, Geophys. J. Int., 105 (2), 429–465, 1991. Kreemer, C., W. Holt, and A. Haines, An integrated global model of present-day plate motions and plate boundary deformation, Geophys. J. Int., 154, 8–34, doi:10.1046/j.1365-246X.2003.01917.x, 2003. Lachenbruch, A., and P. Morgan, Continental extension, magmatism and elevation; formal relations and rules of thumb, Tectonophysics, 174 (1-2), 39–62, 1990. Ledo, J., A. Jones, and I. Ferguson, Electromagnetic images of a strike-slip fault- The Tintina fault Northern Canadian, Geophys. Res. Lett., 29 (8), 66–1, doi:10.1029/2001GL013408, 2002. Lewis, T., C. Lowe, and T. Hamilton, Continental signature of a ridgetrench-triple junction: Northern Vancouver Island, J. Geophys. Res., 102 (B4), 7767–7781, 1997. Lewis, T., R. Hyndman, and P. Fl¨ uck, Heat flow, heat generation, and crustal temperatures in the northern Canadian Cordillera: thermal control of tectonics, J. Geophys. Res., 108 (B6), 2316, doi:10.1029/2002JB002090, 2003. Mathews, W., Neogene Chilcotin basalts in south-central British Columbia: geology, ages, and geomorphic history, Can. J. Earth Sci., 26 (5), 969–982, 1989. 99  Bibliography Mazzotti, S., H. Dragert, J. Henton, M. Schmidt, R. Hyndman, T. James, Y. Lu, and M. Craymer, Current tectonics of northern Cascadia from a decade of GPS measurements, J. Geophys. Res, 108, 2554, 2003. Minster, J., T. Jordan, P. Molnar, and E. Haines, Numerical Modelling of Instantaneous Plate Tectonics, Geophys. J. Int., 36 (3), 541–576, 1974. Nataf, H., Seismic Imaging of Mantle Plumes, Annu. Rev. Earth. Pl. Sc, 28 (1), 391–417, 2000. Parrish, R., Cenozoic thermal evolution and tectonics of the Coast Mountains of British Columbia 1. Fission track dating, apparent uplift rates, and patterns of uplift, Tectonics, 2 (2), 601–632, 1983. Polet, J., and D. Anderson, Depth extent of cratons as inferred from tomographic studies, Geology, 23 (3), 205–208, 1995. Price, R., and D. Carmichael, Geometric test for Late CretaceousPaleogene intracontinental transform faulting in the Canadian Cordillera, Geology, 14 (6), 468–471, 1986. Richards, M., R. Duncan, and V. Courtillot, Flood Basalts and Hot-Spot Tracks: Plume Heads and Tails, Science, 246 (4926), 103–107, 1989. Ritter, J., M. Jordan, U. Christensen, and U. Achauer, A mantle plume below the Eifel volcanic fields, Germany, Earth and Planetary Science Letters, 186 (1), 7–14, 2001. Rohr, K., and L. Currie, Queen Charlotte Basin and Coast Mountains; paired belts of subsidence and uplift caused by a low-angle normal fault, Geology, 25 (9), 819–822, 1997. Rohr, K. M. M., M. Scheidhauer, and A. M. Trehu, Transpression between two warm mafic plates: The Queen Charlotte Fault revisited, J. Geophys. Res., 105 (B4), 8147–8172, 2000. Saltzer, R., and E. Humphreys, Upper mantle P wave velocity structure of the eastern Snake River Plain and its relationship to geodynamic models of the region, J. Geophys. Res., 102 (B6), 11,829–11,842, 1997. Shragge, J., M. G. Bostock, C. G. Bank, and R. M. Ellis, Integrated teleseismic studies of the southern Alberta upper mantle, Can. J. Earth Sci., 39 (3), 399–411, doi:10.1139/E01-084, 2002.  100  Bibliography Sigloch, K., N. McQuarrie, and G. Nolet, Two-stage subduction history under North America inferred from multiple-frequency tomography, Nature Geoscience, 1 (7), 458, 2008. Simons, F., A. Zielhuis, and R. van der Hilst, The deep structure of the Australian continent from surface wave tomography, Lithos, 48 (1-4), 17– 43, doi:10.1016/S0024-4937(99)00041-9, 1999. Souther, J., J. Clague, and R. Mathewes, Nazko cone: a Quaternary volcano in the eastern Anahim Belt, Can. J. Earth Sci., 24 (12), 2477–2485, 1987. Stacey, R. A., Plate tectonics, volcanism and the lithosphere in British Columbia, Nature, 250 (5462), 133–134, doi:10.1038/250133a0, 1974. Thorkelson, D., Subduction of diverging plates and the principles of slab window formation, Tectonophysics, 255 (1-2), 47–63, 1996. Thorkelson, D., and R. Taylor, Cordilleran slab windows, Geology, 17 (9), 833–836, 1989. Tilmann, J. T., The seismic structure of the upper mantle beneath Hawaii, Ph.D. thesis, Camridge, 1999. Van Der Lee, S., and A. W. Frederiksen, Surface-wave tomography applied to the north american upper mantle, Geophysical Monograph, 157, 67–80, 2005. VanDecar, J., Upper-mantle structure of the cascadia subduction zone from non-linear teleseismic travel-time inversion, Ph.D. thesis, University of Washington, Seattle, 1991. VanDecar, J. C., and R. S. Crosson, Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least-squares, B. Seismol. Soc. Am., 80 (1), 150–169, 1990. Wickens, A. J., The upper mantle of southern British Columbia, Can. J. Earth Sci., 14 (5), 1100–1115, 1977. Wickens, A. J., and G. G. R. Buchbinder, S-wave residuals in Canada, B. Seismol. Soc. Am., 70 (3), 809–822, 1980. Wolfe, C., I. Th Bjarnason, J. VanDecar, and S. Solomon, Seismic structure of the Iceland mantle plume, Nature, 385 (6613), 245–247, 1997. 101  Bibliography Yorath, C. J., and R. D. Hyndman, Subsidence and thermal history of Queen Charlotte Basin, Can. J. Earth Sci, 20, 135–159, 1983. Zandt, G., and E. Humphreys, Toroidal mantle flow through the western US slab window, Geology, 36 (4), 295–298, 2008. Zielhuis, A., and G. Nolet, Shear-wave velocity variations in the upper mantle beneath central Europe, Geophys. J. Int., 117 (3), 695–715, doi: 10.1111/j.1365-246X.1994.tb02463.x, 1994.  102  Chapter 5  Conclusion 5.1  Key Contributions  This thesis focused both on the development of an innovative deconvolution scheme and on the understanding of the western Canada lithospheric structure through receiver functions and traveltime tomography studies. In this document we presented and successfully applied a new deconvolution technique that provides a more accurate estimate of the so-called teleseismic P -Green’s function from broadband three-component seismograms than the more traditional receiver function technique. To the best of our knowledge, this new technique offers the first means of extracting meaningful pure P -mode scattering from near-horizontal structures, such as the continental Moho. This thesis also featured a receiver function study that revealed the teleseismic signature of a Proterozoic subduction zone in the Wopmay orogen and provided important insight on the structure and origin of continents. Finally, we presented the first P - and S -velocity models that spanned the entire western Canada from the craton to the coast. These velocity models provided valuable insights on the structure of the uppermantle beneath Western Canada. This study permitted the 1) localization and characterization of the transition from Phanerozoic to Proterozoic in northwestern Canada, 2) the exploration the northern Cascadia region, and the mantle beneath the Anahim hotspot. The next sections present the details of the main contributions included in this thesis.  103  Chapter 5. Conclusion  5.1.1  Improved Green’s functions for passive source structural studies  This chapter presented the application of a innovative deconvolution technique that allows the extraction of novel information from the P -component of the teleseismic Green’s function. It featured results for three stations of the Canadian National Seismic Network, namely, ULM (Lac du Bonnet, Manitoba), WHY (Whitehorse, Yukon) and EDM (Edmonton, Alberta), each representing a different geologic environment. Deconvolution results at each of the stations, showed evidence for the principal, first-order scattered Moho phases and, in particular, the P pM p multiple. Station ULM exhibited a very simple crustal structure that is reasonably well modelled as a layer over a half-space. The crust-mantle boundary at station WHY, as portrayed through PM s and P pM p, was revealed to be rather more complex, representing a thin, high-velocity lid. Intra-crustal scattering dominates the Green’s functions recovered at station EDM. The sedimentary-basement contact produced prominent arrivals at early times followed by signals from a pronounced low-velocity layer that has been previously interpreted from seismic reflection and receiver functions to represent a serpentinized sliver upthrust during Proterozoic subduction.  5.1.2  The teleseismic signature of fossil subduction: Northwestern Canada  This chapter presented receiver functions calculated from data acquired between June 2003 and September 2005, at 20 broadband, three-component seismometers deployed along the MacKenzie-Liard Highway in Canada’s Northwest Territories as part of the joint IRIS-Lithoprobe CANOE transect. Teleseismic receiver functions computed from ∼250 earthquakes clearly reveal the response of the ancient subduction zone. On the radial component, the suture was evident as a direct conversion from the Moho, the depth of which increases from ∼30 km to ∼50 km over a horizontal distance of ∼70 km before its signature disappears. The structure was still better defined on the transverse component where the Moho appears as the up104  Chapter 5. Conclusion per boundary of a 10 km thick layer of anisotropy that can be traced from 30 km to at least 90 km depth. The seismic response of this layer was characterized by a frequency dependence that can be modelled by upper and lower boundaries that are discontinuous in material properties and their gradients, respectively. Anisotropy was characterized by a ±5% variation in shear velocity and hexagonal symmetry with a fast axis that plunges at an oblique angle to the subduction plane. The identification of this structure provided an unambiguous connection between fossil subduction and fine-scale, anisotropic mantle layering. Previous documentation of similar layering below the adjacent Slave province and from a range of Precambrian terranes across the globe provided strong support for the thesis that early cratonic blocks were stabilized through processes of shallow subduction.  5.1.3  Body-wave tomography of western Canada  This chapter presented a tomographic inversion of the teleseismic P and S delay times that took advantage of the deployment of broadband stations, in several recent portable seismic experiment to produce a P and S velocity map of western Canada with improved resolution over previous studies. It presented images of the major large-scale upper-mantle velocity structure beneath the whole western Canada and provided constraints on the structure and evolution of this region. The resolution afforded by our model was best in southwestern British Columbia, and along the CANOE (north-western Alberta, southern Yukon and Northwest Territories) and BATHOLITHS (Northwestern BC) arrays where the station density is the highest, and fair elsewhere. This study focused on three distinct features which are the transition from Phanerozoic to cratonic mantle in northwestern Canada, the complex tectonic environment at the northern terminus of the Cascadia subduction zone where the plate boundary changes from convergent to transform, and the Ahanim hotspot. This study revealed 1) that the main transition from Phanerozoic to cratonic mantle in northwestern Canada occurs at the Cordilleran deformation front over a distance of ∼50 km 2) evidence of the presence of subducted material beyond the north-  105  Chapter 5. Conclusion ern edge of the slab 3) that a-2% low-velocity zone underlies the Anahim volcanic belt reaching a depth of ∼400 km beneath Nazko cone.  5.2  Limitations and Future work  The limitations of the work presented in this documents are detailed project by project in the next sections. Future directions to mitigate the limitations encountered or research opportunities are also presented.  5.2.1  Improved Green’s functions for passive source structural studies  The main limitation of the new deconvolution technique presented in this thesis concerns the estimation of the source spectrum. The technique currently implemented utilizes a source spectrum estimate that is basically a smoothed version of the seismogram amplitude spectrum.  Although,  this method is successful in extracting significant information from the P component of the teleseismic Green’s function, significant arrivals are consistently missing from estimates (e.g. P sM p phase see Figure 2.1). To mitigate this problem, a better source estimation scheme has to be developed and implemented in the existing framework. The recent breakthroughs in the field of large polynomial factorization and in the development of robust optimization methods suitable for this problem may provide a means of obtaining a more accurate estimate of the source, and consequently of the teleseismic Green’s function.  5.2.2  The teleseismic signature of fossil subduction: Northwestern Canada  A central objective of this project was to establish the relationship between the fossil subduction found beneath the Wopmay orogen and layers (H, X and L) documented by Bostock and Cassidy (1997); Bostock (1998) beneath the Slave craton. As previously discussed, the signature of layers documented by Bostock (1998) and that of AM bear a close similarity and sug106  Chapter 5. Conclusion gest common, generic origin related to plate convergence episodes. However, it is unfortunate that our teleseismic profiles were not capable to establish a clear geometrical (continuous or discontinuous) relationship between AM and H, X or L due to the highly variable and intermittent characteristics of the signals at stations A16 and A17. Also, although we clearly identified AM with paleosubduction, the limited observation of similar structures precludes us from finding a conclusive mechanism for its formation. More research will have to be carried out to fully comprehend the processes involved in the formation of AM.  5.2.3  Body-wave tomography of western Canada  The unequal distribution of stations across the study area was the main limitation we encountered throughout this project. The station density varied widely between the different region. The station density was high along CANOE, BATHOLITHS array and in southwestern British Columbia and low to null elsewhere. Several region in central British Columbia and along the northern coast, did not afford sufficient resolution and, consequently, were not studied. Another important limitation concerns the inversion method adopted in this study. The VanDecar (1991) technique while producing good estimates of the spatial variation within the studied area, cannot provide reliable estimates of absolute seismic velocities. The quasi-vertical nature of the body-wave ray paths considerably reduced our capability to accurately estimate the depth extent of the observed anomalies; tomography consistently suffers from significant amount of smearing along the ray paths.  107  Bibliography Bostock, M. G., Mantle stratigraphy and evolution of the Slave province, J. Geophys. Res., 103, 183–200, 1998. Bostock, M. G., and J. F. Cassidy, Upper mantle stratigraphy beneath the southern Slave craton, Can. J. Earth Sci., 34, 577–587, 1997. VanDecar, J., Upper-mantle structure of the cascadia subduction zone from non-linear teleseismic travel-time inversion, Ph.D. thesis, University of Washington, Seattle, 1991.  108  Appendix A  Tomography Resolution Test This appendix presents the input and output of the resolution test for results shown in Figures 4.8 to 4.14 of Chapter 4. The resolution afforded by the data through the inversion procedure detailed Chapter 4 was assessed using a standard checkerboard model. synthetic traveltimes were produced for the event-station ray geometry of the real data and inverted using the inversion scheme and parameters described in sections 4.2.3 and 4.2.4.  109  −2  −1  % of P velocity 136 oW  tin au lt 128 W  aF o  nt  an ro ler n f rdil atio Co form de  120 W  o  o 112 W  144 o W  45 oN  50 oN  N  65 o  −2  −2  −1  0  1  2  60 oN  55 oN  144 o W 136 oW  % of P velocity 136 oW  400 km  128oW  tin au  aF lt 128oW  o 120 W  nt o 120 W  an ro ler n f rdil atio Co form de  110  o 112 W  o 112 W  Figure A.1: Input P resolution test model for the depth slices presented in Figure 4.8.  144 o W  45 oN  0  1  Tin  2  300 km  o 120 W  nt an ro ler n f rdil atio Co form de  60 oN  lt au 128oW  o 112 W  −1  0  % of P velocity  Tin  50 oN  N  % of P velocity 136 oW  45 oN  1  2  lt au  65 o  −2  −1  50 oN  200 km  nt an ro ler n f rdil atio Co form de  55 oN  144 o W  45 oN  0  aF tin Tin  1  2  55 oN  60 oN  65 oN  aF tin Tin  50 oN  55 oN  60 oN  65 oN  100 km  Appendix A. Tomography Resolution Test  136 oW  tin au lt 128 W  aF o  nt  an ro ler n f rdil atio Co form de  120 W  o  o 112 W  144 o W  45 oN  50 oN  −1 136 oW  % of P velocity −0.5  0  0.5  1  136 oW  400 km  128oW  tin au  aF lt 128oW  o 120 W  nt o 120 W  an ro ler n f rdil atio Co form de  111  o 112 W  o 112 W  Figure A.2: Output of the resolution test for the P data-set for the depth slices presented in Figure 4.8.  −0.5  % of P velocity  144 o W  45 oN  0  0.5  Tin  50 oN  N  65 o  −1  60 oN  55 oN  144 o W  Tin  55 oN  300 km  o 120 W  nt an ro ler n f rdil atio Co form de  60 oN  lt au 128oW  o 112 W  % of P velocity  −0.5  0  0.5  1  lt au  N  136 oW  45 oN  50 oN  200 km  nt an ro ler n f rdil atio Co form de  65 o  −0.5  % of P velocity  144 o W  45 oN  aF tin Tin  0  0.5  55 oN  60 oN  65 oN  aF tin Tin  50 oN  55 oN  60 oN  65 oN  100 km  Appendix A. Tomography Resolution Test  −2  −1  % of S velocity 136 oW  tin au lt 128 W  aF o  nt  an ro ler n f rdil atio Co form de  120 W  o  o 112 W  144 o W  45 oN  50 oN  N  65 o  −2  −2  −1  0  1  2  60 oN  55 oN  144 o W 136 oW  % of S velocity 136 oW  400 km  128oW  tin au  aF lt 128oW  o 120 W  nt o 120 W  an ro ler n f rdil atio Co form de  112  o 112 W  o 112 W  Figure A.3: Input S resolution test model for the depth slices presented in Figure 4.9.  144 o W  45 oN  0  1  Tin  2  300 km  o 120 W  nt an ro ler n f rdil atio Co form de  60 oN  lt au 128oW  o 112 W  −1  0  % of S velocity  Tin  50 oN  N  % of S velocity 136 oW  45 oN  1  2  lt au  65 o  −2  −1  50 oN  200 km  nt an ro ler n f rdil atio Co form de  55 oN  144 o W  45 oN  0  aF tin Tin  1  2  55 oN  60 oN  65 oN  aF tin Tin  50 oN  55 oN  60 oN  65 oN  100 km  Appendix A. Tomography Resolution Test  −1 136 oW  % of S velocity  −0.5  tin au lt 128 W  aF o  nt  an ro ler n f rdil atio Co form de  120 W  o  o 112 W  144 o W  45 oN  50 oN  −1  136 oW  136 oW  % of S velocity −0.5  0  0.5  1  60 oN  55 oN  N  65 o  −1  400 km  128oW  tin au  aF lt 128oW  o 120 W  nt o 120 W  an ro ler n f rdil atio Co form de  113  o 112 W  o 112 W  Figure A.4: Output of the resolution test for the S data-set for the depth slices presented in Figure 4.9.  144 o W  45 oN  0  0.5  Tin  1  60 oN  300 km  nt an ro ler n f rdil atio Co form de o 120 W  144 o W  −0.5  0  % of S velocity  Tin  50 oN  N  lt au 128oW  o 112 W  45 oN  0.5  1  lt au  65 o  −1 136 oW  % of S velocity  −0.5  50 oN  200 km  nt an ro ler n f rdil atio Co form de  55 oN  144 o W  45 oN  0  aF tin Tin  0.5  1  55 oN  60 oN  65 oN  aF tin Tin  50 oN  55 oN  60 oN  65 oN  100 km  Appendix A. Tomography Resolution Test  Appendix A. Tomography Resolution Test  Tintina Fault  64 oN 62 oN 60 oN A1 58 oN 56 oN 54 oN 144 o W 13 o 6 W  0 100  ]  200 300 400  A2  o W o 128oW 120 W 112  3 500  1.5  600  0  700 A1 (60.5N, 138.7W)  A2 (60.5N, 108.7W)  −1.5  % of P velocity  Depth [km  Cordilleran front  −3  Figure A.5: Input P resolution test model for the cross sections presented in Figure 4.10.  Tintina Fault  0 100  300 400  64 oN 62 oN 60 oN A1 58 oN 56 oN 54 oN 144 o W 13 o 6 W  A2  o W o 128oW 120 W 112  3 500  1.5  600 700 A1 (60.5N, 138.7W)  0 A2 (60.5N, 108.7W)  −1.5  % of P velocity  ] Depth [km  200  Cordilleran front  −3  Figure A.6: Output of the resolution test for the S data-set for the depth slices presented in Figure 4.10.  114  115  600 B1 (49N, 126W)  500  400  300  200  100  B2 (52.3N, 120W)  Depth [km]  500  600 D1 (52.3N, 132W)  500  600 C1 (51.5N, 129W)  C2 (54N, 124.5W)  400  300  200  100  0  A2 (51.8N, 117W)  400  300  200  100  0  600 A1 (48N, 124W)  500  400  300  200  100  Depth [km]  D2 (55N, 127W)  Figure A.7: Input P resolution test model for the cross sections presented in Figure 4.13.  Depth [km]  0  Depth [km]  0  Appendix A. Tomography Resolution Test  116  600 B1 (49N, 126W)  500  400  300  200  100  B2 (52.3N, 120W)  Depth [km]  500  600 D1 (52.3N, 132W)  500  600 C1 (51.5N, 129W)  C2 (54N, 124.5W)  400  300  200  100  0  A2 (51.8N, 117W)  400  300  200  100  0  600 A1 (48N, 124W)  500  400  300  200  100  Depth [km]  D2 (55N, 127W)  Figure A.8: Output of the resolution test for the P data-set for the depth slices presented in Figure 4.13.  Depth [km]  0  Depth [km]  0  Appendix A. Tomography Resolution Test  Appendix A. Tomography Resolution Test  F2  Location of cross−section  56 oN  54 oN  E1  52 oN  M  ilb  an ke R a I in Itlcgacbo h hw N a uz az ko  E2  50 oN  % of P velocity −2  0  2 F1  48 oN 136 oW  132 oW  128oW  o 124 W  o 120 W  o 116 W  Nazko Cone  0  100  Depth [km]  200  300  400  500  600  700  E1 (52.3N, 131.1W)  E2 (53.5N, 116.4W)  Nazko Cone  0  100  Depth [km]  200  300  400  500  600  700  F1 (48.4N, 123.4W)  F2 (57.4N, 124.1W)  Figure A.9: Input P resolution test model for the cross sections presented in Figure 4.14.  117  Appendix A. Tomography Resolution Test  F2  Location of cross−section  56 oN  54 oN  M  ilb  an ke R a I in Itlcgacbo h hw N a uz az ko  E2  E1  52 oN  50 oN  % of P velocity −2  0  2 F1  48 oN 136 oW  132 oW  128oW  o 124 W  o 120 W  o 116 W  Nazko Cone  0  100  Depth [km]  200  300  400  500  600  700  E1 (52.3N, 131.1W)  E2 (53.5N, 116.4W)  Nazko Cone  0  100  Depth [km]  200  300  400  500  600  700  F1 (48.4N, 123.4W)  F2 (57.4N, 124.1W)  Figure A.10: Output of the resolution test for the P data-set for the depth slices presented in Figure 4.14.  118  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052560/manifest

Comment

Related Items