Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Mercier, Jean-Philippe
Doctor of Philosophy - PhD
University of British Columbia
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 Pp_Mp. 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 Lithoprobe-IRIS 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 Cascadia 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.
23010851 bytes
application/pdf
Improved teleseismic Green’s functions and western Canada mantle structure and evolution by Jean-Philippe Mercier B.Eng., École Polytechnique de Montréal, 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 PpMp. The second part presents teleseismic receiver functions from 20 broadband three-component seismometers deployed along the MacKenzie-Liard High- way in Canada’s Northwest Territories as part of the joint Lithoprobe- IRIS 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 unambigu- ous connection between fossil subduction and fine-scale, anisotropic mantle layering found beneath cratons. Previous documentation of similar layer- ing below the adjacent Slave province provides strong support for the thesis that early cratonic blocks were stabilized through processes of shallow sub- duction. 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 Casca- ii 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 . . . . . . . . . . . . . . 12 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Station ULM (Lac du Bonnet, Manitoba) . . . . . . . 18 2.3.2 Station WHY (Whitehorse, Yukon) . . . . . . . . . . 22 iv Table of Contents 2.3.3 Station EDM (Edmonton, Alberta) . . . . . . . . . . 22 2.4 Discussion and Concluding Remarks . . . . . . . . . . . . . . 24 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3 The teleseismic signature of fossil subduction: Northwestern Canada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Tectonic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.6 Modelling and Numerical Simulation . . . . . . . . . . . . . 42 3.6.1 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.6.2 One-way modeling . . . . . . . . . . . . . . . . . . . . 45 3.7 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 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Data and Method . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2.2 Model parametrization . . . . . . . . . . . . . . . . . 70 4.2.3 Traveltime inversion . . . . . . . . . . . . . . . . . . . 70 4.2.4 Choice of inversion parameters . . . . . . . . . . . . . 72 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3.1 Resolution Test . . . . . . . . . . . . . . . . . . . . . 73 4.3.2 Raw average station residuals . . . . . . . . . . . . . 75 v Table of Contents 4.3.3 P and S Upper Mantle Velocity structure beneath western Canada . . . . . . . . . . . . . . . . . . . . . 77 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.4.1 Cordillera/craton transition . . . . . . . . . . . . . . 81 4.4.2 Subsurface geometry of Northern Cascadia . . . . . . 86 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.1 Key Contributions . . . . . . . . . . . . . . . . . . . . . . . . 103 5.1.1 Improved Green’s functions for passive source struc- tural studies . . . . . . . . . . . . . . . . . . . . . . . 104 5.1.2 The teleseismic signature of fossil subduction: North- western Canada . . . . . . . . . . . . . . . . . . . . . 104 5.1.3 Body-wave tomography of western Canada . . . . . . 105 5.2 Limitations and Future work . . . . . . . . . . . . . . . . . . 106 5.2.1 Improved Green’s functions for passive source struc- tural studies . . . . . . . . . . . . . . . . . . . . . . . 106 5.2.2 The teleseismic signature of fossil subduction: North- western Canada . . . . . . . . . . . . . . . . . . . . . 106 5.2.3 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 . . . . . . . . . . . . . . . . . . . . . . 7 2.2 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 . . . . . . . . . . . . . . . . 18 2.5 Estimate of the SV-component of the Green’s function for station ULM . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Estimate of the P-component of the Green’s function for sta- tion ULM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.7 Estimate of the P and SV-component of the Green’s function for station WHY . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.8 Estimate of the P and SV-component of the Green’s function for station EDM . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Map of western Canadian geological provinces illustrating the distribution of the broad-band three-component CANOE sta- tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 re- sponse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.8 Subsurface model of the Wopmay Orogen showing the geom- etry of the dipping layer . . . . . . . . . . . . . . . . . . . . . 47 3.9 Comparison of synthetic and observed receiver functions . . . 48 3.10 Comparison of synthetic and real response . . . . . . . . . . . 50 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 braod- band 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.9 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. . . . . . . . . . . . . . . . . . . . . 89 4.14 Vertical cross sections within the P velocity model. . . . . . 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 frère Jean-Pascal, sans eux, cette thèse n’aurait sans doute pas vu le jour. Leur amour, conseils et encouragements ont été d’une valeur ines- timable tout au long de ma vie. Je ne peux passer sous silence la contribution de Karine, mon épouse, qui a partagé au cours des quatre dernières années les aleas de la vie de doctorant. Je lui serai toujours reconnaissant de amour et du support dont elle a fait preuve. xii Dedication À Karine, Claudette, Michel, Jean-Pascal et Lysanne. 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 decon- volution 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 collabora- tion with co-principal investigators on the CANOE experiment James Ga- herty, 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 col- laboration with researchers from several different institutions and published in or submitted to scientific journals. This thesis presents my major con- tributions 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 de- convolution problem is regarded by many seismologists as a key issue and an important step toward the accurate study of the structure and compo- sition of the subsurface. In the teleseismic context, deconvolution aims to unveil the earth response to a quasi-planar, impulsive wave field illuminat- ing seismic discontinuities from below. This signal is often referred to as the teleseismic Green’s function. Nowadays, in teleseismic studies, the de- convolution 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 conse- quently 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 func- tion, 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 docu- ment successfully extracts information from the P -component of the Green’s function without suffering from the problems encountered by the aforemen- tioned alternative deconvolution techniques. The second paper included in this thesis presents an important con- tribution 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 signa- ture of a Paleo-proterozoic subduction zone that had previously been docu- mented 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 pres- ence 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 re- ceiver 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 Seismo- logical 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 broad- band 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 Trans- verse 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 Hi- malaya, 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 valu- able information on crust and lithospheric structure. In the mid-1960’s, Phinney (1964) proposed that the spectral ratio of horizontal/vertical com- ponent 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) in- vestigated 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 struc- ture. Shortly thereafter, Langston (1979) introduced the time-domain “re- ceiver 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 for- ward modeled synthetic seismograms for geologically realistic models. They demonstrated that receiver functions could supply information complemen- tary to that produced using active-source seismic methods. Most early teleseismic crustal studies employing receiver functions fo- cused on characterization of the forward-scattered connected PDs phase (where a given discontinuity is designated by “D”) and considered free- surface reflected multiples (or “back scattering”) as a source of noise (see Figure 2.1 for details). The study of the PDs phase in isolation suffers from at least two important limitations. First, there is a significant trade-off be- tween 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 modu- lus 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 PpDs and PsDs). This approach has been subsequently generalized and rendered more robust through waveform stacking by Zhu and Kanamori (2000). The amplitudes of the multiples PpDs and PsDs also afford some sensitivity to density contrast but remain insensitive to compressional mod- ulus. Most recently, efforts have been directed toward exploiting receiver func- tions in the delineation of 2D and 3D Earth structure using migration/inversion methods like those routinely employed in hydrocarbon exploration (Reve- naugh, 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. Incor- poration 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 half- space 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 com- binations thereof. Information on subsurface P-impedance structure is also present within the recorded wavefield but its extraction therefrom will de- mand that we move beyond the classic “receiver function” to a more fun- damental 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 litho- spheric 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 quan- tity 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 assump- tions 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 minimum- phase 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 func- tions correlogram. In theory, as explained by Baig et al. (2005), there are methods for achieving this zero-phase/non-zero-phase component separa- tion. 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 al- ternative blind deconvolution method that we describe below. 2.2.1 Blind Deconvolution In the approach adopted here we shall consider a three-component teleseismic- P 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 ab- sorptive 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. Con- sider, for example, the cross spectrum of two components of a single three- component seismogram, P (ω)V ∗(ω) = |S(ω)|2GP (ω)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 displace- ment 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 en- ergy 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 func- tion 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 as- sume that the amplitude spectrum of the individual Green’s function com- ponents is constant in frequency and, thus, that the source power spec- trum, |S(ω)|2, is equivalent to the amplitude spectrum of the seismogram cross-correlogram, |P (ω)V ∗(ω)|. In synthetic tests comprising simple, hor- izontally layered models (Baig et al., 2005), we have gauged the accuracy of this approach by examining the quality of reproduction of the recov- ered Green’s function cross-correlogram (i.e. the inverse Fourier transform of GP (ω)GV ∗(ω)). This estimate does generally possess many of the fea- tures of the true cross-correlogram especially at negative lags which are dominated by interaction of the incident P-pulse with the first-order P- to-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 re- production 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 com- 10 Chapter 2. Improved Green’s functions for passive source structural studies −4 −3 −2 −1 0 1 2 3 4 Time [s] −4 −3 −2 −1 0 1 2 3 4 −4 −3 −2 −1 0 1 2 3 4 (a) (b) (c) 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 Gaus- sian filter (e(t 2/t2c), 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 ∗(ω))}, where F−1 denotes inverse Fourier transformation and φ(P (ω)V ∗(ω)) 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 opti- mal 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 record- ings and employ the multichannel spectral deconvolution approach origi- nally presented by Andrews (1986). We shall consider a collection of three- component 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 im- provement 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 S o u rc e S p ec tr u m E st im a te S p ec tr u m E st im a te G re en ’s F u n ct io n Short tc Optimal tc Long tc Seismogram Spectrum 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(ω)|]T the three-component Green’s function, and Sm(ω) the mth source. Transformation of the system in equation 2.2 into the log-spectral do- main 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 es- timates, S̃m, obtained by the approach of the previous section, to yield a highly overdetermined system:  I3 1 · · · 0 0 1 · · · 0 ... ... . . . ... I3 0 · · · 1 0 0 · · · 1   log{G(ω)} log{|S1(ω)|} ... log{|SM (ω)|}  =  log{P1(ω)} log{ ∣∣∣S̃1(ω)∣∣∣} ... log{PM (ω)} 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 three- component 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 trans- forms a time series to minimum-phase and exploiting the inherent minimum- phase 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(ω). Accord- ingly, we may therefore define an all-pass filter that rearranges the phase of source m to minimum phase as eiφ S m(ω) = K{Pm(ω)} Pm(ω) = K{Sm(ω)} 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 al- lows us to define two more all-pass filters, eiφ V (ω) and eiφ H(ω), that transform the non-minimum phase Green’s function components, GV (ω) and GH(ω), to their minimum-phase counterparts, that is: eiφ V (ω) = K{eiφSm(ω)Vm(ω)} eiφSm(ω)Vm(ω) = K{K{Sm(ω)}GV (ω)} K{Sm(ω)}GV (ω) = K{GV (ω)} GV (ω) , (2.7) and similarly, eiφ H(ω) = K{eiφSm(ω)Hm(ω)} eiφSm(ω)Hm(ω) = K{GH(ω)} GH(ω) . (2.8) Note that eiφ V (ω) and eiφ H(ω) must be independent of source m as they are ultimately defined solely in terms of the Green’s function. It is then straight- forward 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: GV (ω) = e−iφ V (ω)K{|GV (ω)|} and, (2.9) GH(ω) = e−iφ H(ω)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 AV (ω)eiφ V (ω) = [ 1 M M∑ m=1 e−iφ V m(ω) ]−1 and, (2.11) AH(ω)eiφ H(ω) = [ 1 M M∑ m=1 e−iφ H m(ω) ]−1 , (2.12) where eiφ V m(ω) and eiφ H 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 suc- cessfully 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 Seismo- graph Network (CNSN), namely ULM, WHY, and EDM, that clearly reveal the pure P Moho multiple PpMp. 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 layer- over-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 crust- mantle 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 prop- erties 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 , PMs and PpMs are positive whereas those of PpMp, PsMp and PsMs are negative. The amplitudes of the different first-order scattered phases, excluding PsMp, are comparable and represent approximately 5% of the incident P amplitude for this value of horizontal slowness or ray parameter. The amplitude of PsMp 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 PpMp 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 slow- ness 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 iden- tifying 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 PMs conversion and the first-order Moho multiples for the simple crustal model described in Figure 2.1. We note that the PMs 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 Slowness [s/km] Ti m e [s] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0 2 4 6 8 10 12 14 16 18 20 PpMp PMs PsMp + PpMs PsMs SV-component only P-component only P- and SV-components 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), PMs arrives between 4 and 5 s after direct P followed by PpMp near 10 s on the P component. The kinematic analogues, PsMp and PpMs follow near 14 s on separate components with the final first-order scattered phase PsMs 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. Sta- tion 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, ap- proximately 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 Eng- dahl , 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 ta- per 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 co- herent signals with timing and move-out corresponding to the three ex- pected, first-order scattered phases from the Moho (PMs, PpMs, PsMs). Figure 2.5b shows the result of stacking these time series along the PMs, PpMs and PsMs 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. PpMs and PsMs magnitudes are significantly smaller than PMs and approximately one third those predicted from syn- thetics obtained using an isotropic, non-dissipative model. This discrepancy may be due to energy loss through incoherent scattering and/or intrinsic ab- sorption, 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 PpMp 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 Ti m e [s] Slowness [s/km] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 2 4 6 8 10 12 14 16 18 20 (a) −0.1 0 0.1 −0.1 0 0.1 Time [s] 0 5 10 15 20 25 30 −0.1 0 0.1 (b) A m p li tu d e R el a ti v e to P PpMs PMs PsMs PpMs PMs PsMs 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 PMs, PpMs and PsMs. 20 Chapter 2. Improved Green’s functions for passive source structural studies Ti m e [s] Slowness [s/km] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 2 4 6 8 10 12 14 16 18 20 (a) −0.1 −0.05 0 0.05 0.1 Time [s] 5 10 15 20 25 30 −0.1 −0.05 0 0.05 0.1 (b) A m p li tu d e R el a ti v e to P P pMp PsMp? PpMp 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 PpMp and PsMp. 21 Chapter 2. Improved Green’s functions for passive source structural studies negative peak is observed. In contrast, the weaker, reflected-converted PsMp 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 PpMp 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 PMs 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 oppo- site move-out and smaller amplitude, interpreted to be the PpMs phase, is observed at 14 s showing a similar waveform. PsMs, which is expected at ∼ 18 s, is not readily apparent on the image. On the P-component image, PpMp is clearly defined at 10 s, and like PMs and PpMs, is composed of two distinct pulses of opposite polarity. This bipolar PMs 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 south- west 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 seismo- grams 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 Ti m e [s] Slowness [s/km] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 5 10 15 20 25 (a) Ti m e [s] Slowness [s/km] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 5 10 15 20 25 (b) PMs PpMp PpMs PsMs 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 dis- play the usual first-order scattered Moho phases described in preceding sec- tions, 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.5- 17 km depth. The sedimentary package (denoted by “S”) is evident as coherent energy immediately following direct P on both components. On the SV-component, PSs energy associated with the sediments dominates the Green’s function for 1.5 s, while its duration on the P-component as PpSp 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 sig- nal PpLp of opposite polarity on the P image near 4 s and PLs on the SV images at approximately 2 s. Eaton and Cassidy (1996) combined receiver function observations with Lithoprobe crustal-scale reflection data to in- terpret the layering as serpentinized former oceanic lithosphere thrust up into mid-crustal levels during an episode of Proterozoic subduction. Our observation of a coherent PpLp 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 re- sponse than possible using the conventional receiver function approach. In 24 Chapter 2. Improved Green’s functions for passive source structural studies Ti m e [s] Slowness [s/km] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 5 10 15 20 25 (a) Ti m e [s] Slowness [s/km] 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 5 10 15 20 25 (b) PpMs PsMs PpMp PMs PLs PSs PpLp PpSp 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. PSs, PLs, PpSp and PpLp). 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, un- ambiguously identified and documented the receiver-side PpMp 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 PpMp phase in many instances) are incorrectly mapped to the source and are thus ab- sent 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 Na- tional 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 PpMp 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 PMs and PpMp, is revealed to be rather more complex, representing a thin, high-velocity lid. Intra-crustal scattering dom- inates the Green’s functions recovered at station EDM. The sedimentary- basement contact produces prominent arrivals at early times followed by signals from a pronounced low-velocity layer that has been previously inter- preted from seismic reflection and receiver functions to represent a serpen- tinized 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 similar- ity of earthquakes of different sizes, 37 ed., 259-268 pp., American Geo- physical 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 normal- ization 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 two- dimensional 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 three- component 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 re- ceiver 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 Seismo- logical 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, Geophysi- cal 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 Trans- verse Ranges, California, Science, 268, 1888–1892, 1995. Sherwood, J. W. C., and A. W. Trorey, Minimum-phase and related prop- erties 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 coher- A 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 geo- physical 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) medium- aperture seismic array shows a well developed mantle stratigraphy. Three major layered structures, defined by discontinuous boundaries between az- imuthally 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 sub- duction episodes. This discovery prompted speculation on a possible con- tinuation 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 Incorpo- rated Research Institutions for Seismology (IRIS)-Lithoprobe CAnadian NOrthwest Experiment (hereafter referred to as CANOE) which comple- ments the earlier SNORCLE near-vertical reflection (Cook et al., 1998) and refraction/wide-angle reflection (Fernández-Viejo and Clowes, 2003) pro- files. 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 220˚ 220˚ 230˚ 230˚ 240˚ 240˚ 250˚ 250˚ 260˚ 260˚ 50˚ 50˚ 55˚ 55˚ 60˚ 60˚ 65˚ 65˚ 70˚ 70˚ 0 500 km Phanerozoic Proterozoic Archean Co rd ill er a W op m ay Sl av e Th el on He ar ne Tr an s− Hu ds on Edmonton Whitehorse Yellowknife Fort Nelson Leg C Leg B Leg A 237˚ 237˚ 238˚ 238˚ 239˚ 239˚ 240˚ 240˚ 241˚ 241˚ 242˚ 242˚ 243˚ 243˚ 244˚ 244˚ 60˚30' 60˚30' 61˚00' 61˚00' 61˚30' 61˚30' 62˚00' 62˚00' 62˚30' 62˚30' 0 50 km A0 8 A 09 A1 0 UB C0 1 UB C0 2 A1 1 UB C0 3 UB C0 4 A1 2 UB C0 5 UB C0 6 A1 3 UB C0 7 UB C0 8 A1 4 UB C0 9 UBC10 A15 A16 A17 Nahanni Terrane Fort Simpson Terrane Hottah Terrane Great bear magmatic arc Fort Simpson Fort Providence (a) (b) 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 Yel- lowknife (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.03Ga in age (Bowring and Williams, 1999). Its assembly was complete by 2.6Ga, 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.84Ga (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.91Ga as a magmatic arc. It collided with the Slave province during the Calderian Orogeny (1.90- 1.88Ga), shortening and displacing the sedimentary rocks of the Coronation margin. Subduction of an oceanic plate beneath the Hottah terrane between 1.88 and 1.84Ga 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.745Ga (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-2Hz band) followed by a visual inspection, resulted in a data set that comprises 1428 fair-to-good quality, broadband, three-component recordings with an aver- age 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 sub- horizontal discontinuities requires accurate removal of the earthquake source time function. This operation is conventionally performed through a spec- tral 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 se- ries 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 in- verse 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 repre- sented and laterally coherent signals, even if weak, are readily apparent. Receiver functions are ordered from left to right along the resulting pro- file based on the geographic position of the stations and, for individual sta- tions, on the sampling point of the P -S conversions projected onto the west- east 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 reflec- tion 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 dis- plays 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 (b) west east Surface Discontinuity at depth A0 8 A0 9 A1 0 UB C0 1 UB C0 2 A1 1 UB C0 3 UB C0 4 A1 2 UB C0 5 UB C0 6 A1 3 UB C0 7 UB C0 8 A1 4 UB C0 9 UB C1 0 A1 5 A1 6 A1 7 Tr an sv er se C om po ne nt Ti m e [s] Distance along the east/west direction [km] 0 50 100 150 200 250 300 350 0 5 10 15 (a) Figure 3.3: (a) Illustration of the approach adopted to display receiver func- tions 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 re- verse 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 dis- tance (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 Fernández- Viejo 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− V 2p p2 − √ 1− V 2p p2 ) , (3.1) 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 occur- rence 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 ar- ray 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 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada A0 8 A0 9 A1 0 UB C0 1 U BC 02 A1 1 UB C0 3 UB C0 4 A1 2 UB C0 5 UB C0 6 A1 3 UB C0 7 UB C0 8 A 14 UB C0 9 UB C1 0 A1 5 A1 6 A1 7 Ft Sim p. Ft Pro v. Radial Component Time [s] 0 5 10 15 0 20 40 60 80 10 0 12 0 14 0 Estimated Depth [km] Uncorrected Transverse Component Time [s] 0 5 10 15 0 20 40 60 80 10 0 12 0 14 0 Estimated Depth [km] Corrected Transverse Component Time [s] D is ta nc e al on g th e ea st /w es t d ire ct io n [km ] 0 50 10 0 15 0 20 0 25 0 30 0 0 5 10 15 0 20 40 60 80 10 0 12 0 14 0 Estimated Depth [km] (a) (b) (c) TM 1 TM 2 TM S AM AM A1 A2 F ig ur e 3. 4: R ad ia l co m po ne nt in (a ), tr an sv er se co m po ne nt in (b ), an d po la ri ty co rr ec te d (s ta ti on s A 10 th ro ug h A 15 on ly ) tr an sv er se co m po ne nt in (c ). R ed /b lu e co lo rs co rr es po nd to po si ti ve /n eg at iv e am pl it ud es ,r es pe ct iv el y. In (c ), th e co rr ec ti on is ap pl ie d to ac co un t fo r th e po la ri ty cr os s- ov er in si gn al A M do cu m en te d in se ct io n 3. 6. 3. 5. R es po ns es ar e fil te re d be tw ee n 0. 2 an d 1. 5 H z. 39 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada On the same component, a clear positive (red) pulse labeled TM (Tele- seismic 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 shal- lower, 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 am- plitude 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 conver- sion 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 grad- ual 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 ge- ometry is well defined (i.e. roughly between kilometers 150 and 220), a dominantly 1-θ (i.e. 360◦) periodicity is observed. This behaviour is illus- trated 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 Ti m e [s] Back Azimuth [°] + − − AM 50 100 150 200 250 300 350 0 2 4 6 8 10 12 14 Figure 3.5: Transverse component azimuthal response for stations UBC03, UBC04, and A12 in the frequency band 0.2-1.5Hz. Receiver functions are aligned along signal AM and stacked into 100 back-azimuth bins. Transpar- ent 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-θ am- plitude 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 rea- son 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 gen- erated 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 homo- geneous crust (Poisson solid, Vp = 6.6 km/s, ρ = 2900 kg/m3), ii) a homoge- neous 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 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada F ig ur e 3. 6: C ar to on s sh ow in g (a ) ve lo ci ty el lip so id co rr es po nd in g to m od el le d he xa go na l an is ot ro py w it h fa st sy m m et ry ax is , (b ) de fin it io n of tr en d an d pl un ge of sy m m et ry ax is , an d (c ) or ie nt at io n w it h re sp ec t to su tu re m od el . 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 ob- servations. 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 com- ponents. 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 = 7km/s to 8.1 km/s over a distance of ∼150 km. These velocities are broadly consistent with the low- velocity mantle wedge reported by Fernández-Viejo and Clowes (2003) and Fernández-Viejo et al. (2005). 45 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada Transverse Component 0 100 200 300 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 Sy nt he tic R es po ns e Radial Component 6 8 10 Sy nt he tic Ti m e [s] 6 8 10 O bs er va tio ns 0 100 200 300 6 8 10 Back Azimuth [°] Am pl itu de R el at ive to P 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 Vs ρ Anisotropy trend plunge strike dip (km/s) (km/s) (kg/m3) (±%) (◦) (◦) (◦) (◦) 8.1 4.5 3300 5 40 70 0 20 Table 3.1: Model parameters 46 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada 0 25 50 75 100 125 D e p t h [ k m ] 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 500m. Although this simplified model does not account for all the features within the observed receiver functions it provides a good geometrical tem- plate of the suture suitable for modelling the main arrivals. 47 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada R ad ia l C om po ne nt Synthetic Seismograms Time [s] 0 5 10 Real Seismograms Time [s] 10 0 15 0 20 0 25 0 0 5 10 Tr an sv er se C om po ne nt 10 0 15 0 20 0 25 0 − 0. 25 − 0. 2 − 0. 15 − 0. 1 − 0. 05 00. 05 0. 1 0. 15 0. 2 0. 25 D is ta nc e al on g th e pr of ile [k m] Amplitude Relative to P F ig ur e 3. 9: Sy nt he ti c (t op ) an d ob se rv ed (b ot to m ) re ce iv er fu nc ti on s fo r st at io ns A 10 to A 15 . L ef t pa ne ls sh ow th e ra di al re ce iv er fu nc ti on ; ri gh t pa ne ls sh ow th e tr an sv er se re ce iv er fu nc ti on . T he sa m e se t of sl ow ne ss es an d ba ck -a zi m ut hs w as us ed to ge ne ra te th e on e- w ay sy nt he ti c da ta . T he tr ac es fil te re d be tw ee n 0. 2 an d 1. 5 H z. Su rf ac e se di m en ts re ve rb er at io ns se en on th e re al ra di al co m po ne nt ar e no t m od el ed in th e si m ul at io n. N ot e th at th e sy nt he ti cs re pr od uc e th e m ai n fe at ur es of se en in th e re al ra di al an d tr an sv er se co m po ne nt . 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 syn- thetic 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 pos- itive 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 seis- mic 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 (Fernández Viejo et al., 1999; Fernández-Viejo and Clowes, 2003), and teleseismic re- ceiver 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 fil- tered, reflection section (dominant frequency >15Hz) upon the radial- and transverse-component receiver functions (dominant frequency < 1Hz) in Figure 3.11. This is accomplished through a simple scaling between the di- rect P -to-S conversion time, tps, and the 2-way P travel-time, t2p in the form t2p = 2tps √ 1− V 2p p2(√ 3− V 2p p2 − √ 1− V 2p p2 ) , (3.2) 49 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada Sy nt he tic Ti m e [s] 4 6 8 10 12 CA NO E Ti m e [s] Back Azimuth [°] 50 100 150 200 250 300 350 4 6 8 10 12 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 = 6km/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 Chapter 3. The teleseismic signature of fossil subduction: Northwestern Canada A0 8 A0 9 A1 0 UB C0 1 U BC 02 A1 1 UB C0 3 UB C0 4 A1 2 UB C0 5 U BC 06 A1 3 UB C0 7 UB C0 8 A1 4 UB C0 9 UB C1 0 A 15 A1 6 A1 7 Ft. Si mp . Ft. Pr ov . Radial Component Two Way P−Travel Time [s] 0 5 10 15 20 25 30 35 40 Transverse Component Two Way P−Travel Time [s] D is ta nc e Al on g th e Pr of ile [k m] 0 50 10 0 15 0 20 0 25 0 30 0 35 0 40 0 45 0 0 5 10 15 20 25 30 35 40 a ) b) 0 030 3070 7011 0 11 0 15 0 15 0 Estimated Depth [km] Estimated Depth [km] TM 1 TM 2 TM TM AM A1 A2 F ig ur e 3. 11 : (a ) R ad ia l an d (b ) po la ri ty -c or re ct ed tr an sv er se re ce iv er fu nc ti on s up on w hi ch ti m e m ig ra te d, co he re nc y- fil te re d, ne ar -v er ti ca lr efl ec ti on da ta ar e su pe ri m po se d w it h th e po si ti on of th e st at io ns al on g th e pr ofi le sh ow n on to p. T he di sp la y fo rm at is sl ig ht ly di ffe re nt fr om th e on e ad op te d in F ig ur e 3. 4 in th e se ns e th at re ce iv er fu nc ti on ar e pr oj ec te d al on g th e pr ofi le in st ea d of th e w es t- ea st di re ct io n. T o ag re e w it h tw o- w ay P tr av el ti m e of th e re fle ct io n pr ofi le a si m pl e ti m e sc al in g, de sc ri be d in (3 .2 ) w as ap pl ie d. A n es ti m at e of th e de pt h ba se d on a si m pl is ti c tw o la ye rs fe at ur in g a cr us t an d m an tl e w it h ve lo ci ty of 6 an d 8 km /s ,r es pe ct iv el y is pr ov id ed on th e ri gh t ax is . T he te le se is m ic re sp on se s ar e fil te re d be tw ee n 0. 2 an d 1. 5 H z. T he po si ti on of th e st at io ns al on g th e pr ofi le is sh ow n on to p. 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 re- flection/refraction data indeed show high crustal and low mantle velocities along this stretch of profile (Fernández-Viejo and Clowes, 2003). Further study by Fernández-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 equiv- alently 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 sub- duction event, but the lithologies responsible are difficult to constrain. It is interesting to note that the teleseismic results present evidence for struc- tural 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 sim- ple 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 Yel- lowknife: 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 sub- ducted plate observed to the west. Moreover, he interpreted the layering, and in particular H, as evidence for stranded oceanic crust that had de- veloped 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 Ye llo wk nif e Ti m e [s] 4 6 8 10 12 CA NO E Ti m e [s] Slowness [s/km] 0.045 0.05 0.055 0.06 0.065 0.07 0.075 4 6 8 10 12 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.5Hz 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 re- ceiver 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 impli- cations 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 lo- calized (10-15 km thick) zone of mantle material, AM, may acquire its anisotropy, under the assumption that it represents lattice preferred ori- entation 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 subduc- tion 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 temper- ature deformation of olivine can take place within 2-10 km thick shear do- mains that have been ascribed to transform faults (Prinzhofer and Nicolas, 1980). It would appear reasonable, therefore, that low-temperature, local- ized deformation might also occur in subduction environments, especially in the near vicinity of the plate interface. The location of AM immedi- ately 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 overly- ing 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 under- lying mantle and produces compressive stress (and consequent deformation) over a depth interval below the crust-mantle boundary that is roughly com- mensurate 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 litholo- gies 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 corre- sponding 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 de- velopment 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 olivine- rich harzburgites predominate within the oceanic lithosphere from the Moho to depths of ∼ 25km (Nicolas et al., 1980). These rocks are plastically de- formed at the high temperatures prevailing near the ridge, with particularly large strains developing close to the Moho. Orientations of olivine and py- roxene 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 orthopyrox- ene within peridotite increases. These observations present an alternative interpretation for AM, that is, that it represents a gradient zone of increas- ing 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 spread- ing 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 subduc- tion 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 re- semblance between the structure AM and H. If AM and H represent the same subduction episode and form part of one contiguous structure, the im- plication is that the bulk of mantle lithosphere beneath the western Slave province is Paleoproterozoic (ca. 1.8Ga) 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.6Ga (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 unfor- tunate 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 litho- sphere 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 dif- 58 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 fine- scale (∼10 km thick) anisotropic layering lends further credibility to the thesis that Archean cratons were stabilized through shallow subduction pro- cesses (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 ex- hibiting 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 there- fore signifies cratonic-scale lithospheric underthrusting. There are at least 3 published examples of comparable structures de- fined 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 disconti- nuity 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 Kaap- vaal 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 vari- ations 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 consis- tent 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 Territo- ries provide unambiguous evidence for the signature of fossil subduction in the form of highly anisotropic, shallow-most mantle characterizing the un- derthrust 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 as- sistance in coordination of the field effort. Discussions with Jounada Oueity and Ron Clowes helped to focus our initial analysis. We are grateful to As- sociate 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 subduc- tion, Geophys. Res. Lett., 18, 585–588, 1991. Audet, P., M. G. Bostock, and J.-P. Mercier, Teleseismic waveform mod- elling 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, Dor- drecht, 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 litho- probe Slave-northern cordillera lithospheric evolution (SNORCLE) tran- sect., 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 litho- spheric reflection profiling of the western canadian shield, Tectonics, 18, 1–24, 1999. Davis, W. J., A. G. Jones, W. Bleeker, and H. Grutter, Lithosphere de- velopment in the Slave craton: a linked crustal and mantle perspective, Lithos, 71, 575–589, 2003. Fernández-Viejo, G., and R. M. Clowes, Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay Orogen, nortwest- ern Canada, from a lithoprobe refraction/wide-angle reflection survey, Geophys. J. Int., 153, 1–19, 2003. Fernández Viejo, G., R. Clowes, and J. Amor, Imaging the lithospheric mantle in northwestern Canada with seismic wide-angle reflections, Geo- phys. Res. Lett., 26, 2809–2812, 1999. Fernández-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 discontinu- ities 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 evolu- 63 Bibliography tion, in it Kimberlites and Related Rocks, 358-368 pp., Blackwell, Cam- bridge, 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 three- component 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 upper- most 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 transi- tion from convergent to transform boundary. In this paper we undertake a continental-scale study of the western Canada upper mantle using teleseis- mic 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, surface- wave, 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 body- waves (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 Septem- ber 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 Frederik- sen 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 Cana- dian lithospheric and upper-mantle structure but are either global in cover- age 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 re- gion as performed for example, in the western United States (Humphreys and Dueker , 1994). However, the recent deployment of stations from sev- eral portable experiments (i.e. CANOE, BATHOLITHS, POLARIS BC and POLARIS-Nechako) have contributed to filling some major gaps in sta- tion 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 Data and Method 4.2.1 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 144 oW 136oW 128oW 120 oW 112 oW 45 oN 50 oN 55 oN 60 oN 65 oN BATHOLITHS CANOE CNSN POLARIS OTHERS Figure 4.1: Map of western Canada illustrating the distribution of braod- band and short-period seismic stations used in this study (ANSS), Alaska Tsunami Warning System (ATWS), Alaska Regional Net- work (ARN), BATHOLITHS, CAnadian NOrthwest Experiment (CANOE), Canadian National Seismograph Network (CNSN), and several experiments of the Portable Observatories for Lithospheric Analysis and Research Inves- tigating Seismicity (POLARIS). The P -wave data-set consists of 23, 420 tele- seismic 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 sub- sets based on station location. For each subset, P and S traveltimes were obtained from visual picks of the vertical and transverse component, re- spectively, which were subsequently refined through multichannel cross- correlation (VanDecar and Crosson, 1990). Frequency bands of 0.4Hz to 1.6Hz and 0.05Hz to 0.4Hz 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 slow- ness 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 x Φd + λΦm (4.1) where λ is a free parameter that controls the level of smoothness and flatness of the model, Φd def= ‖ Ax− b ‖2, Φm def= ‖ ∆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 ] , (4.2) x =  ∆sc 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 signal- to-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 Results 4.3.1 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 event- station ray geometry of the actual data. Gaussian noise with standard deviations of 30ms and 200ms was subsequently added to the P and S data- sets, 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 de- pends, 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 mod- els. 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 sys- tematically underestimated and their magnitudes are typically reduced by about 30% compared to the input anomalies. The underestimated ampli- tudes 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 perturba- tion 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 low- velocity structures whereas negative residuals manifest high-velocity anoma- lies. Figures 4.7a, b combine the average station residuals for the four az- imuthal 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 the P (top) and S (bottom) data-sets. The blue and red isosurfaces repre- sent velocity perturbations of the IASP91 reference model of +1% and -1%, respectively. 76 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 suggest- ing fine-scale lateral heterogeneity. Positive residuals are generally observed west of the Cordilleran deformation front whereas negative residuals are ob- served east of that boundary. Southwestern British Columbia, where most of the mainland stations east of Vancouver Island display negative residu- als, 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 S- velocity 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 144 oW 136oW 128oW 120 oW 11 2oW 45 oN 50 oN 55 oN 60 oN 65 oN −1 sec +1 sec 144 oW 136oW 128oW 120 oW 11 2oW 45 oN 50 oN 55 oN 60 oN 65 oN −2.5 sec +2.5 sec 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 Chapter 4. Body-wave tomography of western Canada 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 10 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 20 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 30 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 40 0 km − 2 − 1012 F ig ur e 4. 8: D ep th sl ic es co ve ri ng th e w ho le su rv ey ar ea at 10 0 km (t op le ft ), 20 0 km (t op ri gh t) ,3 00 km (b ot to m le ft ), an d 40 0 km (b ot to m ri gh t) th ro ug h th e P ve lo ci ty m od el . T hi s m od el w as ob ta in ed fr om 10 00 co nj ug at e gr ad ie nt it er at io ns w it hi n 15 H ub er do w nw ei gh ti ng s us in g λ = 35 0 an d µ = 60 . B lu e an d re d co lo rs re pr es en t po si ti ve an d ne ga ti ve ve lo ci ty pe rt ur ba ti on s, re sp ec ti ve ly . 79 Chapter 4. Body-wave tomography of western Canada 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 10 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 20 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 30 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 40 0 km − 2 − 1012 F ig ur e 4. 9: D ep th sl ic es co ve ri ng th e w ho le su rv ey ar ea at 10 0 km (t op le ft ), 20 0 km (t op ri gh t) ,3 00 km (b ot to m le ft ), an d 40 0 km (b ot to m ri gh t) th ro ug h th e S ve lo ci ty m od el . T hi s m od el w as ob ta in ed fr om 10 00 co nj ug at e gr ad ie nt it er at io n w it hi ns 40 H ub er do w nw ei gh ti ng s us in g a λ = 35 0 an d µ = 60 . B lu e an d re d co lo rs re pr es en t po si ti ve an d ne ga ti ve ve lo ci ty pe rt ur ba ti on s, re sp ec ti ve ly . 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 Discussion 4.4.1 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 tran- sition from Phanerozoic mantle to cratonic lithosphere is often characterized by a marked change in physical properties that is typically manifest in to- mographic 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); Ekström et al. (1997); Grand et al. (1997); Si- mons 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 Fred- eriksen, 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 evi- dence 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 (Hynd- man and Lewis, 1999; Lewis et al., 2003), as well as elastic thickness studies (Flück 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 ori- gin 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 ar- ray afford an excellent lateral resolution of structure across the northern Cordillera. They are consequently well positioned to locate and character- ize 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 su- tures 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 Cordilleran frontTintina Fault (A)(B)(C)0 100 200 300 400 500 600 700 De pt h [km ] A1 (60.5N, 138.7W) A2 (60.5N, 108.7W) 144 oW 136oW 128oW 120oW 112 oW 54 oN 56 oN 58 oN 60 oN 62 oN 64 oN A1 A2 % of P velocity −3 −1.5 0 1.5 3 Figure 4.10: Cross section within the P velocity model along the corridor A1- A2 (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 deforma- tion 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 consis- tently 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 temper- ature documented by Hyndman and Lewis (1999) and Lewis et al. (2003) from heat flow measurements, and to a lesser extent by Flück et al. (2003) and Audet et al. (2007) from the study of elastic thickness. This combi- nation 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 north- western Canada with some previous studies can be attributed to the pres- ence 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 ei- ther 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 −1000 −800 −600 −400 −200 0 200 400 600 800 1000 −8 −6 −4 −2 0 2 4 6 8 10 Distance from the cordilleran deformation front [km] % o f p er tu rb at io n Ti nt in a Fa ul t P velocity S velocity Frederiksen 2001 Grand 1997 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 particu- lar. The geometry of its subsurface is only loosely constrained by scattered heat flow, gravity and geochemical data (Lewis et al., 1997). The recent de- ployment of broadband three-component seismic stations on northern Van- couver 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, high- velocity 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 con- straints 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 Au- det et al. (2008) focused on the structure of the subducted Explorer plate 86 Chapter 4. Body-wave tomography of western Canada 132oW 128oW 124oW 120 oW 1 16oW 48oN 50oN 52oN 54oN 56oN 58oN Na zk o Itc haIlg ac hu z Ra inb ow Mi lba nk e P velocity perturbation [%] (E) (D) (F) (G) 2 0 2 EXP JdF PAC WB 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 re- sults, these authors proposed a tectonic model in which the separation of the Explorer plate from the Juan de Fuca plate is described as the conse- quence 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. Bo- stock 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 Chapter 4. Body-wave tomography of western Canada 13 2o W 12 8o W 12 4o W 12 0o W 11 6o W 48 o N 50 o N 52 o N 54 o N 56 o N 58 o N B1 B2 A1 A2 C1 C2 D1 D2 % o f P ve lo ci ty (E ) (D ) (F ) (G ) Lo ca tio n of c ro ss −s ec tio n − 2 0 2 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] A1 (48 N, 12 4W ) A2 (51 .8N , 1 17 W ) (E) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] B1 (49 N, 12 6W ) B2 (52 .3N , 1 20 W ) (E) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] C1 (51 .5N , 1 29 W ) C2 (54 N, 12 4.5 W ) (D) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] D 1 (52 .3N , 1 32 W ) D 2 (55 N, 12 7W ) (G ) F ig ur e 4. 13 : V er ti ca l cr os s se ct io ns w it hi n th e P ve lo ci ty m od el in th e no rt he rn C as ca di a re gi on . 89 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 continu- ous. 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 Casca- dia 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 prin- cipally 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 de- formation 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 west- east 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 3- 1 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 magma- tism in this region. Rohr and Currie (1997) suggested that simple lithospheric- scale 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. As- thenospheric 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 po- tential 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 Yel- lowstone (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 consis- tent with the plume origin (Bevier , 1989). Figure 4.14 presents cross sections perpendicular and parallel to the in- ferred 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 lo- cations 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 136oW 132oW 128oW 124oW 120 oW 1 16oW 48oN 50oN 52oN 54oN 56oN Na zko Itc haIlga chu z Ra inb ow Mi lba nk e E1 E2 F1 F2 % of P velocity (E) (D) (F) (G) Location of cross−section −2 0 2 0 100 200 300 400 500 600 700 D ep th [k m] E1 (52.3N, 131.1W) E2 (53.5N, 116.4W) Nazko Cone (D) (F) 0 100 200 300 400 500 600 700 D ep th [k m] F1 (48.4N, 123.4W) F2 (57.4N, 124.1W) Nazko Cone (E) (F) 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 low- velocity anomaly present beneath Nazko Cone that extends seaward along the Anahim volcanic track suggests a direct relation to the magmatism ob- served 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 be- tween 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 be- neath a continent (e.g. (Dickinson and Snyder , 1979; Johnson and O’Neil , 1984; Forsythe and Nelson, 1985; Thorkelson and Taylor , 1989; Thorkel- son, 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 sug- gestive 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 extend- ing 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 gen- erally 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 character- ized 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 signifi- cant 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 pro- cess 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 deforma- tion 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 Ex- plorer/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 petroge- nesis 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 man- tle 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 north- ern 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 under- thrusting 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 Amer- ica 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 Litho- sphere 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. Ekström, 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. Flück, 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 (1- 2), 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 (3- 4), 175–191, 2001. Grand, S., V. der Hilst R.D., andW. 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 (3- 4), 397–422, doi:10.1016/S0040-1951(99)00068-2, 1999. Hyndman, R., C. Currie, and S. Mazzotti, Subduction zone backarcs, mo- bile 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 electromag- netic 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 ridge- trench-triple junction: Northern Vancouver Island, J. Geophys. Res., 102 (B4), 7767–7781, 1997. Lewis, T., R. Hyndman, and P. Flück, 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 Moun- tains 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 to- mographic studies, Geology, 23 (3), 205–208, 1995. Price, R., and D. Carmichael, Geometric test for Late Cretaceous- Paleogene 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 Let- ters, 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 tele- seismic 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 vol- cano 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 deconvolu- tion 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 deconvo- lution technique that provides a more accurate estimate of the so-called tele- seismic 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 mean- ingful 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 upper- mantle 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 tech- nique 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 PpMp multiple. Station ULM exhib- ited 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 por- trayed through PMs and PpMp, 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 be- tween 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 tran- sect. Teleseismic receiver functions computed from ∼250 earthquakes clearly reveal the response of the ancient subduction zone. On the radial compo- nent, 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 up- 104 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 char- acterized 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 struc- ture 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 sta- tions, 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 veloc- ity 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 CA- NOE (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 northwest- ern 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 northwest- ern 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 cur- rently 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 consis- tently missing from estimates (e.g. PsMp phase see Figure 2.1). To mitigate this problem, a better source estimation scheme has to be developed and im- plemented in the existing framework. The recent breakthroughs in the field of large polynomial factorization and in the development of robust optimiza- tion 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 docu- mented by Bostock (1998) and that of AM bear a close similarity and sug- 106 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 identi- fied AM with paleosubduction, the limited observation of similar structures precludes us from finding a conclusive mechanism for its formation. More re- search 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 lim- itation we encountered throughout this project. The station density varied widely between the different region. The station density was high along CA- NOE, 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 re- liable 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 consis- tently 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 Appendix A. Tomography Resolution Test 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 10 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 20 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 30 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 40 0 km − 2 − 1012 F ig ur e A .1 : In pu t P re so lu ti on te st m od el fo r th e de pt h sl ic es pr es en te d in F ig ur e 4. 8. 110 Appendix A. Tomography Resolution Test 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 10 0 km − 0. 5 00. 5 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 20 0 km − 1 − 0. 5 00. 5 1 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 30 0 km − 0. 5 00. 5 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of P velocity 40 0 km − 1 − 0. 5 00. 5 1 F ig ur e A .2 : O ut pu t of th e re so lu ti on te st fo r th e P da ta -s et fo r th e de pt h sl ic es pr es en te d in F ig ur e 4. 8. 111 Appendix A. Tomography Resolution Test 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 10 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 20 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 30 0 km − 2 − 1012 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 40 0 km − 2 − 1012 F ig ur e A .3 : In pu t S re so lu ti on te st m od el fo r th e de pt h sl ic es pr es en te d in F ig ur e 4. 9. 112 Appendix A. Tomography Resolution Test 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 10 0 km − 1 − 0. 5 00. 5 1 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 20 0 km − 1 − 0. 5 00. 5 1 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 30 0 km − 1 − 0. 5 00. 5 1 14 4o W 13 6o W 12 8o W 12 0o W 11 2o W 45 o N 50 o N 55 o N 60 o N 65 o N Tin tina Fa ult Cor dille ran defo rma tion fron t % of S velocity 40 0 km − 1 − 0. 5 00. 5 1 F ig ur e A .4 : O ut pu t of th e re so lu ti on te st fo r th e S da ta -s et fo r th e de pt h sl ic es pr es en te d in F ig ur e 4. 9. 113 Appendix A. Tomography Resolution Test Cordilleran frontTintina Fault 0 100 200 300 400 500 600 700 De pt h [km ] A1 (60.5N, 138.7W) A2 (60.5N, 108.7W) 144 oW 136oW 128oW 120oW 112 oW 54 oN 56 oN 58 oN 60 oN 62 oN 64 oN A1 A2 % of P velocity −3 −1.5 0 1.5 3 Figure A.5: Input P resolution test model for the cross sections presented in Figure 4.10. Cordilleran frontTintina Fault 0 100 200 300 400 500 600 700 De pt h [km ] A1 (60.5N, 138.7W) A2 (60.5N, 108.7W) 144 oW 136oW 128oW 120oW 112 oW 54 oN 56 oN 58 oN 60 oN 62 oN 64 oN A1 A2 % of P velocity −3 −1.5 0 1.5 3 Figure A.6: Output of the resolution test for the S data-set for the depth slices presented in Figure 4.10. 114 Appendix A. Tomography Resolution Test 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] A1 (48 N, 12 4W ) A2 (51 .8N , 1 17 W ) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] B1 (49 N, 12 6W ) B2 (52 .3N , 1 20 W ) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] C1 (51 .5N , 1 29 W ) C2 (54 N, 12 4.5 W ) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] D 1 (52 .3N , 1 32 W ) D 2 (55 N, 12 7W ) F ig ur e A .7 : In pu t P re so lu ti on te st m od el fo r th e cr os s se ct io ns pr es en te d in F ig ur e 4. 13 . 115 Appendix A. Tomography Resolution Test 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] A1 (48 N, 12 4W ) A2 (51 .8N , 1 17 W ) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] B1 (49 N, 12 6W ) B2 (52 .3N , 1 20 W ) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] C1 (51 .5N , 1 29 W ) C2 (54 N, 12 4.5 W ) 0 10 0 20 0 30 0 40 0 50 0 60 0 Depth [km] D 1 (52 .3N , 1 32 W ) D 2 (55 N, 12 7W ) F ig ur e A .8 : O ut pu t of th e re so lu ti on te st fo r th e P da ta -s et fo r th e de pt h sl ic es pr es en te d in F ig ur e 4. 13 . 116 Appendix A. Tomography Resolution Test 136oW 132oW 128oW 124oW 120 oW 1 16oW 48oN 50oN 52oN 54oN 56oN Na zko Itc haIlga chu z Ra inb ow Mi lba nk e E1 E2 F1 F2 % of P velocity Location of cross−section −2 0 2 0 100 200 300 400 500 600 700 D ep th [k m] E1 (52.3N, 131.1W) E2 (53.5N, 116.4W) Nazko Cone 0 100 200 300 400 500 600 700 D ep th [k m] F1 (48.4N, 123.4W) F2 (57.4N, 124.1W) Nazko Cone Figure A.9: Input P resolution test model for the cross sections presented in Figure 4.14. 117 Appendix A. Tomography Resolution Test 136oW 132oW 128oW 124oW 120 oW 1 16oW 48oN 50oN 52oN 54oN 56oN Na zko Itc haIlga chu z Ra inb ow Mi lba nk e E1 E2 F1 F2 % of P velocity Location of cross−section −2 0 2 0 100 200 300 400 500 600 700 D ep th [k m] E1 (52.3N, 131.1W) E2 (53.5N, 116.4W) Nazko Cone 0 100 200 300 400 500 600 700 D ep th [k m] F1 (48.4N, 123.4W) F2 (57.4N, 124.1W) Nazko Cone Figure A.10: Output of the resolution test for the P data-set for the depth slices presented in Figure 4.14. 118
Thesis/Dissertation
eng
Geophysics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Improved teleseismic Green's functions and western Canada mantle structure and evolution
Text