Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Near-vertical and wide-angle seismic reflection studies of the Moho and sub-crustal lithosphere in NW… Oueity, Jounada 2010

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

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

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

Full Text

Near-vertical and wide-angle seismic reflection studies of the Moho and sub-crustal lithosphere in NW Canada by Jounada Oueity  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Geophysics)  The University of British Columbia (Vancouver) September, 2010  © Jounada Oueity 2010  Abstract High quality, coincident near-vertical incidence (NVI) and refraction/wide-angle seismic reflection data (R/WAR) acquired along a profile in the Northwest Territories are used to study the nature of the Moho and subcrustal reflectors. First, we re-examine distinct subhorizontal reflections on NVI data in the uppermost mantle that were interpreted previously as a separate feature from a relict subducted slab. Using forward and inverse traveltime modeling of both data sets along the crooked line, we investigate the origin of the reflections. Our results demonstrate that the subhorizontal reflectors are the continuation of the relict subducted slab, which extends laterally for 300 km at depths from 35 to 90 km. Its base is the source of the R/WAR reflections. The apparent flattening is an artifact of projecting a 3-D geometry onto a 2-D cross section. The shallowly subducted slab probably contributed to the thickening and stabilization of the sub-crustal lithosphere in the region. Second, we examine the detailed structure of the Moho and propose a possible scenario for its formation and evolution. Strong Moho reflections are observed on the NVI data (shot gathers and stacked section). The WA data are characterized by a ~0.5 s coda trailing the PmP (Moho) phase. For analysis of these observed data, we follow two approaches, forward and inverse modeling. In forward modeling, we calculate wide-angle and near-vertical synthetic seismograms using 1- and 2-D wave propagation algorithms. Comparison between synthetic and observed data for shot gathers was made possible through development of a novel noise-removal technique using the curvelet transform. For the inverse method, we use a statistical analysis approach based on the von Karman autocorrelation function. Our results indicate that the Moho is a finite-thickness (~3 km), heterogeneous transition zone. The heterogeneities can be described by laterally discontinuous layering, lamellae structure with randomly distributed ellipses or a von Karman distribution with a lateral correlation length of 936 m. The transition zone separates the lower crust with a lateral correlation length of 732 m from the uppermost mantle with a correlation length of 261 m. The Moho is interpreted as a thermal/metamorphic front, a regional décollement, or both.  ii  Table of Contents Abstract .................................................................................................................................................. ii Table of Contents .................................................................................................................................. iii List of Tables ........................................................................................................................................ vii List of Figures ..................................................................................................................................... viii Acknowledgements ................................................................................................................................ x Dedication ............................................................................................................................................. xi Statement of Co-authorship .................................................................................................................. xii  Chapter 1 Introduction, Objectives and Background .................................................................... 1 1.1 Introduction .................................................................................................................................. 1 1.2 Theme and research objectives ..................................................................................................... 3 1.3 Background .................................................................................................................................. 4 1.4 Thesis outline ............................................................................................................................... 7 Bibliography ......................................................................................................................................... 10  Chapter 2 Paleoproterozoic subduction in NW Canada from near-vertical and wide-angle reflection data ..................................................................................................................................... 14 2.1 Introduction ................................................................................................................................ 14 2.2 Tectonic setting .......................................................................................................................... 18 2.3 Seismic experiments and data .................................................................................................... 20 2.4 Analysis procedures and results ................................................................................................. 22 2.4.1 NVI reflection data .............................................................................................................. 22 2.4.2 Two-dimensional ray-tracing .............................................................................................. 25 2.4.3 Two-dimensional models .................................................................................................... 27 2.4.4 Three-dimensional forward algorithm ................................................................................. 28 2.4.5 Three-dimensional models from forward algorithm ........................................................... 28 2.4.6 Three-dimensional finite-difference traveltime inversion ................................................... 31 2.4.7 Three-dimensional models from traveltime inversion......................................................... 32 2.5 Discussion .................................................................................................................................. 35 2.5.1 Slab geometry ...................................................................................................................... 35  iii  2.5.2 Origin of the seismic signature of the slab .......................................................................... 38 2.5.3 Flat versus low-angle subduction ........................................................................................ 40 2.5.4 Relation with other upper-mantle reflectors ........................................................................ 43 2.6 Conclusions ................................................................................................................................ 43 Bibliography ......................................................................................................................................... 45  Chapter 3 Enhancing crustal reflection data through curvelet denoising ................................. 50 3.1 Introduction ................................................................................................................................ 50 3.2 Curvelet denoising ...................................................................................................................... 53 3.3 Application to synthetic data ...................................................................................................... 56 3.3.1 Application to a synthetic shot gather ................................................................................. 58 3.3.2 Application to synthetic post-stack data .............................................................................. 60 3.4 Application to crustal reflection data.......................................................................................... 62 3.4.1 Application to a vibroseis shot gather ................................................................................. 65 3.4.2 Application to a stacked seismic section ............................................................................. 67 3.5 Discussion .................................................................................................................................. 70 3.6 Conclusions ................................................................................................................................ 71 Bibliography ......................................................................................................................................... 73  Chapter 4 The nature of the Moho transition in NW Canada from combined near-vertical and wide-angle seismic reflection studies ...................................................................................... 76 4.1 Introduction ................................................................................................................................ 76 4.2 Tectonic background .................................................................................................................. 79 4.3 Seismic experiment and data ...................................................................................................... 81 4.3.1 SNORCLE Line 1 ............................................................................................................... 81 4.3.2 Processing of the reflection data.......................................................................................... 82 4.3.3 Processing of the refraction data ......................................................................................... 84 4.4 Synthetic seismic modeling ........................................................................................................ 88 4.4.1 Methodology ....................................................................................................................... 88 4.4.2 Synthetic models ................................................................................................................. 90 4.5 Discussion ................................................................................................................................ 104 4.5.1 Synthetic models ............................................................................................................... 104 iv  4.5.2 Multigenetic origin for the continental Moho ................................................................... 105 4.6 Conclusions .............................................................................................................................. 109 Bibliography ....................................................................................................................................... 112  Chapter 5 Heterogeneity of the lower crust, Moho and upper mantle from statistical anaylsis of near-vertical and wide-angle seismic reflection data ............................................................. 121 5.1 Introduction .............................................................................................................................. 121 5.2 The Wopmay orogen ................................................................................................................ 124 5.3 Seismic data .............................................................................................................................. 126 5.4 Estimation of power law parameters ........................................................................................ 128 5.5 Two-dimensional finite-difference modeling ........................................................................... 129 5.6 Calibration approach ................................................................................................................ 130 5.6.1 Estimation of von Karman parameters from SNORCLE Line 1 NVI data ....................... 131 5.6.2 Synthesis of velocity models ............................................................................................. 135 5.6.3 Simulation of wide-angle and near-vertical data ............................................................... 138 5.7 Attributes analysis .................................................................................................................... 141 5.7.1 Complex trace envelope .................................................................................................... 141 5.7.2 Envelope energy decay ...................................................................................................... 143 5.7.3 F-X analysis ....................................................................................................................... 143 5.8 Discussion ................................................................................................................................ 146 5.8.1 Fabric in the crust and uppermost mantle .......................................................................... 146 5.8.2 Constraining the correlation length and the upscaling factor ............................................ 149 5.8.3 Character and attributes of wide angle reflection data ...................................................... 151 5.8.4 The Moho transition zone .................................................................................................. 152 5.9 Conclusions .............................................................................................................................. 154 Bibliography ....................................................................................................................................... 156  Chapter 6 Conclusions................................................................................................................... 162 6.1 Main contributions.................................................................................................................... 162 6.1.1 Shallow relict subduction zone .......................................................................................... 162 6.1.2 Application of curvelet denoising ..................................................................................... 163 6.1.3 The nature of the Moho transition ..................................................................................... 164 v  6.1.4 The stochastic lower crust, the Moho and the uppermost mantle...................................... 164 6.1.5 The origin and tectonic evolution of the Moho ................................................................. 165 6.2 Limitations................................................................................................................................ 165 6.3 Future work .............................................................................................................................. 168 6.3.1 Application of curvelet denoising ..................................................................................... 168 6.3.2 Application of quantitative analysis .................................................................................. 168 6.3.3 More and better data leading to 3-D surveys ..................................................................... 169 Bibliography ....................................................................................................................................... 171  Appendix A Wave propagation algorithms ................................................................................. 173 A.1 One dimensional reflectivity method....................................................................................... 173 A.2 Two dimensional visco-elastic finite difference method ......................................................... 174 Bibliography ....................................................................................................................................... 176  Appendix B Supplementary information for Chapter 5............................................................. 177 B.1 Estimation of power law parameters........................................................................................ 177 Bibliography ....................................................................................................................................... 179  vi  List of Tables Table 2.1. Parameters for reprocessing of NVI data. .................................................................... 25 Table 2.2. Statistics of model fitting. ............................................................................................ 27 Table 2.3. Statistics of model fitting. ............................................................................................ 34 Table 3.1. Noise parameters for synthetic models......................................................................... 60 Table 3.2. Data acquisition parameters. ........................................................................................ 65 Table 4.1. Processing parameters for NVI data. ............................................................................ 83 Table 4.2. Average uncertainties. .................................................................................................. 88 Table 4.3. Different combination of all the parameters ............................................................... 101 Table 5.1. Processing parameters for NVI data. .......................................................................... 126 Table 5.2. Velocity models stochastic parameters. ..................................................................... 138  vii  List of Figures Figure 2.1. Tectonic map of the study area. .................................................................................. 15 Figure 2.2 The originally migrated seismic data. .......................................................................... 16 Figure 2.3 Record section for shot 1101.. ..................................................................................... 16 Figure 2.4. a) Location of leg A, part of CANOE experiment. ..................................................... 17 Figure 2.5. Comparison between original and re-processed section of NVI data. ........................ 24 Figure 2.6. Velocity model and ray path diagram. ........................................................................ 26 Figure 2.7. Comparison between observed and synthetic data from 3-D forward modeling.. ...... 30 Figure 2.8. Results for 3-D finite difference inversion for the WAR and NVI data.. ................... 33 Figure 2.9. Comparison of observed and calculated traveltimes for WAR and NVI data. ........... 35 Figure 2.10. 3-D perspective view of the subducted slab. ............................................................. 36 Figure 2.11. The re-interpreted seismic reflection profile ............................................................. 37 Figure 2.12. Schematic illustration of the subducted slab. ............................................................ 38 Figure 3.1. Curvelets in time-space and frequency-wavenumber domains. .................................. 52 Figure 3.2. Discrete curvelet partitioning of the 2-D Fourier plane. ............................................. 55 Figure 3.3. The cooling method with iterative thresholding. ........................................................ 55 Figure 3.4. Signal-to-noise ratio vs. threshold values for white noise. ......................................... 57 Figure 3.5. Noise-free shot gather computed for a crustal-scale velocity model . ........................ 59 Figure 3.6. Noise-free shot gather computed for a synthetic salt dome model. ............................ 61 Figure 3.7. Synthetic post-stack data from the Marmousi model. ................................................. 63 Figure 3.8. Tectonic map of northwestern Canada. ....................................................................... 64 Figure 3.9. Shot gather from SNORCLE Line 1 segment. ............................................................ 66 Figure 3.10. A 60-km-long segment of stacked curvelet denoised seismic section ...................... 69 Figure 4.1. Tectonic map of the study area. .................................................................................. 80 Figure 4.2. The seismic reflection profile...................................................................................... 84 Figure 4.3. Shot gathers after curvelet denoising. ......................................................................... 85 Figure 4.4. Record section for shots 1104, 1105 and 1108.. ......................................................... 86 Figure 4.5. Amplitude and amplitude spectra of the observed PmP phase. .................................. 87 Figure 4.6. Synthetic and observed seismograms for a layered transition. ................................... 92 Figure 4.7. Synthetic seismograms of a layered velocity model. .................................................. 93 viii  Figure 4.8. Comparison between synthetic and observed amplitude of the PmP phase ............... 94 Figure 4.9. Normalized amplitude spectra of synthetic PmP phase .............................................. 94 Figure 4.10. Synthetic near-vertical data for a layered transition model...................................... 95 Figure 4.11. Normalized amplitude spectra of synthetic near-vertical data .................................. 95 Figure 4.12. Synthetic wide-angle seismograms for laterally discontinuous layers. ..................... 97 Figure4.13. Amplitude of the synthetic PmP phase. ..................................................................... 98 Figure4.14. Synthetic near-vertical shot gather for laterally discontinuous layers. ...................... 99 Figure4.15. Synthetic wide-angle seismograms for lamellae structure transition zone. ............. 100 Figure4.16. Amplitude of the synthetic PmP arrivals. ................................................................ 102 Figure4.17. Synthetic near-vertical seismograms for a lamellae transition zone. ....................... 103 Figure4.18. Schematic evolutionary model for the development of the Moho. .......................... 110  Figure 5.1. Tectonic map of the study area. ................................................................................ 125 Figure 5.2. Stacked, unmigrated seismic reflection profile ......................................................... 127 Figure 5.3. Record section for shot 1105. ................................................................................... 128 Figure 5.4. Map of estimated correlation lengths. ....................................................................... 132 Figure 5.5. Map of estimated power law exponent. .................................................................... 134 Figure 5.6. Map of average minimum least squares residuals..................................................... 135 Figure 5.7. Tri-modal velocity fields. .......................................................................................... 137 Figure 5.8. Observed and synthetic wide-angle seismograms..................................................... 139 Figure 5.9. Observed and synthetic near-vertical seismograms. ................................................. 140 Figure 5.10. Computed PmP envelope curves ............................................................................ 142 Figure 5.11. Best-fitting exponential curve. ................................................................................ 144 Figure 5.12. Envelope decay curves. ........................................................................................... 144 Figure 5.13. Trace-balanced F-X spectra maps. .......................................................................... 145 Figure 5.14. Trace-balanced peak frequency curves ................................................................... 146 Figure 5.15. Interpretation of the correlation length map............................................................ 148 Figure 5.16. Composite velocity model and corresponding synthetic seismograms. .................. 153  ix  Acknowledgements First and foremost, I would like to express my gratitude to my thesis advisor Ron Clowes. Ron has been a great supervisor providing overwhelming support and scientific insight. He guided and encouraged me since the first day I arrived at UBC. His reassuring and calm personality helped me enormously through the tough times and set a great example for me to follow in life. I would also like to thank Michael Bostock for his frequent input and enthusiastic support of my work. Thanks to Jim Mortensen and Felix Herrmann for serving on my supervisory committee and for their invaluable feedback as well. I have had the pleasure to collaborate closely with Vishal Kumar and Stefan Carpentier. Working with Vishal and Stefan introduced me to whole new and exciting fields and the outcome was more than gratifying. Many thanks to Phil Hammer who always inspired me in my research. Many of the ideas, which I developed in this thesis, were the result of countless stimulating and exciting discussions with Phil. I have had the privilege of crossing paths with many spectacular people in the department of Earth and Ocean Sciences at UBC. I would particularly like to acknowledge Ali Vaghri for a positive perspective in life, Kim Welford for all her help to get me started with my research, Goran Markovic, Peyman Moghaddam, Tom-Pierre Frappe, Brendan Smithyman, Rafael Sanguinetti, Dave Moore, Andrew Hamilton, Cecilia Li, Francis Jones, Sara Harrison, and Sukhi Hundal for making all these years I spent in the department a fantastic experience. Finally and most importantly, I wish to thank my parents, Naila and Hassan. This thesis would have not been written without their never-ending love and support. I wish my graduation will bring them pride and pleasure and will compensate for living so far away from each other. For making our life so wonderful, I thank my amazing wife Bannoura who shared the frequent frustrations and the rewards in this journey. I am forever grateful for her love, support and the enormous sacrifices she had to make.  x  ‫‪Dedication‬‬  ‫اﻟﻰ أﻣﻲ و أﺑﻲ‬ ‫ﻧﺎﺋﻠﺔ وﺣﺴﻦ‬  ‫‪xi‬‬  Statement of Co-authorship Chapter 2 was published with Ron Clowes, my supervisor, at UBC. Based on results from one of his students, Michael Bostock suggested the research program. I performed the research, synthetic modeling and interpretation. The finite-difference traveltime inversion algorithm was provided by Barry Zelt. Michael and Jean-Philippe Mercier assisted in several discussions about receiver functions. Jim Mortensen assisted in the geologic interpretation. I wrote the manuscript with input from Ron. Michael reviewed the manuscript at the final stage. Chapter 3 was prepared for publication with Vishal Kumar, now at CGGVeritas in Singapore, Ron Clowes and Felix Herrmann at UBC. Felix identified the research program. Vishal performed the technical research, algorithm development and I assisted him in synthetic modeling. I contributed the application to and interpretation of real data. I was fully involved in writing the manuscript in close collaboration with Vishal and Ron, and am the corresponding author with respect to the journal submission. Chapter 4 was prepared for publication with Ron Clowes at UBC. Ron and I identified the research program. I performed the research, synthetic modeling and interpretation. Jim Mortensen assisted in the geologic interpretation. Ron assisted in many discussions regarding the synthetic modeling. I wrote the manuscript with input from Ron. Chapter 5 was prepared for publication with Stefan Carpentier at the Swiss Federal Institute of Technology in Zurich and Ron Clowes at UBC. I identified the research program and performed the research. Stefan provided the algorithms and guidance. Stefan also assisted with the interpretation and in many discussions about stochastic analysis. I wrote the manuscript with input from Stefan and Ron.  xii  Chapter 1  Introduction, Objectives and Background 1.1 Introduction Seismology is the study of seismic, or elastic, waves that travel through the Earth in order to learn about its properties in the subsurface. For seismological studies of the whole Earth or large regions of it, naturally occurring earthquakes are usually the source of the seismic waves. For regional seismic studies of the crust and uppermost mantle and for exploration seismology in sedimentary basins, active sources (explosives, vibrating trucks, marine airguns, etc.) generate the seismic waves. The minute ground motions that result from the passage of the seismic waves are detected by arrays of seismometers or geophones and preserved on seismic recording systems, current ones being fully digital. The primary objective of active-source seismology is to extract information about the type and distribution of rocks and their properties in the subsurface. The primary parameter used to achieve this objective is the traveltime of the seismic waves from source to receiver but variations in amplitude, frequency, phase and wave shape can provide important additional information. Within active-source seismology, two complementary techniques, near-vertical incidence (NVI) reflection and refraction/wide-angle reflection (R/WAR), are principally used. As the name implies, in the NVI reflection method, the seismic energy propagates initially downward and then some of this energy is redirected upwards due to physical property changes in the subsurface, the paths of the seismic waves being near vertical. In other words, the offsets between the source and receivers are small relative to the depth of investigation. Dense shot (50-500 m) and receiver (25-100 m) spacings are typical in crustal studies. In the R/WAR method, the seismic energy also propagates away from the source but the location of the seismometers typically extends from near the source to offsets that are much larger (often 10x or more) than the planned depth of investigation.  1  For crustal studies, this geometry typically requires more sparse shot (20-100 km) and receiver (1-5 km) spacings. The seismic waves are refracted and reflected due to velocity and related variations in the subsurface; as a result of the geometry, their travel paths are much more horizontal than vertical. Between the two methods, near-vertical data have the higher resolution due to their frequency content (typically ~10-80 Hz for near-vertical and ~3-15 Hz for wide-angle data). For comparison, the vertical resolution of near-vertical signals in the lower crust is ~100-200 m based on the ¼ λ (wavelength) rule and the horizontal resolution is ~3 km based on the first Freznel zone (Widess, 1973; Sheriff and Geldart, 1995). On the other hand, the vertical resolution of wide-angle data is ~1.5-2.0 km and the horizontal resolution is ~30 km based typically on perturbation studies of model parameters (e.g., Zelt and Smith, 1992). As a result of these contrasting recording geometries and frequency contents, wide-angle and near-vertical data have complementary strengths: the former provide an estimate of the seismic velocity distribution and gross structures, whereas the latter provide detailed structural images. Following the discovery of a velocity discontinuity at the crust-mantle boundary in 1909 (Mohorovičić, 1910), albeit using earthquake recordings, seismic studies of the crust used the R/WAR technique almost exclusively for many decades (e.g., Steinhart and Meyer, 1961; German Research Group for Explosion Seismology, 1964). It has remained an important method for crustal and lithospheric studies to the present, as evidenced by continuing studies of this type throughout the world. During the 1950s, German and Russian scientists demonstrated that NVI reflected energy from within the crust could be used to map structural features with strong velocity contrasts (e.g., Dohr, 1957; Dohr and Fuchs, 1967; Zverev, 1967).  However, in the 1960s a group of  geophysicists at the University of Alberta in Edmonton were the first to demonstrate the efficacy and value of the NVI seismic reflection technique, as applied by industry for petroleum exploration, for detailed studies of Earth’s crust (Kanasewich and Cumming, 1965; Clowes et al., 1968). Following on this success, the Consortium for Continental Reflection Profiling (COCORP) was established in 1974 at Cornell University in Ithaca, NY to apply NVI reflection technology to map crustal structure at locations throughout 2  the United States (Oliver, 1986; Brown et al., 1986). The method is now used worldwide, both at sea and on land, and is indispensable for detailed studies of crustal/lithospheric structure. In this thesis, data from both NVI reflection and R/WAR studies along the same profile in the Northwest Territories are utilized to achieve more and better information about the Moho, or crust-mantle transition, and the uppermost mantle than was interpreted in the original studies of the data sets. The opportunities and advantages of using the two sets of complementary data are at the core of the research herein.  1.2 Theme and research objectives As noted in the Introduction, this thesis focuses on acquiring new results from existing NVI reflection and R/WAR data that were acquired as part of LITHOPROBE’s SNORCLE transect in northwestern Canada (Clowes, 1993). Subsequent to the interpretation of these data (Cook et al., 1999; Fernandez Viejo and Clowes, 2003) and new information about the sub-crustal lithosphere from a UBC study of teleseismic data along the same profile (Mercier et al., 2008), my supervisor and I recognized that an alternate, and likely improved, interpretation could be derived through consideration of the joint data set. Consequently, the objective of the first part of my study is a reanalysis and reinterpretation of these data as they applied to better definition of structure in the subcrustal lithosphere. Thereafter, my supervisor and I determined that detailed quantitative studies of the Mohorovičić discontinuity using the combined data sets should represent the balance of my thesis work. The basis of this decision was two-fold. First, a number of qualitative studies of Moho characteristics had been undertaken using LITHOPROBE data (e.g., Hammer and Clowes, 1997; Cook, 2002; Eaton, 2006) but none had addressed the analysis quantitatively. Indeed, only a few similar studies had been carried out by other researchers and these relied primarily on wide-angle data (e.g., Carbonell et al., 2000). Secondly, the Moho is one of the most prominent worldwide discontinuities on Earth, yet its nature and characteristics still are not well known or understood.  3  Is the Moho a finite-thickness, laterally and vertically variable transition zone that can be represented by a lamellae structure and does it have a stochastic heterogeneity distribution? Or is it an interface that represents a change in the scale length of heterogeneity between the lower crust and upper mantle? What are the origins for the deep seismic reflections both at Moho levels and in the uppermost mantle? Does the seismic wavefield recorded on the near-vertical and wide-angle seismic data hold the keys for such questions? These are some of the questions for which we are seeking answers throughout this thesis. Unraveling the detailed structure of this transition zone and the way in which it forms and evolves may provide relevant information that further helps Earth scientists understand some major tectonic processes such as crustal growth, continental accretion and delamination. However, the primary focus in this thesis is quantitative means by which characteristics of the Moho transition zone can be established, not the tectonic processes that might derive from these characteristics. Nevertheless, some discussion of such processes is included. Since this thesis is written in the form of a collection of publications and manuscripts intended for publication, there will inevitably be some overlap between the chapters, particularly in the introductory information that is provided.  1.3 Background Much of our knowledge about the physical structure and the chemical composition of the Earth’s lower crust and upper mantle is based on the interpretation of seismic data, simply because of the difficulty of directly sampling significant parts of them. Unlike samples obtained from either xenoliths or exhumed lower crustal cross sections (e.g., Ivrea-Verbano zones in northern Italy), from which assessment of their applicability to the lower crust/upper mantle as a whole is difficult, seismic methods provide a means of indirectly measuring their structure and composition in situ across large sections.  4  The existence of a discontinuity at the crust-mantle boundary was first recognized by Mohorovičić (1910) based on the interpretation of regional earthquake recordings (i.e., teleseismic data) in Eastern Europe. Mohorovičić deduced that P-wave arrivals often lay on two separate travel-time branches, which he attributed to a rapid change in seismic velocities from about 5.65 km/s to over 7.5 km/s depth near 54-km depth in that area. The discontinuity has since been known as the Mohorovičić discontinuity, abbreviated as ‘Moho’ (Jarchow and Thompson, 1989; Cook et al., 2010). Subsequently, wide-angle surveys have been used to delineate the Moho globally with a greater detail than earthquake studies (e.g., Mooney et al., 1998; Perry et al., 2002). These surveys showed that there is significant variation in terms of velocity contrast and thickness of the discontinuity especially within continental regions. In order to provide a consistent definition, Steinhart (1967) proposed that the ‘refraction Moho’ is the depth at which Pwave velocities increase rapidly or discontinuously to 7.6-8.6 km/s. In the absence of a steep velocity gradient, the refraction Moho is interpreted as the level at which the Pwave velocity exceeds 7.6 km/s. As mentioned previously, it was not until the late 1960s that useful near-vertical reflections from the lower crust became available. The Moho was expected to provide one of the strongest near-vertical reflections, yet reflections from the vicinity of the refraction Moho have proven to be highly variable and, at times, elusive features. In some regions there are no reflections, in others there are distinct, multicyclic bands of nearly continuous reflections, and generally variable characteristics between these two extremes (e.g., Hammer and Clowes 1997; Cook, 2002). However, the one dominant feature of virtually all good reflection profiles, regardless of the detailed interpretation, is that the crust is substantially more reflective than the mantle. As a result, the ‘reflection Moho’ has been defined as the deepest, high amplitude, broad reflection or group of reflections that can usually be followed in a piecewise continuous manner at travel times comparable to the crustal thickness (e.g. Klemperer et al., 1986). In most cases where modern data are available, the seismic refraction and reflection Moho appears to correspond in depth to within the resolution of the two methods (Mooney and Meissner, 1992; Cook et al., 2010). 5  From a petrological perspective, the crust-mantle boundary represents either 1) a compositional change from mafic rocks above the transition to peridotitic rocks below (e.g., Green and Ringwood, 1972) or 2) a phase change from mafic granulite above to eclogitic rocks below (e.g., Ito and Kennedy, 1971). As the technological advances in both geophysical (i.e., seismic) and geological (i.e., petrologic) experiments allowed for a better understanding of the transition, it became clear that the Moho (geophysical representation) and the crust-mantle transition (geological change) may not always represent the same boundary (e.g., Griffin and O’Reilly, 1987; Xu et al., 1996). However, for the purpose of this thesis the two terms, the crust-mantle boundary and the Moho have been used interchangeably to refer to the boundary detected remotely by near-vertical and wide-angle seismic methods. The complex nature of the near-vertical signal reflected from Moho levels and the rich coda following the onset of the wide-angle Moho arrivals (PmP phase) can be explained by a model in which the crust-mantle boundary is a finite thickness, heterogeneous transition zone. The transition consists of thin, high- and low-velocity layers (lamellae) that are anastomosing rather than continuous and flat (e.g., Braile and Chang, 1986; Carbonell et al., 2002). Another model that can explain the aforementioned seismic characteristics is based on the heterogeneity distribution found in the lower crustal surface exposures and borehole logs. These heterogeneities are not confined to one typical scale length, but extend over a wide range of scales (fractal) and seem random in their distribution. The stochastic crustal heterogeneities usually exhibit spatially selfaffine, or power-law scaling, often described by the von Karman autocovariance function (Goff and Jordan, 1988; Holliger and Levander, 1992; Holliger, 1996). The von Karman function is parameterized by 2 fractal parameters: correlation length, which is an upper limit for the scale invariance in the heterogeneity, and Hurst number, which is an exponent that controls the degree of scale invariance in heterogeneity below the correlation length scale (Carpentier and Roy-Chowdhury, 2007). The Moho transition is also the site where probably one of the most pervasive and extensive tectonic feature of the Earth, “lithospheric subduction”, occurs. Subduction zones are the geodynamic systems that produce island arcs, the building blocks of the 6  continental crust. They deliver raw materials to the subduction factory, where sediments, oceanic lithosphere, and seawater return to and reequilibrate with ambient mantle, initiating melting and incidentally creating continental crust (Stern, 2002). Thus, subduction zones represent one part of our planet’s recycling system. Such processes associated with subduction impact the Moho significantly, re-shaping it and affecting its evolution over time. Traditionally, two types of subduction zones have been identified based on the age of the subducting lithosphere. Mariana type (old lithosphere), which is characterized by steep angles of dip and Chilean type (young lithosphere), which is characterized by shallow dipping slabs. However, more recent studies indicate that the main parameter controlling the angle of dip is the overriding plate absolute motion (e.g., Lallemand et al., 2005; Heuret and Lallemand, 2005). One important aspect about shallow subduction zones, which is the theme of the first part of this thesis, is that they may have played a major role in the generation and stabilization of new continental crust (e.g., Helmstaedt and Schulze, 1989; Abbott, 1991).  1.4 Thesis outline In Chapter 2, we re-examine near-vertical incidence and refraction/wide-angle reflection seismic data recorded as part of LITHOPROBE studies in the Paleoproterozoic-Archean domains of Canada`s Northwest territories. The seismic data reveal strong, subhorizontal features in the uppermost mantle from one of the best imaged relict subduction zones within very old lithosphere. Previous interpretations of the seismic data did not correlate such features with the dipping slab. By proper consideration of the crooked line 3D acquisition geometry and careful use of 2D ray tracing and 3D forward and travel-time inversion algorithms, we propose a new interpretation for the deep subhorizontal reflections. We also investigate the origin of these reflections, their relation with other upper-mantle reflectors, the characteristic features of shallow subduction zones and finally construct a 3D perspective for the east-dipping, subducted slab. Chapter 3 arises from the need to increase the signal-to-noise ratio of reflection data to assist quantitative studies of such data. It presents a new method for the suppression of incoherent noise, present in the seismic signal, based on the curvelet transform. We make 7  use of the parsimonious representation of seismic data in the curvelet domain to perform the noise attenuation while preserving the coherent energy. First, we demonstrate the applicability of the method on noisy synthetic data then we extend this to deep crustal, near-vertical seismic data (same as in Chapter 2) where the signal-to-noise ratio is low. The curvelet denoising is applied to both pre-stack shot gathers and post-stack data. Comparing the results with that of similar coherency filtering techniques (e.g., F-X deconvolution), curvelet denoising outperforms the latter by attenuating incoherent noise with minimal harm to the signal. In Chapter 4, we examine the detailed structure of the Moho and the processes associated with its formation and evolution. Sharp, multicyclic bands of reflections at Moho levels have been imaged on near-vertical data beneath the Great Bear magmatic arc. They show a remarkably flat Moho and do not reflect the complex tectonic history of the Wopmay orogen. The wide-angle data is characterized by a prominent PmP phase followed by long coda. We calculate the synthetic seismograms of 1- and 2-D velocity models using 1-D reflectivity and 2-D finite-difference wave propagation algorithms to determine the extent to which these models can replicate the wavefield characteristics noted above. Two end models that can explain the observed wavefield are suggested for the Moho transition, which formed as a combination of thermal and mechanical processes. Chapter 5 extends the studies in Chapter 4 by utilizing an objective statistical approach to further constrain the nature of the Moho transition zone. It represents an inverse approach to the problem, in contrast to the forward modeling approach of Chapter 4. The von Karman stochastic distribution of heterogeneities conserved in the recorded near-vertical seismic data beneath the Great Bear magmatic arc is used to extract the stochastic parameters. The estimated lateral correlation lengths and the Hurst number representing the aforementioned parameters are then used to construct tri-modal 2-D stochastic velocity fields for the lower crust, the Moho transition zone and the uppermost mantle. To further constrain the stochastic parameters, synthetic wide –angle and nearvertical seismograms are calculated from the velocity fields and compared to the observed seismic data using a 2-D finite-difference wave propagation algorithm. Our 8  preferred velocity model contains three stochastic regions; the lower crust, the Moho transition zone, and the uppermost mantle. The Moho transition is characterized by the highest correlation length values, whereas the uppermost mantle is characterized by the lowest values.  9  Bibliography Abbott, D. 1991. The case for accretion of the tectosphere by buoyant subduction. Geophysical Research Letters, 18, p. 585– 588. Braile, L.W., and Chiang, C. S., 1986. The continental Mohorovičić discontinuity results from near-vertical and wide-angle seismic reflection studies. In Barazangi, M., Brown, L., eds., Reflection Seismology. A Global Perspective: Washington DC, American Geophysical Union Geodynamic Series, 13, p. 257-272. Brown, L., Barazangi, M., Kaufman, S., and Oliver, J., 1986. The first decade of COCORP: 1974-1984. In Barazangi, M., Brown, L., eds., Reflection Seismology: A Global Perspective. Washington DC, American Geophysical Union Geodynamic Series, 13, p. 107-120. Carbonell, R., Gallart, J., and Perez-Estaun, A., 2002. Modeling and imaging the Moho transition: the case of the southern Urals. Geophysical Journal International, 149, p. 134-148. Carpentier, S. and Roy-Chowdhury, K., 2007. Underestimation of scale lengths in stochastic fields and their seismic response: a quantification exercise. Geophysical Journal International, 169, p. 547–562. Clowes, R.M. (ed.), 1993. LITHOPROBE Phase IV Proposal – Studies of the Evolution of a Continent. LITHOPROBE Secretariat, The University of British Columbia, Vancouver, B.C., 290 p. Clowes, R.M., Kanasewich, E.R., and Cumming, G.L., 1968. Deep crustal seismic reflections at near-vertical incidence. Geophysics, 33, p. 441-451. Cook, F.A., 2002. Fine structure of the continental reflection Moho. Geological Society of America Bulletin, 114, p. 64-79. Cook, F. A., White, D. J., Jones, A. G., Eaton, D. W., Hall, J., and Clowes, R. M., 2010. How the crust meets the mantle: Lithoprobe perspectives on the Mohorovičić discontinuity and crustmantle transition. Canadian Journal of Earth Sciences, 47, p. 315-351. Dohr, G., 1957. Ein Beitrag der Reflexionsseismik zur Erforschung des teiferen Untergrundes. Geologische Rundschau, 46, p. 17-26. (Translated by T. Gretener.) Dohr, G., and Fuch, K., 1967. Statistical evaluation of deep crustal reflections in Germany. Geophysics, 32, p. 951-967. 10  Eaton, D. W., 2006. Multi-genetic origin of the continental Moho: insights from LITHOPROBE. Terra Nova, 18, p. 34-43, doi: 10.1111/j.1365-3121.2005.00657. Fernandez Viejo, G., and Clowes, R.M., 2003. Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay orogen, northwestern Canada, from a LITHOPROBE refraction/wide-angle reflection survey. Geophysical Journal International, 153, p. 1-19. German Research Group for Explosion Seismology, 1964. Crustal structure in western Germany. Ztschr. f. Geophys., 30, p. 209-234. Goff, J. and Jordan, T., 1988. Stochastic modeling of seafloor morphology: inversion of sea beam data for second-order statistics. Journal of Geophysical Research, 93, p.13589–13608. Green, D. H., and Ringwood, A. E., 1972. A comparison of recent experimental data on the gabbro-garnet granulite-eclogite phase transition. Journal of Geology, 80, p. 277-288. Hammer, P. T. C., and Clowes, R.M., 1997. Moho reflectivity patterns - a comparison of Canadian LITHOPROBE transects. Tectonophysics, 269, p. 179-198. Helmstaedt, H., and Schulze, D. J. 1989. Southern African kimberlites and their mantle sample: Implications for Archean tectonics and lithosphere evolution. In Kimberlites and Related Rocks, Blackwell , Cambridge, Mass., p. 358– 368. Heuret, A., and Lallemand, S. 2005. Plate motions, slab dynamics and back-arc Deformation. Physics of the Earth and Planetary Interior, 149, p. 31-51. Holliger, K. and Levander, A.,1992, A stochastic view of lower crustal fabric based on evidence from the Ivrea Zone. Geophysical Research Letters, 19, p. 1153–1156. Holliger, K., 1996. Upper-crustal seismic velocity heterogeneity as derived from a variety of pwave sonic logs. Geophysical Journal International, 125, p. 813–829. Ito, K., and Kennedy, G.C., 1971, An experimental study of the basalt-garnet granulite-eclogite transition, in Heacock, J.G., ed., The Earth beneath the continents. American Geophysical Union Geophysical Monograph 14, p. 303–314. Jarchow, C.M., and Thompson, G.A., 1989. The nature of the Mohorovicic discontinuity. Annual Review of Earth and Planetary Sciences, 17, p. 475-506. Kanasewich, E.R., and Cumming, G.L., 1965. Near-vertical-incidence seismic reflections from the ‘Conrad’ discontinuity. Journal of Geophysical Research, 70, p. 3441-3446.  11  Klemperer, S.L., Hauge, T. A., Hauser, E. C., Oliver, J. E., and Poterr, C. J., 1986. The Moho in the northern Basin and Range province, Nevada, along the COCORP 400N seismic reflection transect. Geological Society of America Bulletin, 97, p. 603-618. Lallemand, S., Heuret, A., and Boutelier, D. 2005. On the relationships between slab dip, backarc stress, upper plate absolute motion, and crustal nature in subduction zones. Geochemistry Geophysics Geosystems, doi:10.1029/2005GC000917. Mercier, J.P., Bostock, M.G., Audet, P., Gaherty, J.B., Garnero, E.J., Revenaugh, J. 2008. The teleseismic signature of fossil subduction: Northwestern Canada. Journal of Geophysical Research., 113, B04308, doi: 10.1029/2007JB005127. Mohorovičić, A., 1910, Das Beben vom 8. X. 1909, Jahrbuch des meteorologischen Observatoriums in Zagreb, 9/4, p. 1-63. Mooney, W. D. , and Meissner, R., 1992. Multi genetic origin of crustal reflectivity: a review of seismic reflection profiling of the continental lower crust and Moho. In Arculus, D. M., and Kay, R. W., eds, Continental Lower Crust, Elsevier, p. 45-79. Mooney, W. D., Laske, G., and Masters, T. G., 1998. CRUST 5.1: A global crustal model at 5˚ x 5˚. Journal of Geophysical Research, 103, p. 727-747. Oliver, J., 1986. A global perspective on seismic reflection profiling of the continental crust. In Barazangi, M., Brown, L., eds., Reflection Seismology. A Global Perspective: Washington DC, American Geophysical Union Geodynamic Series, 13, p. 1-3. Perry, H. K. C., Eaton, D. W. S., and Forte, A. M., 2002. LITH5.0: a revised crustal mode for Canada based on Lithoprobe results. Geophysical Journal International, 150, p. 285-294. Sheriff, R. E., and Geldart, L. P., 1995. Exploration Seismology, Cambridge University Press, New York, NY, 592 p. Steinhart, J. S., 1967, Mohorovičić discontinuity, In Runcorn, S. K., ed., International Dictionary of Geophysics, Volume: 2, Pergamon Press, New York, p. 991-994. Steinhart, J.S., and Meyer, R.P., 1961. Explosion studies of continental structure. Carnegie Inst. Wash. Publ. 622, Washington. The Kirby Lithographic Company, Inc. 409 p. Stern, R., J. 2002. Subduction zones. Reviews of Geophysics, 40, doi:10.1029/2001RG000108.  12  Telford, W.M., Geldart, L.P., and Sheriff, R.E., 1990. Applied Geophysics (2nd edition), Cambridge University Press, 792 p. Zelt, C. A., and Smith, R. B., 1992. Seismic traveltime inversion for 2-D crustal velocity structure. Geophysical Journal International, 108, p. 16-34. Zverev, S.M. (ed.), 1967. Problems in Deep Seismic Sounding. New York: Consultants Bureau, 166 p.  13  Chapter 2  Paleoproterozoic subduction in NW Canada from near-vertical and wideangle reflection data1 2.1 Introduction Westward growth of the continental lithosphere in northwestern Canada during the Paleoproterozoic era is the result of continental breakup and a series of major orogenic episodes, as established by geological studies (e.g., Hoffman 1989). The tectonic processes included formation of a passive margin on the western side of the Archean Slave craton (e.g., Coronation Supergroup), accretion of exotic terranes (e.g., Hottah and Fort Simpson terranes) and development of an Andean-style magmatic arc through subduction (e.g., Great Bear magmatic arc); see Fig. 2.1. As part of the multidisciplinary studies associated with LITHOPROBE’s Slave-NORthern Cordillera Lithospheric Evolution (SNORCLE) transect, geophysical studies were undertaken to provide information in the third dimension – depth. These included acquisition of near-vertical incidence (NVI) seismic reflection (Cook et al. 1998, 1999), seismic refraction/wide-angle reflection (R/WAR; Fernandez Viejo and Clowes 2003) and, more recently, broadband teleseismic data along the only road in the region (Fig. 2.1; Mercier et al. 2008), plus magnetotelluric data along the road and throughout much of the Slave craton (Jones et al. 2005). This contribution re-examines the NVI and R/WAR data sets in light of the new data and interpretation from the teleseismic studies. 1  A version of this chapter has been published: Oueity, J., and Clowes, R. M., 2010. Paleoproterozoic  subduction in NW Canada from near-vertical and wide-angle reflection data. Canadian Journal of Earth Sciences, 47, p. 35-52.  14  NVI data revealed very distinct reflections from within the upper mantle (Fig. 2.2; Cook et al. 1998, 1999). Near the Fort Simpson-Hottah boundary, a parallel pair of reflectors at ~9 and 11 s two-way traveltime dip eastward from the Moho (~33 km) over a distance of ~60 km to times of 14 and 17 s. They can be correlated with a similar pair of reflectors that reach times of 21 and 23 s about 80 km farther downdip. The reflectors have been interpreted as the top and bottom of a slab of subducted oceanic crust preserved since the ca. 1.84 Ga convergence of the Fort Simpson terrane with the Hottah terrane, two components of the Wopmay orogen (Cook et al. 1998, 1999). Fifty kilometers eastward, another pair of reflectors is observed at 20 and 22 s. They extend continuously farther east at about the same times for about 100 km before dipping shallowly to times of 25 and 28 s (depths of 85-100 km).  Figure 2.1. Tectonic map of the study area and location of SNORCLE Line 1. Stars identify shot locations for the wide-angle experiment; some are identified by numbers. Inverted triangles represent the teleseismic station locations. Sediments of the Phanerozoic Western Canada Sedimentary Basin overlie the Precambrian domains west of the long dashed line. Short dashed black lines show political boundaries. CS - Coronation Supergroup, SD - Sleepy Dragon, GSLsz - Great Slave Lake shear zone, AB - Alberta, BC British Columbia, NWT - Northwest Territories, YK - Yukon. Inset shows location of map within Canada.  15  Figure 2.2 The originally processed, coherency filtered, migrated seismic data along SNORCLE Line 1, modified from Cook et al., (1999). Feature denoted by the letter (S) refers to upper mantle reflections interpreted as the top and bottom of a relict subducted slab.  From the refraction experiment, wide-angle reflections below the Moho (the PM1P phase) were identified on many of the data sections (Fig. 2.3; Fernandez Viejo and Clowes 2003). After establishing a crustal velocity model, simple 2-D traveltime modeling that assumed a horizontal reflector showed that the PM1P phase could be generated from a laterally extensive upper mantle reflector at a depth of ~75 km, consistent with the general depth from which the NVI mantle reflections were generated (see Fig. 12 in Fernandez Viejo and Clowes 2003).  Figure 2.3 Record section for shot 1101. Data are bandpass-filtered from 3-13 Hz, and displayed as tracenormalized amplitudes. Identified phases are described in the text. Enlargement of phase PM1P is shown in the inset; grey circles show picks.  16  Figure 2.4. a) Location of leg A, part of CANOE experiment (Mercier et al., 2008), which is coincident with Line 1. Triangles represent the teleseismic station locations. YK - Yellowknife, NB - Nahanni Butte. b) Polarity-corrected transverse receiver functions. Red and blue denote positive and negative polarity respectively. The image is superimposed on top of the reflection data (Cook et al., 1999). Dashed lines represent the interpreted slab from the reflection data. The slab lies directly above an anisotropic layer representing the shallowmost part of mantle material and associated with positive and negative polarity (modified from Mercier et al., 2008).  In the original interpretation of the NVI data, the pair of sub-horizontal reflections at 20 and 22 s were not considered part of the subducted slab. They represented a different feature, one associated with imbrication of the mantle lithosphere (Cook et al. 1998, 1999). While this is a valid interpretation, we note that the reflections become subhorizontal approximately where the road along which they were recorded veers from an east-west direction to one that is almost north-south. Recent results from analysis of the teleseismic data clearly identify the dipping slab but indicate that its dip might become more shallow at greater depths (Fig. 2.4; Mercier et al. 2008). These points raise the question: do the sub-horizontal reflectors represent the same slab as the dipping ones? The crooked line acquisition geometry, which involved 3-D aspects, allowed us to reanalyze the NVI and R/WAR datasets using both 2-D and 3-D procedures. These include 3-D forward wavefield modeling using the Kirchhoff integral formula, and 3-D finite-  17  difference traveltime inversion. Through application of those procedures, we investigate the relation between the descending slab and the subhorizontal reflectors, the continuity of the slab and its three-dimensional geometry and structure.  2.2 Tectonic setting The Wopmay orogen in northwest Canada is an amalgamation of domains that developed during the Calderian orogeny between 1.92 Ga and 1.84 Ga (Fig. 2.1; Hoffman 1989). The orogen comprises several major north-south trending tectonic elements, the Coronation Supergroup, the Hottah terrane, the Great Bear magmatic arc and the Fort Simpson terrane. On the east, the Coronation Supergroup is a continental passive margin depositional prism developed between 1.92 and 1.90 Ga following the cessation of lithospheric accretion in the Slave province. The Coronation Supergroup was detached from its basement, thickened and thrust eastward over the western edge of the Slave craton as a result of the Hottah terrane collision (Hoffman and Bowring 1984). The Coronation margin deposits record the transition from the opening of an ocean basin to subsequent collision of its margin. To the west, the Hottah terrane is an arc-bearing microcontinent that collided with the Coronation margin ca. 1.90 Ga (Hildebrand and Bowring 1999). The Hottah terrane consists  of  Paleoproterozoic  sedimentary  and  intermediate  volcanic  rocks  metamorphosed to amphibolite grade and intruded by calc-alkaline granitic plutons during the period 1.940-1.902 Ga (Hildebrand et al. 1987). The accretion of the Hottah terrane onto the Slave craton is considered to be the result of west-dipping ocean basin subduction since no coeval arc magmas are found on the western Slave. To the west, another basin with its associated oceanic plate was converging against the western Hottah, generating an east-dipping subduction zone below the Hottah terrane and the Slave craton and leading to the formation of the Great Bear magmatic arc. The 100 km-wide Great Bear magmatic arc, which consists of volcanosedimentary sequences and plutonic rocks, lies between the Coronation Supergroup and the Hottah terrane (Hildebrand et al. 1987). The rocks of the arc were deposited on both the Hottah  18  terrane and deformed rocks of the Coronation margin during the period between 1.875 to 1.843 Ga (Hildebrand et al. 1987; Gandhi et al. 2001). Aeromagnetic data indicate that rocks of the Great Bear arc extend from north to south beneath the younger Proterozoic and Paleozoic cover for a total length of about 900 km (Coles et al. 1976). Seismic reflection data from SNORCLE Line 1 showed that the Great Bear arc is about 2.0 – 4.5 km thick in the region of the survey and does not seem to disrupt the older structure beneath it (Cook et al. 1999). In their geochronological study, Gandhi et al., (2001) concluded that the source for magmatism in the southern part of the arc, where Line 1 is located, is local. The only explanation for the intact older structure beneath the arc (as seen on the reflection data) is that the arc magmas migrated through the crust along preexisting zones of weakness without disrupting the structure (Cook et al., 1999). After the cessation of magmatic activity, the entire zone was affected by right-lateral oblique slip resulting in a system of northeast-trending faults (Gandhi et al. 2001). A subsequent event of compression within the Great Bear – Hottah – Slave assembly is attributed either to the collision of the Fort Simpson terrane on the west side of the Hottah terrane (Hildebrand et al. 1987) or the docking of the Nahanni domain (Hoffman 1989). The timing for the collision of the Fort Simpson terrane with the Hottah terrane occurred after the formation of Fort Simpson magmatic rocks ca. 1.845 Ga (Villeneuve et al. 1991) and pre-1.71 Ga mafic intrusions of sedimentary rocks (Abbott et al. 1997). Based on the geometry of reflectors seen on the seismic data, Cook et al. (1999) suggested that the upper crust of the Fort Simpson terrane was thrust over an accretionary wedge formed by the collision. At the same time, the lower crust of the Fort Simpson terrane was attached to the subducting crust below the wedge and was carried into the mantle. Following the terminal collision of the Fort Simpson with the Hottah terrane a period of lithospheric extension was responsible for the formation of the Fort Simpson basin. Seismic reflection data reveal that the basin had at least 20 km of subsurface relief (Cook et al. 1999). The entire plate assembly was affected by a number of crustal extension phases that continued to Neoproterozoic time (Thompson et al. 1987). The Proterozoic rocks of the Wopmay orogen along SNORCLE Line 1 are overlain by Phanerozoic sedimentary rocks  19  that thicken from their feather edge near the eastern boundary of the Great Bear magmatic arc to about 1000 m over the Nahanni domain (Cook et al. 1999).  2.3 Seismic experiments and data The LITHOPROBE SNORCLE transect investigated the various processes involved in the westward growth of North America from the Archean to the present (Cook and Erdmer 2005). The first NVI reflection experiment, the 725-km-long Line 1 acquired in 1996, crossed the southwest Slave craton and Wopmay orogen, along highways from east of Yellowknife to west of Fort Simpson (Fig. 2.1). Acquisition involved 4 to 5 vibroseis trucks with a 90 m shot point spacing, 60 m receiver spacing and a total of 404 channels/record. Sweep frequencies are 10-80 Hz and the sample rate is 4 ms. Complete field acquisition parameters are given in Cook et al. (1999). In the original analysis, the reflection data were processed following standard procedures (Cook et al. 1999). The NVI data have relatively high frequency content (vibroseis sweeps between 10 and 80 Hz) and generate high resolution crustal seismic images, from which structure and other characteristics are interpreted. In theory, vertical resolution is given by λ/4, where λ is the wavelength. Horizontal resolution is given by the radius of the first Freznel zone, r = λh / 2 , where h is the depth. However, in reality both resolutions are strongly influenced by the overlying velocity structure and thus are not as good as theory predicts. At 70 km depth where the seismic velocity is ~8.3 km/s (Fernandez Viejo and Clowes 2003) and the dominant frequency associated with the upper mantle reflectors is ~25 Hz, the nominal vertical resolution is ~80 m and the nominal horizontal resolution is ~3.4 km. The processed seismic images showed highly reflective crust, a generally transparent upper mantle and a relatively flat Moho across the different tectonic regimes (Fig. 2.2). Within the upper mantle, a number of reflectors are recognized. The most prominent ones (S) have been interpreted as the top and bottom of a relict subducted slab resulting from the collision between the Fort Simpson terrane and the Hottah terrane (Cook et al. 1998, 1999).  20  In 1997, the Line 1 R/WAR experiment was undertaken (Fig. 2.1). A total of 594 seismographs, deployed with an average spacing of 1.2 km, recorded signals from 13 explosive shots ranging in size between 400 and 10,000 kg. Two of these shots were located 200 km off either end of the line to provide information about the deep lithospheric mantle structure. Further details of the field acquisition can be found in Fernandez Viejo and Clowes (2003). Interpretation of the R/WAR data provided information on the seismic velocity structure, but the resolution associated with such interpretation is much less than that associated with NVI data. Based on procedures discussed by Fernandez Viejo and Clowes (2003), the estimated uncertainty in velocities in the upper mantle is 0.15 km/s and the associated lateral resolution is > 80 km. For the WAR data, traveltime modeling indicates vertical resolution at upper mantle depths in the order of 2 to 3 km. Interpretation of the R/WAR data showed average crustal velocities of 6.4-6.5 km/s and upper mantle velocities ranging from 8.1 km/s at 33 km depth to 8.4 km/s at 75 km depth. Along the east-west section of the line, high lower crust velocities (~7.1 km/sec) are located in the transition zone between the Fort Simpson and the Hottah terrane (Fernandez Viejo and Clowes 2003). In the same region, the upper mantle is characterized by anomalously low velocities (~7.5 km/s). Fernandez Viejo and Clowes (2003) attributed these characteristics to the effects of serpentinization of mantle peridotites fluxed by water derived from dehydration of the subducted slab interpreted from the NVI data. Coherent secondary phases at offset distances greater than 250 km on many shot gathers (e.g., PM1P in Fig. 2.3) were interpreted as wide-angle reflections from horizontal upper mantle reflectors at about 75 km depth for the PM1P phase and 100 km depth for another phase. A normal-moveout-corrected and stacked section of shot gathers on which the secondary phases were observed showed subhorizontal reflections between 22 s and 25 s that extended laterally over ~150 km below the Hottah terrane and Great Bear magmatic arc, and more quasi-continuous beneath regions to the east and west (Fernandez Viejo et al., 1999). The stacked wide-angle reflections below the Hottah and Great Bear terranes corresponded remarkably well with the subhorizontal reflections observed on the Line 1 NVI data. 21  The broad-band, three-component teleseismic CAnada NOrthwest Experiment (CANOE) took place between 2003 and 2005 (Mercier et al. 2008). Northeasterly leg A of the array coincided partly with Line 1. Over the Wopmay orogen, stations were spaced more closely than elsewhere to elucidate further the subducted slab interpreted from the NVI data (Fig. 2.1). Teleseismic methods generate images and information that are complementary to those derived from active-source methods. However, since they are characterized by lower frequencies (0.2–2.0 Hz), they inherently have lower resolution. Mercier et al. (2008) generated teleseismic receiver function images for the Leg A stations that were superimposed on the NVI data showing the subducted slab (e.g., Fig. 2.4). Both radial and transverse components clearly reveal the response of the subducted slab. Based on modeling and numerical simulations of both components, a 10-km-thick anisotropic layer was identified. On the transverse component, it can be traced from 30 km to at least 90 km depth over a horizontal distance of ~150 km (Fig. 2.4). However, the dipping feature does not coincide with the NVI-interpreted subducted crust. Instead, it parallels the bottom of the subducted crust imaged by the NVI data. As such, the layer is interpreted as a lid of anisotropic shallowmost mantle material (Mercier et al. 2008).  2.4 Analysis procedures and results 2.4.1 NVI reflection data The original processing of the NVI reflection data used parameters optimal for crustal images. The results were of high quality and enabled a detailed interpretation of crustal structure (Cook et al. 1999). They also allowed interpretation of Paleoproterozoic crustal rocks within the upper mantle, these being associated with a subducted crustal slab that occurred as the Fort Simpson terrane collided with the Hottah terrane. In an attempt to provide higher resolution images associated with the slab, we reprocessed the data with modified parameters (Table 2.1). In particular, bin sizes were enlarged to 90 m by 2000 m to increase the nominal fold and the minimum trace offset included in any bin was 2500 m to avoid masking effects of ground roll at near offsets. Figure 2.5 shows the  22  original (a) and new (b) seismic sections from 15 to 30 s (approximately 50 to 110 km depth). The new image has better signal-to-noise ratio and the lateral continuity of the reflectors is more pronounced. Another important observation is the presence of a deeper band of reflectivity below 25 s (reflections D in Fig. 2.5c), which extends subhorizontally along the entire section. The two dark bands at ~60 km distance and from 23 to 30 s are an artifact of processing; we were not able to eliminate them without compromising the seismic image. Based on these results we provide two possible interpretations regarding the base of the fossil crustal slab; the top surface is the same for both (Figs. 2.5b, 2.5c). In the first case, we consider that the long-dash line (A – B in Fig. 2.5c) represents the base of the slab. This is based on similar reflectivity within the interpreted slab and NNEdipping reflectivity that appears to truncate sub-horizontal reflectivity between 230 and 250 km distance. This interpretation also maintains an approximately constant thickness for the slab across the seismic section. In the second case, we consider that the reflectivity at A and C (Fig. 2.5b) represents the base of the slab (short-dash line, A – C in Fig. 2.5c) because they are the strongest reflections sub-parallel to and below the top of the slab. From C to the end of the section, a weak dipping reflection sub-parallel to the top of the slab can be discerned, as indicated by the short-dash line. Approximately at the location of C, this interpreted lower interface merges with zone of sub-horizontal reflectivity, D. One difficulty with this interpretation is that it involves an increase in the thickness of the slab east of A from ~10 km to ~15 km, a feature for which we don’t have an explanation. For this reason, we prefer the interpretation following A – B as representing the base of the slab, but the data are certainly allowable of either case. The deep band of reflectivity (D in Figs. 2.5b, 2.5c) probably represents another reflective zone within the upper mantle. However, it is difficult to associate this zone with any tectonic event because there is no obvious tie to shallower features. At a depth of ~100 km, this band is coincident with a wide-angle reflector interpreted from the R/WAR survey along Line 1 (horizon J in Fig. 12b of Fernandez Viejo and Clowes 2003). This feature is discussed later in the paper.  23  Figure 2.5. Comparison between original and re-processed section of NVI data. Inset shows the location of seismic section (grey line) along SNORCLE Line 1 (black line). YK - Yellowknife, NB - Nahanni Butte. (a) Original coherency filtered, migrated data (Cook et al., 1999). (b) The re-processed data (this study) with modified parameters (see text) showing very distinct upper mantle reflectors. Data are coherency filtered; no migration. The two dark bands at about 60 km distance and from 23 to 30 s are processing artifacts. (c) Two possible interpretations: east-dipping slab of crustal rocks with constant thickness (long dashed line) or with variable thickness (short dashed line). Deeper reflections may represent a second mantle reflector. Features identified by letters are explained in the text.  24  Parameter Crooked line geometry  Static correction Trace balance Bandpass filter Deconvolution Automatic gain control NMO correction on common midpoint gathers Stack Inverse Q filter F-X deconvolution Coherency filter Plot  Value 90 m × 2000 m bins 2.5 km minimum offset included 13 km maximum offset included Average amplitude scaled to be 1.0 5 - 35 Hz Filter length = 150 ms; gap length of 2.0 ms 500 ms window Velocities from 5.8 - 8.4 km/s nominally 305 fold Q = 500 upper crust, Q = 1000 lower crust/upper mantle sliding window 0.1 ms/meter maximum slowness  Table 2.1. Parameters for reprocessing of NVI data.  2.4.2 Two-dimensional ray-tracing We implemented a 2-D forward ray-tracing algorithm (RAYINVR; Zelt and Smith 1992) to enable a 2-D interpretation regarding the geometry of the subducted slab, particularly to determine whether a dipping structure was compatible with the picked traveltime data. The upper mantle phase, PM1P, was identified and picked at offsets larger than 200 km on six shot gathers (e.g., Fig. 2.3). Average uncertainties associated with the picks were calculated by comparing the energy in a limited time window before and after the pick. The ratio of the two energy values is taken as the signal-to-noise ratio (SNR). A traveltime uncertainty, based on the varying SNRs’ values, is then assigned to each pick (e.g. Table 2.2). Because RAYINVR is a 2-D algorithm and Line 1 has a crooked geometry that involves 3-D aspects, we used 2 shots, for which the geophones and the corresponding shots were aligned approximately along a straight line oriented ENE-WSW, in our raytracing analysis. The starting 2-D velocity model was based on the model generated by Fernandez Viejo and Clowes (2003). The synthetic traveltimes were then calculated and compared to the observed traveltimes and the subcrustal velocity model was altered accordingly at each iteration step. This process was continued until we achieved  25  Figure 2.6. Velocity model (numbers in km/s inside the model) and ray path diagram for phases PmP, Pn, and PM1P for shot points 1101 and 1110 and comparison of picked (vertical bars) and calculated (solid lines) traveltimes for seismic phases. (a) Black solid line represents SONRCLE Line 1 and dashed line represents the profile onto which WAR data were projected to interpret the subducted slab. Stars represent shot locations. YK - Yellowknife, NB - Nahanni Butte. In (b) and (c) the PM1P phase represents a wideangle reflection from a subducted slab.  26  acceptable root-mean-square residuals ( Trms ) between the observed and calculated traveltimes for PM1P (Table 2.2). Other phases (PmP, Pn) were also picked and modeled in order to achieve good control over the velocities of the lower crust and upper mantle. Uncertainty tests (Zelt and Smith 1992) for selected nodes in our velocity model give velocity uncertainties of 0.1, 0.15 and 0.15 km/s for the upper crust, lower crust and upper mantle, respectively. Shot Point 1101 1110  Phase PmP Pn PM1P PmP Pn PM1P  Pick Uncertainty (s) 0.12 0.15 0.12 0.12 0.14 0.13  Trms (s) 0.04 0.07 0.05 0.04 0.03 0.07  Table 2.2. Statistics of model fitting for each phase from shots 1101 and 1110. Average uncertainties associated with the picks and the root-mean-square traveltime residual (Trms) with respect to the observed data are included.  2.4.3 Two-dimensional models On shot 1101, the PM1P phase was identified and picked between 240-310 and 375-405 km offset (Fig. 2.3). A good match between the observed and calculated traveltimes ( Trms = 0.05 s) was achieved using a 2-D model in which the upper mantle reflector is represented by an east-dipping interface between 32 km and 90 km depth and with an apparent dip angle of ~8° (Fig. 2.6a). The slope of the calculated traveltime curve is similar to that of the picked traveltimes and is within the estimated traveltime uncertainties. The PM1P phase on shot 1110 is very distinct between offsets of 370 and 440 km. Using the same geometry for the deep reflector as that used for shot 1101 resulted in a slightly larger traveltime misfit ( Trms = 0.07 s). In this case, the slope of the calculated traveltime curve is slightly less than that of the picked traveltimes (Fig. 2.6b), although generally within the traveltime uncertainties. Changing the angle of dip and/or the velocities could not improve the fit for shot 1110 while maintaining that for shot 1101.  27  Projection of the 3-D shot-receiver geometry onto a 2-D profile might explain some of the difficulties in obtaining better fits for shots from opposite directions. Comparing the simple 2-D ray-tracing velocity model with the subducted slab imaged on the NVI data (Fig. 2.2), the results strongly suggest that the PM1P phase is a wideangle reflection from the slab, which represents a velocity contrast in the upper mantle. By projecting the model, which is oriented WSW-ENE, onto the east-west direction (dip direction), we obtain a value of ~13° for the dip angle. This interpretation is consistent with the results from the teleseismic data. 2.4.4 Three-dimensional forward algorithm In order to reproduce the seismic signature of the upper mantle reflections, we used a 3-D forward modeling algorithm (Eaton and Clarke 2000). The program generates synthetic seismograms using the Kirchhoff integral formula, which allows for the computation of the elastic wavefield parameters at any point in a medium based on the wavefield and its normal derivatives on a closed control surface of arbitrary shape. The algorithm can solve for only 2 layers at a time. Thus, the top layer in our velocity model represents the crust and that part of the upper mantle that is above the subducted slab, whereas the bottom layer represents the slab itself. For the top layer, an average velocity value ( v1 = 7400 m/s) was used based on the previous 2-D velocity model. The bottom layer velocity was assigned a value v 2 = 8400 ± 420 m/s (±5% of the upper mantle velocity). The corresponding densities are ρ1 = 2900 kg/m3 and ρ 2 = 3400 kg/m3. A number of interfaces with variable strikes and dips representing the subducted slab were tested and the resulting synthetic seismograms were then compared with the observed seismic profile. The algorithm is relatively fast in terms of computation time; however, it does not account for laterally varying velocities or for vertical velocity gradients. 2.4.5 Three-dimensional models from forward algorithm The regional orientation of the main tectonic domains (Hottah and Fort Simpson) crossed by Line 1 is north-south. They were accreted onto the western margin of North America by a series of east-dipping subduction zones. Hence, we modeled the subducted slab by a 28  number of dipping planes with their strikes ranging between N20°W and N20°E. Only the middle section of Line 1, where the road changes its direction from west-east to almost north-south, was examined (Fig. 2.7a). Synthetic data generated by dipping planes with northwest strikes did not match the observed data either in terms of geometry or traveltimes. Thus, models representing the subducted slab were limited to north/northeast-striking planes only. NVI data show an east-dipping upper mantle reflector, interpreted as the top of the slab, west of the highway junction where Line 1 is oriented west-east (Figs. 2.7a, 2.7b). The reflector reaches a depth of ~22 s and then appears to rise to ~20 s before it becomes flat. East of the junction, where the line changes orientation to NNE-SSW, this reflector extends subhorizontally over ~80 km distance. Toward the northern end of the line segment, the reflector starts to dip eastward again. Two models were judged to provide synthetic reflections that best matched the observed data. In the first one, the slab was represented by an east-dipping plane that flattens beneath the north-south section of Line 1 and then continues dipping eastward (Fig. 2.7c). The slope of the synthetic reflector generated by this model west of the highway junction is lower than that of the real data. The dipping surface reaches a depth of only ~20 s before becoming flat (Fig. 2.7c). East of the junction, the geometry of the synthetic data is very similar to that of the real data, where the reflector is flat over ~80 km and then dips eastward. However, in terms of traveltimes, the synthetic data east of the junction arrive about 0.5 s later than the observed data, whereas west of the junction they arrive ~1.5 s earlier. Some of the discrepancy will be due to the known lateral and vertical variations in velocity (Fernandez Viejo and Clowes 2003) that could not be accounted for with the 3-D modeling algorithm. In the second model, the slab is represented by a continuously dipping plane toward the east and southeast (Fig. 2.7d). The synthetic data are in a very good agreement with the real data (Fig. 2.7b). The geometry is consistent with the NVI image, including the “dip” near the highway junction, a subhorizontal segment and dipping structure near the end of the section. Even with the limitations of one average velocity above the slab, the traveltimes are remarkably consistent with those of the data section. Model 2 29  Figure 2.7. Comparison between the observed and synthetic data from 3-D forward modeling. (a) Thick black line represents the middle part of Line 1 along which the synthetic models are calculated. YK - Yellowknife, NB - Nahanni Butte. (b) The observed data along the middle part of SNORCLE Line 1. Solid black line indicates the top of the slab. (c) The reflecting surface is modeled as an east dipping slab that is flattening at the point where the seismic line changes its course from E-W to NNW-SSE and corresponding synthetic seismogram on the right. (d) The reflecting surface is modeled as a continuously east/southeast dipping slab and corresponding synthetic seismogram on the right. Thick black lines as in (b).  30  demonstrates that recording along a crooked line above a continuously dipping slab can generate a complex reflection geometry. The results from the second model indicate that the two sets of reflectors observed on either side of the highway junction are continuous and could represent a continuously dipping subducted slab. This is our preferred interpretation from the forward modeling. The upward trend of the NVI reflections immediately east of the junction and the subsequent flattening prior to continuation of the eastward dip are artifacts resulting from the projection of a 3-D structure recorded along a crooked line onto a 2-D cross-section. 2.4.6 Three-dimensional finite-difference traveltime inversion To further investigate the 3-D expression of the slab and to determine whether it is flattening or continuously dipping, we implemented a 3-D finite-difference traveltime inversion algorithm that can handle 3-D velocity models and also solve for both NVI and R/WAR data. The algorithm was introduced first by Vidale (1988) as a method of computing traveltimes using a finite-difference approximation to the eikonal equation to propagate first-arrival traveltimes in 2-D and 3-D velocity models. Hole and Zelt (1995) modified the algorithm to correctly handle large, sharp contrasts in velocity, which enabled solutions for complex 3-D reflector geometries. In the modified version, the reflector model is allowed to vary smoothly in depth, thereby increasing the accuracy compared to a discretized reflector model. The velocity model is based on the one we obtained from ray-tracing (Fig. 2.6), which showed that lateral velocity changes in the lower crust and upper mantle beneath the roughly north-south segment of Line 1 (Fig. 2.6) are small and within the velocity uncertainties established by ray-tracing (0.15 km/s). Hence, we opted to use a 2.5-D velocity model in which velocities are changing vertically and in the east-west direction but constant in the north-south direction. The dimensions of the gridded velocity model are n x × n y × n z = 276 × 106 × 53 (i.e., ~1.5 million nodes) or 550 × 210 ×104 km with a grid spacing h = 2 km.  31  The calculated traveltimes for a reflecting surface specified in a starting model are compared to the observed traveltimes. The time differences (dt) between the calculated and observed traveltimes are related to small changes in reflector depth (dz) through the dip angle α of the reflector and the velocity above the reflector. Once calculated, (dz) is then added to the reflecting surface, which becomes the new starting surface. The process is repeated until (dz) is within, or close to, the value of the grid spacing. 2.4.7 Three-dimensional models from traveltime inversion We applied the inversion algorithm to the picked traveltimes for both the WAR and NVI data associated with the mantle reflector. In the case of the WAR data, a total of 6 shots (1101, 1103, 1104, 1107, 1109 and 1110), where PM1P is observed, were used (Fig. 2.8a). As a starting model for the subducted slab, we chose a horizontal plane at 70 km depth. The results showed an east- dipping slab in agreement with the reflection data. However, after 3 iterations the rms residuals ( Trms ) and the associated average difference in depth (dz) stopped decreasing and were still too high (Table 2.3). Changing the depth of the starting model (between 50 and 90 km) did not improve the error. Thus, the horizontal starting model was substituted with a plane dipping at 13° and an average depth of 60 km. The inversion results showed a slab dipping toward the east between ~80 and ~90 km depth and at an angle of 11° (Fig. 2.8b). These values are consistent with the bottom of the slab as interpreted on the NVI data and not with the top of the slab (Fig. 2.5c). Both Trms and dz values were reduced significantly after the 1st iteration (Table 2.3). For the NVI data, three sets of traveltimes were picked, one from the top of the slab and two from its bottom, to incorporate the alternative interpretations shown in Fig. 2.5c. The starting depth model is a dipping plane identical to the previous model used for the WAR data. The resulting slabs from the three sets of inversions are similar in terms of their geometry and angle of dip. Only the absolute depths are different. Figure 2.8c shows the calculated slab depths for the bottom of the slab as shown by our preferred interpretation A – B in Fig. 2.5c. The top of the slab is about 10 km shallower, consistent with the ~2 s two-way travel time difference between the top and bottom of the slab. 32  Figure 2.8. Results for 3-D finite difference inversion for the WAR and NVI data. (a) The colored rectangle represents the part of the slab recovered by the WAR and NVI inversions. White box represents the part of the slab recovered from the WAR inversion only. Red stars represent the location of the shot points used for the inversion of the wide-angle data. YK - Yellowknife, NB - Nahanni Butte. (b) Plan view of the recovered slab from the inversion of the wide-angle data only. Black line is the seismic line and white dots are the midpoints between the shots and the receivers where the depth to the slab was calculated. (c) Plan view of the lower surface of the slab recovered from the NVI data. Note that the two surfaces are not identical since two different datasets (WAR and NVI) were used for the inversion. However, the difference is within the depth uncertainty (2-3 km).  The final step was solving for both the WAR and NVI data simultaneously. First, the WAR data were solved with the NVI picks from the top of the slab. Then the combined data were solved with the NVI picks from the bottom interface, including the two alternative interpretations shown in Fig. 2.5c. For all inversions, the rms traveltime error, Trms , was reduced significantly within the first 3 iterations and did not improve by more than 0.01 s after the 3rd iteration (Table 2.3; Fig. 2.9). On the other hand, the associated average difference in depth (dz) after the 3rd iteration was smaller than the grid 33  spacing; thus no further iterations were attempted. In all cases, the results show a continuously east-dipping slab between 30 and 95 km with an average dip angle of ~16° (Fig. 2.10). The region of the surface where there is very good depth control is limited to the location of the midpoints, whereas the depth to the remainder of the surface is determined by the interpolation part of the algorithm. The 2.5-D velocity model has the same velocity uncertainties as the 2-D velocity model (section 4.2) and a grid spacing h =  Data type WAR Data (Horizontal plane) WAR Data (Dipping plane) NVI Data (Top of Slab) NVI Data (Bottom of Slab A & C) NVI Data (Bottom of Slab A & B) WAR & NVI Data (Top of Slab) WAR & NVI Data (Bottom of Slab A & C) WAR & NVI Data (Bottom of Slab A & B)  1st iteration dz (km) T rms (s) 1.96 13.5 1.66 -10.71 3.07 11.17 4.14 15.31 4.02 14.87 1.01 12.38 1.32 12.1 1.16 11.12  3rd iteration dz (km) T rm s (s) 0.84 2.91 0.27 1.74 0.46 0.38 0.52 -0.81 0.5 -0.73 0.88 1.83 0.61 -1.12 0.36 0.57  Table 2.3. Statistics of model fitting resulting from the finite-difference traveltime inversion algorithm. The left column represents the type of data used by the algorithm to generate the final reflector model based on a starting model. Trms is the root-mean-square traveltime residual and dz is the average difference in depth of the starting model and the recovered model.  2 km. This would result in a depth uncertainty in the modeled slab of 2 – 3 km. The modeled slab resulting from the WAR data and the top NVI picks exhibits a pronounced trough around the centre of the slab (Fig. 2.10b) suggesting either a deeper surface at that location or anomalously low upper mantle velocities. The two models (Figs. 2.10c, 2.10d) resulting from the WAR data and the bottom NVI picks (with the two alternative interpretations) appear more realistic. However, the model with NVI picks along reflections A and B (Fig. 2.5c) has the smoothest surface, the simplest structure (Fig. 2.10d) and the smallest rms traveltime residuals (Table 2.3). Therefore it is our preferred model. Notice that the rms residuals of the 2-D ray tracing (Table 2.2) are much smaller than the rms residuals of the 3-D finite-difference traveltime inversion (Table 2.3).There are two reasons for this: (1) for waves traveling in the upper mantle at 8.4 34  km/s, a 2-km grid spacing used in the finite-difference inversion is equal to ~0.25 s, representing the minimum traveltime misfit that can be achieved; (2) the 2-D forward ray-tracing does not have such limitations. The analytic solution in forward ray-tracing allows for misfits less than traveltime uncertainties by selection of suitable models. However, many such models can be generated and may be geologically unrealistic.  Figure 2.9. Comparison of picked (observed) and calculated traveltimes for WAR and the NVI data from the finite-difference inversion. (a) WAR and NVI data picked from the top of the slab (b) WAR and NVI data picked from bottom of the slab (reflections A – B, Fig. 2.5c). Note the improvement of the data misfit between the 1st and 3rd iteration.  2.5 Discussion 2.5.1 Slab geometry The results from the forward and inverse modeling of both the WAR data and the NVI data strongly support our hypothesis that the subhorizontal reflectors imaged by the NVI data beneath the Wopmay orogen are part of the previously interpreted Paleoproterozoic subducted slab (Fig. 2.11), instead of being a separate entity as suggested by Cook et al. (1999). Below the near-surface position of the Fort Simpson terrane, the two parallel reflectors at Moho level (~32 km depth) represent the underthrust lower crust of the Fort Simpson terrane. Downdip, the other parallel pair of dipping reflectors at approximately 65 and 75 km depth, together with the subhorizontal pair of reflectors further east, represent the top and bottom of a fossil subducted oceanic crust. 35  Figure 2.10. 3-D perspective view of the subducted slab resulting from the inversion of both the wide-angle and the NVI data. (a) The black line indicates the seismic line and the purple dots indicate the locations at which the depth to the slab has been calculated for both the NVI and WAR data. Away from the midpoints, the depth of the slab is calculated by interpolation. (b) The results from inverting the NVI data picked from the top of the slab and the wide-angle data. (c) The results from inverting the wide-angle data and the NVI data picked from the bottom of the slab along reflections A and C (Fig. 2.5c). (d) The results from inverting the wide-angle data and the NVI data picked from the bottom of the slab along reflections A and B (Fig. 2.5c).  The slab extends eastward with relatively constant thickness of ~10 km for another 200 km beneath the Great Bear magmatic arc before it disappears once it reaches depths of 90 - 100 km. This can be explained by slab delamination in which any remainder of the slab has sunk into the mantle, or by reaching a level at which metamorphic minerals were converted to higher grade and became seismically indiscernible from mantle rocks. In modern subduction zones, the depth range of 75-85 km marks the lower stability limit of 36  amphibole minerals in mafic rocks, suggesting a region of active eclogitization and fluid release (ANCORP working group 2003).  Figure 2.11. The re-interpreted seismic reflection profile in the upper mantle along SNORCLE Line 1. The fossil subducted slab is continuously dipping between ~32 km and ~100 km depth over a profile distance of ~400 km (lateral distance of ~300 km). Thin solid line represents the top of the slab imaged by the NVI data only. Thick solid line represents the bottom of the slab imaged by both the NVI and WAR data. Here we have drawn our preferred interpretation following interface A - B in Fig. 2.5c. Thick dashed line delineates the top of an anisotropic layer representing a lid of the uppermost mantle material inferred from the teleseismic data. Thin dashed lines highlight deeper reflectivity that may represent an upper mantle petrological boundary.  Based on the results from the 3-D forward modeling (Fig. 2.7) and traveltime inversion (Fig. 2.10), the slab seems to be continuously dipping but with a varying angle. Figure 2.12 shows our preferred 3-D impression of the descending slab based on the joint solution of both the WAR and the NVI data. The slab has a dip angle of ~14° between 33 km and 62 km depth over a distance of ~60 km before it starts dipping shallowly at an angle of ~7°. The lateral and vertical extent of the slab depicted in Fig. 2.12 is speculative since the farthest midpoints at which the depth to the slab is calculated are only 50 km away from the seismic line. We argue that the apparent flattening of the slab as interpreted on Fig. 2.11 is an artifact due to the projection of the 3-D seismic line onto a 2-D cross section. The beginning of subduction of this oceanic plate probably started after the final accretion of the Hottah terrane around 1.88 Ga and continued until around 1.84 Ga,  37  which marked the last recorded magmatic activity in the Great Bear magmatic arc (Gandhi et al. 2001).  Figure 2.12. Schematic illustration of the subducted slab (solid lines). Long dashed line on the surface of the subducted oceanic crust represents the downdip extent beyond which the depth of the slab is unconstrained (short dashed lines).  2.5.2 Origin of the seismic signature of the slab The modeled slab generated by the traveltime inversion of the WAR data (Fig. 2.8b) corresponds to the bottom of the subducted oceanic crust imaged by the reflection data (Fig. 2.2). This indicates that the strong PM1P phase is being reflected from the lower boundary of the oceanic crust. Such a result can readily be explained if the layer underlying the subducted oceanic crust has a relatively higher velocity. Results from the teleseismic data (Mercier et al. 2008) reveal a 10- to 15-km-thick layer which has average mantle velocities and anisotropy decreasing with depth, and which parallels the bottom of the subducted oceanic crust (Fig. 2.4). The fast axis lies in the same direction as that of  38  the PM1P phase propagation. This anisotropic mantle lid represents the higher velocity layer needed for the strong PM1P reflections. On the other hand, the lack of wide-angle reflections from the top of the subducted crust indicates that it has an impedance (velocity x density) comparable to that of the surrounding mantle (Hansen and Balling 2004). Assuming this is true, we indicate a scenario by which this could be established. As the subducted plate descended, it underwent progressive metamorphism (Stern 2002). The hydrous basalt and gabbro converted to blueschist and then amphibolite. At depths greater than 40 km, the subducted oceanic crust entered the stability field of eclogite converting amphibolite into the latter (Peacock and Wang 1999). Eclogite has velocities and densities that are comparable to those of peridotite, the major constituent of the surrounding mantle (Christensen 1996; Ji et al. 2002). This similarity reduces the impedance contrast and hence explains the absence of PM1P wide-angle reflections from the top of the subducted oceanic crust. Unfortunately, mineralogic profiles of a downgoing plate are difficult to construct because of the large uncertainties regarding the thermal structure of subduction zones and the complex mineralogical changes associated with the wide range of subducted components. Also, the validity of equilibrium models especially for cold material being slowly heated and squeezed is not clear, because phase transformations are likely to be kinetically retarded (Stern 2002). From NVI seismic reflection data recorded in the Baltic Shield and North Sea, 1.9 to 1.8 Ga old fossil subduction zones preserved in the upper mantle also have been interpreted (BABEL Working Group 1990; Balling 2000; Hansen and Balling 2004). Results from these studies suggest that the transformation of basaltic/gabbroic material to eclogite within the subducted slab is likely the main source for the upper mantle reflectivity that was observed. Along SNORCLE Line 1, all three datasets (WAR, NVI and teleseismic) detected a well-defined lower boundary between the subducted oceanic crust and the uppermost mantle. However, only the NVI data were able to detect the top of the oceanic crust. Laboratory petrophysical measurements, field studies and synthetic reflection modeling  39  have attributed prominent reflectivity to the seismic impedance contrasts associated with shear zones, mainly resulting from mylonites (e.g. Christensen and Szymanski 1988; Wang et al. 1989; Shaocheng et al. 1997). A few hundred meters thick shear zone immediately above the subducted oceanic slab is probably too thin to be detected by the low frequency wide-angle and teleseismic data (vertical resolution is around 2 km) but sufficient enough to produce strong reflections when probed by the NVI data (vertical resolution is ~80 m). Hansen and Balling (2004) investigated the seismic signatures (WAR and NVI data) of a number of upper-mantle reflectors through analytical and numerical modeling, including full wavefield modeling. They found that both uppermantle shear zones and relict subduction slabs comprising intermediate-velocity eclogite may contain sufficient impedance contrast to produce near-vertical incidence reflectivity but not enough to generate wide-angle reflections. This is in good agreement with our interpretation. 2.5.3 Flat versus low-angle subduction A large number of parameters that control the dip of a slab have been suggested: age of the downgoing plate, absolute plate motion, slab pull, or subduction polarity. Slab pull has been widely accepted as the primary influence (Vlaar and Wortel 1976; Molnar and Atwater 1978). However, recent statistical analyses, which examined these parameters along 159 oceanic subduction transects, found that a global correlation exists between the overriding plate absolute motion, slab dip and back-arc deformation (Lallemand et al. 2005; Heuret and Lallemand 2005). Shallow slab dips (< 50°) are associated with backarc shortening where the overriding plate is moving trenchwards. This can be readily understood since the dip angle has an effect on the degree of interplate coupling through the area of friction between converging plates, which increases when the angle of dip decreases (Gutscher and Peacock 2003). Unfortunately, all these observations are only for active subduction zones since any reconstruction of slab dip through time is very speculative (Lallemand et al. 2005). Although the low-angle of the descending slab imaged beneath the Hottah terrane and the Great Bear arc raises the possibility of a flat subduction setting, the main  40  characteristics of such a subduction zone are missing. Gutscher et al. (2000) examined 10 current regions where flat subduction has been identified (Peru, Central Chile, Ecuador, Northwestern Colombia, Costa Rica, Mexico, Cascadia, Southern Alaska, SW Japan and Western New Guinea). The one element they found involved in nearly all cases was the presence of anomalously thick (15 - 20 km) oceanic crust with its buoyancy effect being the source for flat subduction. On the other hand, the subducted oceanic crust beneath the Hottah terrane is only 8 - 10 km thick as imaged by the NVI data. Also, most of these flat subduction zones did not develop a magmatic arc and in the case where the magmatic arc is present, the trench – arc gap is in excess of 400 km. Comparing it to the Hottah subduction, the distance between the Great Bear magmatic arc and the trench is ~ 200 km, which is near the upper limit for the trench – arc gap (166 ± 60 km) in a normal subduction zone (Gill 1981). Taking into consideration the shortening associated with shallow dip angle, the original distance between the Great Bear magmatic arc and the trench could have been larger. However, a decrease from an original 400 km arc-trench distance, characteristic for a flat subduction, to the present 200-km gap requires 200 km of shortening. Such a substantial shortening is only observed where two continental plates are colliding such as in the Himalaya (Schelling 1992) or the Alps (Menard et al. 1991). Thus, although the original arc-trench distance may have been larger, it is unlikely that it was large enough to be characteristic of a flat subduction. Another element that is common among many flat subduction zones is the presence of anomalous adakitic geochemical signatures found in relic or present-day arc magmas (Defant and Drummond 1993; Morris 1995). No adakitic magmatism was found in the Great Bear magmatic arc (Gandhi et al. 2001). Thus, there is no compelling evidence that the relict subduction below Hottah terrane is a flat subduction. It would more appropriately be regarded as low-angle or shallow subduction, similar to other relict subduction zones interpreted from NVI reflection data in the Archean Superior Province (Calvert et al. 1995) and Paleoproterozoic Baltic Shield and North Sea region (BABEL Working Group 1990; Balling 2000; Hansen and Balling 2004). These relict subduction zones provide direct evidence that plate tectonics has been active since at least 2.6 Ga. All are characterized by low dip angles ranging between 15˚and 25˚. Such observations attest to a proposal that  41  shallow/flat subduction may have been more common in Earth’s early history. Other examples of shallow subduction zones are found in southwest Japan and Cascadia in western North America, but these are much younger and magmatic arcs have been generated. Assuming that the relict slab in the Wopmay orogen is the result of shallow subduction, what are some of the implications that could be associated with such a feature? Vlaar (1986) and Helmstaedt and Schulze (1989) argue that Archean cratons were stabilized through shallow subduction processes. Helmstaedt and Schulze (1989) further suggest that subcrustal lithosphere is generated and stabilized by stacking of subducted slabs. In considering three mechanisms by which the ancient continental lithosphere could be formed and stabilized, Abbott (1991) argues that basal addition of oceanic slabs during buoyant subduction processes is the most applicable one. For the Archean, she argues that a hotter mantle would simultaneously cause the production of more thick oceanic crust and more young oceanic plates, relative to the present. Hotter mantle temperatures also would increase the depletion of the oceanic lithospheric mantle making it more buoyant. As a result, flat or shallow subduction of buoyant slabs would have been more common than at present and the stacking of these would increase the thickness of the lithosphere and its preservation. Furthermore, buoyant subduction of cooler oceanic lithosphere could allow the continental lithosphere to thicken more rapidly because it would protect the base of the lithosphere from thermal erosion by the hotter Archean mantle. Most of these arguments are made in relation to generation of continental lithosphere beneath Archean regions in which diamonds are found. To a certain extent, similar arguments could be made about the effects of shallow subduction during Paleoproterozoic tectonics, such as we interpret for the Wopmay orogen, although the overall effects might be lessened somewhat due to cooler mantle temperatures relative to the Archean, but ones still higher than at present. With this perspective at hand, we suggest that our interpreted shallow subducted slab that extends laterally for about 300 km represents an addition to, and thickening of, the Paleoproterozoic sub-crustal mantle and may have contributed to stabilization of the lithosphere beneath the Wopmay orogen.  42  2.5.4 Relation with other upper-mantle reflectors The interpreted ~10-km-thick anisotropic layer immediately below the bottom of the subducted oceanic crust (Mercier et al. 2008) is characterized by a sharp pulse followed by a more diffuse one (Fig. 2.4). This signature is similar to that of another upper-mantle feature at 70-80 km depth beneath the Slave craton from which teleseismic converted S waves were identified by Bostock (1998). Mercier et al. (2008) explained this similarity as being due to a common signature of fossil subduction recorded in two distinct episodes of plate convergence, rather than an indication of the same feature. Unfortunately, the lack of upper mantle reflectivity on the NVI data beneath the Hottah-Slave boundary makes it difficult to establish a geometric relation between the subducted slab below the Great Bear magmatic arc and the layered structure beneath the Slave craton. Reprocessing of the remainder of the NVI data beneath the Slave craton to identify any possible weak reflectors would be needed in order to investigate any geometric relation between the relict slab and the anisotropic layer imaged by the teleseismic data. The deep band of reflectivity observed beneath the slab at 25 – 30 s (Figs. 2.5b, 2.5c) correlates well with an upper mantle reflector identified by the previous wide-angle reflection study (Fernandez Viejo and Clowes 2003). Because there is no obvious relation between these reflections and other features at or near the surface, it is difficult to associate them with any tectonic events. Two possible explanations come from mantle xenolith studies in the Slave craton and the Kaapvaal craton (Kopylova et al. 1998; Helmstaedt and Schulze 1989). Kopylova et al. (1998) showed that the transition from spinel-garnet peridotite to garnet-only peridotite occurs at around 120 km depth and temperatures of about 800°C. Thus, these reflections may represent a petrological boundary in the upper mantle indicating a transition to a garnet-only peridotite. Alternatively they may represent a subcrustal imbricate structure or tectonic wedging as proposed by Helmstaedt and Schulze (1989).  2.6 Conclusions Our results from 2-D and 3-D forward and inverse modeling of WAR and NVI data suggest that a parallel pair of reflectors in the upper mantle as observed on NVI data 43  along SNORCLE Line 1 are an extension of the Paleoproterozoic slab interpreted by an earlier study (Cook et al. 1999), rather than a separate feature. The apparent flattening of the slab as imaged by the NVI data is not real; rather, it is an artifact resulting from the projection of a 3-D seismic line geometry onto a 2-D cross section. The relict slab is the result of shallowly east-dipping subduction of oceanic lithosphere, which once separated the Fort Simpson and Hottah terranes, beneath the Hottah terrane. It extends from the Fort Simpson – Hottah boundary to the Hottah – Slave boundary, where images of it disappear, a lateral distance of ~300 km. Its depth ranges from 35 to 100 km. The subduction zone was active between 1.88-1.84 Ga and is responsible for the formation of the Great Bear magmatic arc. Refraction/wide-angle reflection, near-vertical incidence and teleseismic data reveal a distinct interface between the base of the subducted oceanic crust and the uppermost mantle. However, the top of the slab is imaged only by the NVI data, suggesting a finer-scale impedance contrast. This relict subduction zone does not exhibit any of the main characteristics of modern flat subduction zones but it does have a low-angle, continuously dipping slab similar to SW Japan and Cascadia subduction zones. Our revised interpretation of a laterally extensive, preserved, and shallowly subducted slab beneath a ca. 1.8 Ga orogenic belt indicates that such subduction may play a role in the generation and stabilization of new continental lithosphere, as inferred by others (e.g., Abbott 1991; Helmstaedt and Schulze 1989; Mercier et al. 2008).  44  Bibliography Abbott, D. 1991. The case for accretion of the tectosphere by buoyant subduction. Geophysical Research Letters, 18, p. 585– 588. Abbott, G., Thorkelson, D., Creaser, R., Bevier, M.L., and Mortensen, J. 1997. New correlations among Proterozoic successions and intrusive breccias in the Ogilvie and Wernecke Mountains, Yukon. In SNORCLE Transect and Cordilleran Tectonics Workshop Meeting, University of Calgary, Calgary, Alta., March 7-9, 1997. Compiled by F.A. Cook, and P. Erdmer. LITHOPROBE Secretariat, The University of British Columbia, Vancouver, B.C. Lithoprobe Report No. 56, p. 188–197. Aitken, J.D. 1993. Cambrian and lower Ordovician – Sauk sequence, In Sedimentary Cover of the Craton in Canada. Geological Survey of Canada, Geology of Canada. Compiled by D.F. Stott and J.D. Aitken. 5, p. 96-124. ANCORP Working Group, 2003. Seismic imaging of a convergent continental margin and plateau in the central Andes (Andean Continental Research Project 1996 (ANCORP’96)). Journal of Geophysical Research, 108, 2328, doi: 10.1029/2002JB001771. BABEL Working Group, 1990. Evidence for early Proterozoic plate tectonics from seismic reflection profiles in the Baltic Shield. Nature, 348, p. 34–38. Balling, N. 2000. Deep seismic reflection evidence for ancient subduction and collision zones within the continental lithosphere of north-western Europe. Tectonophysics, 329, p. 259–291. Bostock, M.G. 1998. Mantle stratigraphy and evolution of the Slave province. Journal of Geophysical Research, 103, p. 183-200. Calvert, A.J., Sawyer, E. W., Davis, W. J., and Ludden, J. N. 1995. Archean subduction inferred from seismic images of a mantle suture in the Superior Province. Nature, 375, p. 670-674. Christensen, N. I. 1996. Poisson’s ratio and crustal seismology. Journal of Geophysical Research, 101, p. 3139-3156. Christensen, N. I., and Szymanski, D. L. 1988. Origin of reflections from the Brevard Fault Zone. Journal of Geophysical Research, 93, p. 1087-1102.  45  Coles, R.L., Haines, G.V., and Hannaford, W. 1976. Large-scale magnetic anomalies over western Canada and the Arctic: a discussion, Canadian Journal of Earth Sciences, 13, p. 790– 802. Cook, F. A., and Erdmer, P. 2005. Introduction to special issue on the LITHOPROBE Slavenorthern Cordillera lithospheric evolution (SNORCLE) transect. Canadian Journal of Earth Sciences, 42, p. 865-867. Cook, F.A., van der Velden, A.J., Hall, K.W., and Roberts, B. 1998. Tectonic delamination and subcrustal imbrication of the Precambrian lithosphere in northwestern Canada mapped by LITHOPROBE. Geology, 26, p. 839-842. Cook, F.A., van der Velden, A.J., and Hall, K.W. 1999. Frozen subduction in Canada's Northwest Territories: LITHOPROBE deep lithospheric reflection profiling of the western Canadian Shield. Tectonics, 18, p. 1-24. Defant, M. J., and Drummond, M. S. 1990. Derivation of some modern arc magmas by melting of young subducted lithosphere. Nature, 347. p. 662-665. Eaton, D. W., and Clarke, G. 2000. A Kirchhoff integral method to model 3-D elastic-wave scattering. In SEG Technical Program Expanded Abstracts, p. 2472-2475. Fernandez Viejo, G., and Clowes, R.M. 2003. Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay orogen, northwestern Canada, from a LITHOPROBE refraction/wide-angle reflection survey. Geophysical Journal International, 153, p. 1-19. Fernandez Viejo, G., Clowes, R.M., and Amor, J.R. 1999. Imaging the lithospheric mantle in northwestern Canada with seismic wide-angle reflections. Geophysical Research Letters, 26, p. 2809-2812. Gandhi, S.S., Mortensen, J.K., Prasad, N., and van Breemen, O. 2001. Magmatic evolution of the southern Great Bear continental arc, northwestern Canadian Shield: geochronological constraints. Canadian Journal of Earth Sciences, 38, p. 767-785. Gill, J. B. 1981. Orogenic Andesites and Plate Tectonics. Springer-Verlag, New York Gutscher, M., and Peacock, S. M. 2003. Thermal models of flat subduction and the rupture zone of great subduction earthquakes. Journal of Geophysical Research,108, doi:10.1029/2001JB000787.  46  Gutscher, M., Spakman, W., Bijwaard, H., and Engdahl, E.R. 2000. Geodynamics of flat subduction: Seismicity and tomographic constraints from the Andean margin. Tectonics, 19, p. 814-833. Hansen, T. M., and Balling, N. 2004. Upper-mantle reflectors: modeling of seismic wavefield characteristics and tectonic implications. Geophysical Journal International, 157, p. 664-682. Heuret, A., and Lallemand, S. 2005. Plate motions, slab dynamics and back-arc Deformation. Physics of the Earth and Planetary Interior, 149, p. 31-51. Helmstaedt, H., and Schulze, D. J. 1989. Southern African kimberlites and their mantle sample: Implications for Archean tectonics and lithosphere evolution. In Kimberlites and Related Rocks, Blackwell , Cambridge, Mass., p. 358– 368. Hildebrand, R. S., and Bowring, S. A. 1999. Crustal recycling by slab failure. Geology, 27, p. 1114. Hildebrand, R. S., Hoffman, P.F., and Bowring, S.A. 1987. Tectono-magmatic evolution of the 1.9-Ga Great Bear magmatic zone, Wopmay orogen, northwestern Canada. Journal of Volcanology and Geothermal Research, 32, p. 99-118. Hoffman, P.F. 1989. Precambrian geology and tectonic history of North America. In Bally, A. W., and Palmer, A. R., eds., The Geology of North America: An Overview. Geological Society of America, Boulder , Colorado, v. A, p. 447-512. Hoffman, P. F., and Bowring, S. A. 1984. Short-lived 1.9 Ga continental margin and its Destruction. Wopmay orogen, NW Canada. Geology, 12, p. 68-72. Hole, J.A., and Zelt, B.C. 1995. 3-D finite-difference reflection traveltimes. Geophysical Journal International, 121, p. 427-434. Ji, S., Wang, Q., and Xia, B. 2002. Handbook of Seismic Properties of Minerals, Rocks and Ores. École Polytechnique de Montréal, Montreal, 630 p. Jones, A. G., Lezaeta, P., Ferguson, I. J., Grant, N., McNeice, G., Spratt, J., Farquharson, C., Roberts, B., Wennberg, G., Wolynee, L., and Wu, X. 2005. A 1600 km-long magnetotelluric transect from the Archean to the Tertiary: SNORCLE MT overview. Canadian Journal of Earth Sciences, 42, p. 955-981.  47  Lallemand, S., Heuret, A., and Boutelier, D. 2005. On the relationships between slab dip, backarc stress, upper plate absolute motion, and crustal nature in subduction zones. Geochemistry Geophysics Geosystems, doi:10.1029/2005GC000917. Menard, G., Molnar, P., and Platt, J. 1991. Budget of crustal shortening and subduction of continental crust in the Alps. Tectonics, 10, p. 231-244. Mercier, J.P., Bostock, M.G., Audet, P., Gaherty, J.B., Garnero, E.J., Revenaugh, J. 2008. The teleseismic signature of fossil subduction: Northwestern Canada. Journal Geophysical Research., 113, B04308, doi: 10.1029/2007JB005127. Molnar, P., and Atwater, T. 1978. Interarc spreading and Cordilleran tectonics as alternates related to the age of subducted oceanic lithosphere. Earth Planet Scientific Letters, 41, p. 827857. Morris, P.A. 1995. Slab melting as an explanation of Quaternary volcanism and aseismicity in southeastern Japan. Geology, 23, p. 395-398. Peacock, S. M. 2002. Thermal structure and metamorphic evolution of subducting slabs. In J. Eiler and G. A. Abers, eds., Washington, D.C., American Geophysical Union Geophysical Monographs Series, 138, p. 7-22. Peacock, S., M., and Wang, K. 1999. Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan. Science, 286, p. 937-939. Schelling, D. 1992. The tectonostratigraphy and structure of eastern Nepal Himalaya. Tectonics, 11, p. 925-943. Shaocheng, J., Long, C., Martignole, J., and Salisbury, M. 1997. Seismic reflectivity of a finely layered, granulite-facies ductile shear zone in the southern Grenville Province (Quebec). Tectonophysics, 279, p. 113-133. Stern, R., J. 2002. Subduction zones. Reviews of Geophysics, 40, doi:10.1029/2001RG000108. Thompson, R.I., Mercier, E., and Roots, C.F. 1987. Extension and its influence on Canadian Cordilleran passive margin evolution. In Continental Extension Tectonics, Geological Society of London, Special Publication, p. 409-417. Vidale, J. 1988. Finite difference calculation of travel times. Bulletin of Seismological Society of America, 78, p. 2062-2076.  48  Villieneuve, M.E., Theriault, R.J., and Ross, G.M. 1991. U-Pb ages and Sm-Nd signature of two subsurface granites from the Fort Simpson magnetic high, northwest Canada. Canadian Journal of Earth Sciences, 28, p. 1003-1008. Vlaar, N. J. 1986. Archean global dynamics. Geology, 65, p. 91– 101. Vlaar, N. J., and Wortel, M. J. R. 1976. Lithosphere aging, instability and subduction. Tectonophysics, 32, p. 331-351. Wang, C. Y., Okaya, D. A., Ruppert, C., Davis, G. A., Guo, T. S., Zhong, Z., and Wenk, H. R. 1989. Seismic reflectivity of the Whipple Mountain shear zone in Southern California. Journal of Geophysical Research., 94, p. 2989-3005. Zelt, C.A., and Smith, R.B. 1992. Seismic traveltime inversion for 2-D crustal velocity structure. Geophysical Journal International, 108, p. 16 – 34. Zhao, D. 2001. New advances of seismic tomographic and its applications to subduction zones and earthquake fault zones. Island Arc, 10, p. 68-84.  49  Chapter 3  Enhancing crustal reflection data through curvelet denoising2 3.1 Introduction In near-vertical incidence (NVI) seismic reflection data, recorded wavefronts (i.e., reflections) arise from the interaction of the incident wavefield with inhomogeneities in the Earth’s subsurface. These wavefronts can become contaminated with various types of noise due to acquisition and processing problems (Olhovich, 1964). The noise sources can be separated into two categories. The first category of noise is associated with the experiment itself. Such noise comprises any unexpected perturbation of the recording environment during data acquisition. A geophone may malfunction or the recording systems may have glitches creating erratic noise in the seismic record. Wind motion or cable vibrations can generate random noise. Outside factors, such as the activities of mammals and/or drilling rigs for marine acquisition, might also contaminate seismic records. Such noise sources can create coherent energy in the data that may be misinterpreted as true signal. The second category of seismic noise derives from modeling or processing uncertainties (Tarantola, 1987). Such uncertainties occur when the physical description and parameterization of the earth is incomplete, as it inevitably must be. This incomplete description results from the inherent complexity of wavefield propagation in the subsurface, leading to various noisy components in the data. Noise in seismic data corrupts the quality of the signal; this can lead to poor interpretation or misinterpretation. Thus, separation of signal and noise forms an important part of seismic data processing, particularly for crustal data where the signal-to-noise ratio is low. 2  A version of this chapter has been accepted for publication: Kumar, V., Oueity, J., Clowes, R. M., and  Herrmann, F. J. Enhancing crustal reflection data through curvelet denoising, Tectonophysics, doi: 10.1016/ j.tecto.2010.07.017.  50  Because noise attenuation is a major step in seismic processing, many methods have been developed to suppress this incoherent signal component (e.g., Jones and Levy, 1987; Donoho, 1995; Ulrych et al., 1999; Fomel, 2002; Zhang and Ulrych, 2003). Within the seismic exploration industry, predictive deconvolution is one of the most common techniques to attenuate noise. This adaptive method is also known as F-X (frequencyspace) deconvolution or F-X prediction (Abma and Claerbout, 1995). The prediction filter works in a windowed sense, predicting one dip at a time and tends to attenuate all other secondary dips in that window. Another method for incoherent noise suppression is localized slant stack (McMechan, 1983; Milkereit, 1987). The localized slant stack method (also known as linear radon transform) is based on assumptions of linear events in a localized window. The window size is fixed and all linear events in the window can be modeled. However, the method fails to model curved and complicated events in the window. Thus F-X prediction and slant stack methods are not suited for curved events and for events with conflicting dips such as those that occur at caustics and pinch-outs. Recently, curvelets have been introduced to resolve problems associated with retention of complicated events involving directionality, such as multiple dips at the same timespace location or curved events. Curvelets belong to a family of multiscale and multidirectional data expansions (Candés et al., 2004). They are redundant transforms that preserve energy through both the forward and inverse operations. They are contained within strictly localized anisotropic wedges in Fourier space and decay rapidly in physical space where they act as little plane waves, exhibiting oscillations in one direction and smoothness in perpendicular directions (Fig. 3.1; Candés et al., 2006). Because high frequencies are more directional, the number of wedges increases with frequency, which leads to fine-scale curvelets that are increasingly anisotropic (needlelike) and hence increasingly directional selective. As shown in Fig. 3.1, each curvelet is associated with position, frequency and angle. The curvelets differ from representation in the F-K domain in the sense that curvelets are able to decompose conflicting dips into localized functions having different frequency and wavenumber bands, while F-K functions are global in nature and specific dipping events are not isolated (Fig. 3.1). In seismic data, the coherent components of signal 51  differ from the incoherent noise in terms of structure that can be followed across at least a few traces. For example, dipping events map to the coefficients of specific angles in the curvelet domain whereas noise maps over all angles. Thus the coherent energy leads to significant coefficients that can be separated from the diffused noise by sparsitypromotion (Neelamani et al., 2008; Hennenfent and Herrmann, 2006). Curvelets also have been applied to other seismic processing issues such as ground roll removal (Yarham and Herrmann, 2008 ), coherent noise removal (Neelamani et al., 2008) and multiple separation (Wang et al., 2008).  Figure 3.1. (a) Five curvelets in time-space domain and (b) their representation in the frequencywavenumber domain. Each curvelet function is spatially localized because its amplitude rapidly decays to zero outside a certain region. The curvelets are localized in wedges in the f-k domain. The orientation of a curvelet in the time-space domain is perpendicular to its wedge in the f-k domain. (c) Curvelet (top row) and f-k (bottom row) decomposition of two events with conflicting dips. Note the localization of curvelet functions and global nature of the f-k functions used to decompose the image. (Adapted from Neelamani et al. 2008).  52  Because noise tends to dominate the signal in crustal reflection data, we propose a curvelet-based denoising procedure that is designed to promote sparsity (i.e., to minimize the number of coefficients required) and to preserve coherent energy. First, we introduce our algorithm based on sparsity promotion in the curvelet domain and demonstrate its applicability on noisy synthetic data. This is followed by application of the algorithm to a real shot gather and a stacked section from a crustal reflection experiment in Lithoprobe’s SNORCLE transect in the Northwest Territories, Canada (Cook et al., 1999).  3.2 Curvelet denoising The seismic denoising problem corresponds to estimating the unknown noise-free data, m, from the noisy data given by: y = m + n,  (3.1)  where y is the known noisy data and n is zero-centered Gaussian noise (Hennenfent and Herrmann, 2006). Depending on the processing applied to the data beforehand, the noise can be white or colored. Given our main goal of recovering m without loss of coherent reflection information, we expand m in the curvelet domain; i.e., we have y = CHx + n,  (3.2)  where x is the curvelet transform vector and CH is the curvelet transform synthesis operator. Here, the symbol H denotes conjugate transpose that is equivalent to the inverse for this choice of curvelet transform. The forward and inverse curvelet transform operators are taken from Curvelab codes available on www.curvelet.org. First, we detect the curvelet coefficients whose absolute values exceed a threshold level selected such that the residue, i.e., the difference between the noisy data and the estimate for m, contains as much noise as possible without touching the coherent reflection energy. Because the curvelet transform is redundant, the use of a threshold level does not lead to a sparse curvelet vector. Therefore, we include a second iterative stage during which additional noise is removed while preserving the reflection amplitudes. Notice that this procedure differs from most common denoising approaches that only apply one of these two stages.  53  However, our requirement of not removing coherent energy calls for this combined approach. After detecting the set of curvelet coefficients that survive the threshold, we solve the following constrained optimization problem: subject to ||y − Ax||2 ≤ ε, where  (3.3)  represents the estimated curvelet transform coefficient vector, and ε is  proportional to the noise level (Elad et al., 2005). During the solution of this optimization problem, we restrict ourselves to the detected curvelet coefficients by restricting the curvelet synthesis matrix, i.e., A = CHR =  , which corresponds to compounding the  inverse curvelet transform CH with a diagonal matrix R that has ones at the entries of the detected coefficients and zeros otherwise. The noise level is also set in accordance with the noise estimate from the detection stage; i.e., ε =  .  Empirically, we obtained the best results, i.e., results that look best in the so-called ’eye-ball norm’, by choosing the curvelet transform to have 5 scales, 16 angles at the second coarsest level and 64 angles at the finest (fifth) scale. Figure 3.2 shows the F-K tiling for our choice of curvelet transform. In essence, the curvelet transform entails a tiling of the frequency plane into multiple-angle and dyadic multiscale wedges (Fig. 3.2). Following the procedure summarized above, the estimate for the denoised data, , is found amongst the detected coefficients by minimizing the L1 norm. Following the seminal work by Figueiredo and Nowak (2003), Daubechies et al. (2004) and Elad et al. (2005), we solve the constrained optimization problem (cf. Eq. 3.3) by a series of unconstrained optimization problems for decreasing regularization parameters λ: (3.4)  The parameter λ determines the trade-off between data consistency and the sparsity in the curvelet domain. During the cooling procedure, the trade-off parameter λ is lowered gradually. This allows more curvelet coefficients to enter into the solution until λ is  54  Figure 3.2. Discrete curvelet partitioning of the 2-D Fourier plane into second dyadic coronae (boxes) and sub-partitioning of the coronae into angular wedges (radiating lines). The tiling corresponds to 5 scales, 16 angles at the second coarsest level and 64 angles at the finest (fifth) scale. The angles double every other scale. Numbers in the upper left indicate that if the length of an angular wedge segment is proportional to 2j, where j is an integer, then the end wedge dimension is proportional to 2j/2. (Adapted from Hennenfent and Herrmann 2006).  Figure 3.3. The cooling method with iterative thresholding to solve equation (3). L is the number of inner iterations, K is the number of outer iterations, k is the counter for the outer iterations, and λk are the series of trade-off parameters with λ1 being the largest; other symbols are as explained in the text. Note that this is a two-stage iterative procedure. (Adapted from Herrmann and Hennenfent 2008).  55  sufficiently small such that the stopping criterion ||y − Ax||2 ≤ ε is met. Each Kth subproblem is solved by L iterations of soft thresholding (Figueiredo and Nowak , 2003; Daubechies et al., 2004) that involves the following updates for x: x ← where  (x + AH (y − Ax))  (3.5)  (x) = sgn(x) · max(0, |x| − |λk|), is the elementwise soft thresholding operator  (Donoho, 1995). The descent update, x ← (x + AH(y − Ax)), minimizes the quadratic part of Eq. 3.4 while the soft thresholding promotes sparsity. We solve a series of problems (Eq. 3.4) starting with high values of λ (λ1 is set slightly below the largest curvelet coefficient, i.e., λ1 < ||AHy||∞) and decreasing the value of λ until ||y − Ax||2 ≤ ε. Figure 3.3 summarizes the algorithm. In contrast to F-X deconvolution and the localized slant stack method, the algorithm does not work within the context of a sliding window that has a fixed width in the x-direction. Rather, the entire data are treated as a single image and decomposed into localized functions of different frequencies, angles and positions. All these curvelet functions have a weight associated with them to model the desired coherent data. The weights are estimated by the algorithm.  3.3 Application to synthetic data Prior to the application of the curvelet denoising algorithm to synthetic data, the set of decreasing trade-off parameters (λ1, λ2, …, λk ); L, the number of inner iterations; and K, the number of outer iterations, must be selected. To choose the best threshold level for the restriction in the curvelet domain (Cr), a test in which a graph of signal-to-noise ratio (SNR) versus threshold level is plotted (Fig. 3.4). The SNR is calculated from the following formula SNR =  where m is the true model (noise free data) and  (3.6)  is the estimated model (denoised data).  The test was conducted on a synthetic shot gather computed from a salt dome model as discussed later. The resulting graph for white noise shows a fairly rapid rise to a peak value of SNR followed by a more gradual decrease (Fig. 3.4). This is expected because 56  the initial threshold level discards the noisy component of the data, resulting in an increase of the SNR value. After a certain point (the arrow in Fig. 3.4), the threshold level starts harming the signal components causing a decrease in SNR value. Based on this test, the threshold value corresponding to the peak value of SNR is selected. In practice we compute the parameter ε from the noise estimated according to the procedure explained earlier. However, here we use the exact value of noise level (ε = ||n||2) since it is known in the synthetic example.  Figure 3.4. Signal-to-noise ratio (SNR) vs. threshold values for white noise for synthetic shot gather from salt dome model (see Fig. 3.6). The arrow points to the peak SNR value from which the corresponding threshold value λ is selected. In the examples in this paper, white or colored noise is used.  The values of K = 2 (two outer iterations) and L = 1 (one inner iteration) were selected after exhaustive testing in order to get the best possible SNR value. We observed that increasing the number of iterations caused a decrease in SNR values because the iterative thresholding started introducing bias in the amplitudes of the estimated curvelet coefficients. This bias corresponds to harming the components that represent coherent energy. Since our chosen criterion is to avoid harming the coherent energy to the maximum extent possible, the number of iterations is limited.  57  3.3.1 Application to synthetic shot gathers To exemplify the capability of the curvelet method to remove incoherent noise while leaving coherent reflectivity relatively untouched, the algorithm is applied to two synthetic shot gathers, one derived from a crustal velocity model and one from a basin salt dome model. The synthetic shot gather of Fig. 3.5a is computed using the finite difference algorithm on a model comprising the deterministic 2-D crustal velocity model from Lithoprobe’s SNORCLE experiment (Fernandez Viejo and Clowes, 2003) augmented with a 3-km-thick complex Moho transition zone described by Oueity and Clowes (2010). The latter was added to provide some detail to the simple crustal model. Most of the high amplitude phases from 0 – 4 s do not represent reflected energy. One shallow reflection at 1.2 s (due to a known geological boundary in the model) and three deeper reflected phases are generated. Those at ~3.4 s and ~8.0 s are due to small step increases in velocity in the crustal model. Discontinuous reflectivity from 11.0 to 11.5 s results from the insertion of the Moho transition zone. To generate a shot gather with random noise, zero-centered white Gaussian noise with a standard deviation σ = 50, yielding a noise level ε  11642 is added to the shot gather (Fig. 3.5b; Table 3.1). The  curvelet denoising algorithm (Fig. 3.3) is applied to the noisy data to generate a gather in which the noise has been removed (Fig. 3.5c). Comparison of Fig. 3.5c with Fig. 3.5a shows that the incoherent noise has been effectively removed and much of the coherent reflectivity has been retained. This is also illustrated in Fig. 3.5d which shows primarily noise with a very small component of the high amplitude phases at less than 4s. Fig. 3.5e shows the reconstruction error of the recovered model; i.e., Fig 3.5(a) – Fig 3.5(c). The reconstruction error shows that the algorithm is not perfect and numerical errors do creep in when applying the algorithm. However the signal-to-noise ratio (SNR) of 15.9 of the recovered model and the norm of the removed noise (ε  11627 ) show that the gather is  reconstructed to an accurate level (see also Table 3.1). The shot gather computed from the crustal model shows limited complexity because such models do not include detailed velocity information due to the nature of the interpretation procedure. To further illustrate the effectiveness of curvelet denoising, a  58  Figure 3.5. (a) Noise-free shot gather computed using the finite-difference algorithm on a crustal scale velocity model as described in text. (b) Shot gather with added Gaussian white noise. (c) Shot gather of (b) following application of the denoising algorithm. (d) Noise removed by the curvelet denoising algorithm; i.e., (b) minus (c). (e) Difference between the initial model and recovered model; i.e., (a) minus (c). The reconstruction error clearly shows some of the high amplitude phases before 4 s and weak remnants of the deeper reflectivity, implying that the algorithm is not 100% amplitude preserving.  59  more complex shot gather is computed from an industry-related velocity model (Fig. 3.6). The model consists of a high-velocity salt dome surrounded by sedimentary layers and a water bottom (Herrmann et al., 2007). Many reflections, some overlapping, are evident in the gather. Zero-centered white Gaussian noise with a standard deviation σ = 50, yielding a noise level ε  12800, is added to the shot gather (Fig. 3.6c; Table 3.1). The curvelet  denoising algorithm (Fig. 3.3) is applied to the noisy data to generate a gather in which the noise has been removed (Fig. 3.6d). Comparison of Fig. 3.6d with Fig. 3.6b shows that the incoherent noise has been effectively removed and the coherent reflectivity has been almost completely retained. Fig. 3.6e represents the noise that was removed during the curvelet denoising procedure (Table 3.1); no coherent reflectivity is observed. To further illustrate the effectiveness of the method, Fig. 3.6f shows the difference between the initial synthetic data (Fig. 3.6b) and those after denoising (Fig. 3.6d). Some weak coherent arrivals from the reflection curves and the splatter of errors within the data window that are associated with curvelets indicate that the algorithm does not perfectly preserve amplitude information. However, Fig. 3.6d demonstrates visually that the denoised gather is almost identical with the original synthetic data; see Table 3.1. Noisy data  Recovered  Standard  Norm of  Norm of  (SNR)  model  deviation of  noise added  noise  (SNR)  noise added  removed  (σ) Crustal shot gather  2.9  15.9  50  11642  11627  Salt dome shot  3.6  14.3  50  12820  12600  4.3  14.9  50  4675  4617  gather Marmousi model  Table 3.1. Noise parameters for synthetic models.  3.3.2 Application to synthetic post-stack data For the post-stack synthetic study, we chose to use data from an exploration-oriented model because crustal velocity models derived from refraction and wide-angle reflection data do not provide the detail that leads to adequate complexity to effectively illustrate the efficacy of curvelet denoising. This is illustrated also by Fig. 3.5. Thus, the synthetic 60  Figure 3.6. (a) Subset of the 2-D salt dome velocity model (after Herrmann et al. 2007). (b) Noise-free synthetic shot gather computed using finite-difference modeling on the velocity model. (c) Shot gather with added Gaussian white noise filtered between 5 and 60 Hz. (d) Shot gather of (c) following application of the denoising algorithm. (e) Noise removed by the curvelet denoising algorithm; i.e., (c) minus (d). (f) Difference between the initial model and recovered model; i.e., (b) minus (d). The reconstruction error shows some curvelets implying that the algorithm is not 100% amplitude preserving.  61  post-stack data were generated from a subset (5.0 to 6.6 km distance and 1.4 to 3.0 km depth) of the industry-standard Marmousi model (Versteeg, 1994). The model is a complex, synthetic 2D acoustic model that has been applied in a number of industry tests and comparisons. Figure 3.7a shows the normal incidence reflectivity derived for our application. For simplification and the purpose of this test experiment, we assumed that the depth axis was the time axis with a sampling rate of 4 ms; 400 traces were generated across the model. Near-surface complexities, a prominent pinch-out associated with a reservoir and curved reflections with contrasting dips are prominent. Colored Gaussian noise was added to the section (Fig. 3.7b; Table 3.1); some of the weaker reflectivity is masked by the noise. The curvelet denoising algorithm as used in the shot gather tests was applied to provide a denoised section (Fig. 3.7c). Comparing Figs. 3.7a and 3.7c, the random background noise has been removed and almost all of the coherent reflective features have been retained. This is shown more explicitly by considering the small data subset at the tip of the pinch-out (rectangles in Figs. 3.7a to 3.7c). Figures 3.7d, e, f and g show the synthetic reflectivity, the reflectivity contaminated with noise, the reflectivity recovered after denoising and the difference between Figs. 3.7d and 3.7f, respectively. Figure 3.7f demonstrates the extent to which the noise has been removed, leaving a section that is very similar to the original data. All the major reflectivity, including the pinch-out, is preserved. However, the weak reflection at 1.05 s in Fig. 3.7d is not recovered, although an equivalently weak dipping reflection at 1.3 s is visible on the denoised section of Fig. 3.7f. The reconstruction error (Fig. 3.7g) demonstrates that virtually none of the desired coherent energy is missing from the denoised data section; also see Table 3.1.  3.4 Application to crustal reflection data As part of multidisciplinary studies in Lithoprobe’s Slave-Northern Cordillera Lithospheric Evolution (SNORCLE) transect, near-vertical incidence crustal reflection data were recorded along a 725 km long profile over Paleoproterozoic and Archean domains in the Northwest Territories of Canada (Fig. 3.8). Cook et al. (1999) provide information on the survey and a complete interpretation. For the seismic data acquisition  62  Figure 3.7. (a) Synthetic post-stack data, noise free, from a subset of the Marmousi model (see text). The small rectangle shows the location of the data enlargement displayed in (d). (b) Data of (a) with added Gaussian white noise filtered between 5 and 60 Hz. The small rectangle shows the location of the data enlargement displayed in (e). (c) Data of (a) following denoising with the curvelet algorithm. The small rectangle shows the location of the data enlargement displayed in (f). (d) Data enlargement from rectangle in (a). (e) Data enlargement from rectangle in (b). (f) Data enlargement from rectangle in (c). Note the extent to which the curvelet denoised data in (f) has recovered the features of the noise-free data in (d). (g) Difference between the initial model and the recovered model; i.e., (d) minus (f).  63  parameters, refer to Table 3.2. Processing to a stacked seismic section followed standard procedures used in Lithoprobe studies (Cook et al. 1999). For the present example, data from a 60 km segment of the profile over the Paleoproterozoic domains were used (Fig. 3.8).  Figure 3.8. Tectonic map of northwestern Canada showing location of SNORCLE Line 1 (heavy black line outlined in white). The white rectangle highlights the part of the line from which data for the present study were taken. Sediments of the Phanerozoic Western Canada Sedimentary Basin overlie the Precambrian domains west of the long dashed line. The thin black line with teeth identifies the Cordilleran deformation front. Short dashed black lines show political boundaries. CS - Coronation Supergroup, SD - Sleepy Dragon Complex, GSLsz - Great Slave Lake shear zone, AB - Alberta, BC - British Columbia, NWT Northwest Territories, YK - Yukon. Inset shows location of map within Canada.  For the application of curvelet denoising to real data, the same number of scales and angles that were used for the synthetic data were retained (Fig. 3.2). Choosing the value of the noise level, ε, is a challenge for real data. In this study, the noise removed for a series of threshold values and associated noise levels was observed. From these experiments, the result that caused reasonable noise attenuation, and for which the difference between the observed and denoised data appeared as random noise, was 64  selected; this is known as the best threshold result. Based on this result, the restricted curvelet transform  and ε were defined. In other words, this procedure allows only  those curvelet coefficients that survive the best thresholding to enter the solution. As demonstrated by the examples that follow, choosing the parameters in this manner caused minimum harm to coherent reflected energy. Parameter  Value  Number of vibrators  4 or 5  Sweep frequencies  10 - 80 Hz  Sample rate  4 ms  Sweep length  20 s  Record length  32 s (uncorrelated)  Number of channels/record  404  Receiver spacing  60 m  Vibrator point spacing  90 m  Geophone frequency  10 Hz  Nominal fold  134  Table 3.2. Data acquisition parameters.  3.4.1 Application to a vibroseis shot gather In standard processing, crustal reflection studies rely heavily on the redundancy (multifold aspect) of the data to improve signal-to-noise ratios and generate good images for interpretation; nominal fold for the SNORCLE data set was 134. Normally, single shot gathers do not show strong reflectivity, particularly for the lower crust and Moho. However, the quality of the data recorded along the SNORCLE profile was exceptional and on a few shot gathers coherent reflectivity to the base of the crust was observed (Fig. 3.9a). Nevertheless, this raw shot gather is still contaminated with bad traces, ground roll and random and other noise that perturb or mask many features, especially deeper reflections. As shown, the shot gather also is characterized by a gradual decrease in reflectivity with time, which is due to weak reflections from the deeper levels. The coherent reflection across the gather at about 11 s is related to the crust-mantle transition or Moho. 65  Figure 3.9. Shot gather from SNORCLE Line 1 segment highlighted in Fig. 7. (a) Noisy data without AGC. GR – ground roll, M – Moho reflection. (b) Noisy data with AGC. (c) Result of applying curvelet denoising (inversion approach) on the shot gather with AGC. Moho reflection, M, is clearer. (d) Noise removed [i.e., (b) minus (c)] by curvelet denoising. Note that ground roll, GR, is removed as noise.  To deal with high amplitude ground roll and bad traces (left in to demonstrate the efficacy of curvelet denoising), the amplitudes of the shot gather were balanced with an automatic gain control (AGC) using a 500 ms window; see Fig. 3.9b. Following the same procedures as used on the synthetic data, the curvelet denoising algorithm was applied to the shot gather with AGC; results are shown in Fig. 3.9c. Reflectivity within the shot gather is greatly improved. Many more coherent reflections, including dipping ones, are observed and none of the desired coherent energy was removed by the process. Some of the “bad trace” information was left, but continuity across the gather is greatly improved. Figure 3.9d is a plot of the noise removed by curvelet denoising. Note that not only the 66  random noise but also the ground roll of Fig. 3.9a and much of the vertical “bad trace” energy have been removed. 3.4.2 Application to a stacked seismic section To demonstrate the application of curvelet denoising on a stacked seismic section, a 60 km long segment of the SNORCLE profile was extracted from the original stacked section (Fig. 3.10a; location in Fig. 3.8). Standard Lithoprobe processing applied to the data to generate the section is described in Cook et al. (1999). The stacked unmigrated seismic section (Fig.3.10a) shows reflectivity throughout the crust. The crust-mantle transition is identified at 10.8-11.8 s (34-38 km, arrows 2 and 3). Reflectivity in the mantle below this level shows a decrease in amplitude. As observed in Fig. 3.10a, the clarity of the recorded seismic data is degraded by the low signal-to-noise ratio, which causes a fuzzy or blurred image and precludes detailed interpretation of the reflectivity. Due to high background noise, the crust-mantle transition becomes more diffuse on the southwestern end of the segment (arrow 3 in Fig. 3.10a). Some northeast-dipping reflectivity can also be identified on the southwestern end of the section at 4.0 and 6.5 s but can only be followed over very short distances (arrow 1 in Fig.3.10a). Following the same procedures as before, curvelet denoising was applied to the stacked data of Fig. 3.10a; see Fig. 3.10b. To demonstrate the performance of the algorithm at removing random noise and not the desired signal, a difference image derived from the stacked section and the curvelet denoised data was generated; see Fig. 3.10d. As observed in Figs. 3.10a and 3.10b, the improvement in image quality is clear when comparing the two data sets before and after curvelet denoising. The denoised seismic image (Fig. 3.10b) shows a highly reflective crust with dipping reflective features (arrow 1 and between block arrows), a relatively flat Moho (arrow 2) with an offset of about 0.8 s at about CDP 1100 to a deeper position (arrow 3) and a transparent upper mantle. The lower crust contains a series of discrete, northeast-dipping reflections within a zone that is about 3.0 s thick (about 10 km between block arrows) and extends to the Moho. Below about 9.5 s at the southwestern end of the section, the reflectivity becomes sub-horizontal  67  68  Figure 3.10. (a) A 60-km-long segment of the original stacked seismic section from SNORCLE Line 1 (Fig. 7) that is contaminated with incoherent noise. Arrow 1 identifies northeast-dipping reflectivity; arrow 2 identifies clear Moho reflections; arrow 3 shows diffuse Moho reflection. (b) Result of applying curvelet denoising algorithm on the stacked section of (a). Block arrows bound a region of dipping reflectivity; arrows 2 and 3 identify same features as in (a) but with more clarity. (c) Result of applying F-X prediction algorithm on the stacked section of (a); arrows 1, 2 and 3 identify same features as in (a) but with more clarity. However, dipping reflectivity as seen in (b) is not clearly observed. (d) Noise removed by curvelet denoising. Arrows refer to same features as in (a) and (b). (e) Noise removed by F-X prediction. Arrows 1, 2 and 3 refer to same features as in (a) and (c); arrows 4 identify subtle coherent reflectivity that has been removed by the F-X processing.  69  to about the position of Moho offset. These events were not visible on the original stacked image. The Moho is characterized by a sharp, narrow band of reflections (200 300 ms) that are piece-wise continuous, except for the offset, over the 60-km-long segment. The difference plot between the original data and the curvelet result (Fig. 3.10d) shows that the noise removed by curvelet processing is mostly incoherent in nature. Thus, the curvelet denoising technique achieves excellent noise attenuation without removing any coherent events.  3.5 Discussion The shot gather example provided in section 3.4.1 (Fig. 3.9) demonstrates the efficacy of the curvelet denoising algorithm for removing unwanted coherent energy as well as random noise. Such results may enhance the ability of crustal seismologists to examine reflectivity on individual shot gathers, for example the Moho arrivals, instead of just on stacked data, which is the current procedure. To this end, curvelet denoising has been applied to many shot gathers from SNORCLE Line 1, these being used in a quantitative study of Moho reflectivity to further understand the characteristics of the crust-mantle transition (Oueity and Clowes, 2010). The study has demonstrated the nature of lateral and vertical heterogeneities within the crust-mantle transition and enabled inferences about the tectonic processes involved in development of the Moho. Application of curvelet denoising to the 60-km segment of SNORCLE stacked, unmigrated data (Fig. 3.10a) generated a significantly improved image (Fig. 3.10b) that would be better for interpretation than the former. A revised, and improved, interpretation of Line 1 would require application of the denoising procedure to the entire 725-km-long seismic line, an undertaking that is beyond the scope of this contribution. As noted in the Introduction, predictive deconvolution is a common technique to attenuate noise. To generate an image for comparison with the curvelet denoised image (Fig. 3.10b), an F-X prediction algorithm was applied to the stacked data. The F-X predicted image (Fig. 3.10c) shows most of the features observed on the curvelet denoised image, but the northeast-dipping reflectivity and that from the Moho transition are not as clearly expressed. Within the prominent Moho reflectivity, individual and 70  coherent reflectivity across tens of traces, but discontinuous across the section, are delineated better on the denoised section of Fig. 3.10b compared to those on the F-X section of Fig. 3.10c. In addition, the noise removed by F-X prediction (Fig. 3.10e) contains some coherent energy (arrow 4). Therefore, we claim that results from curvelet processing are in many aspects superior to those from F-X prediction. Currently, interpretation of seismic sections (stacked, or stacked and migrated) in crustal reflection studies make use of coherency filtering to provide an improved image for interpretation and display in publications (Kong et al., 1985; Milkereit and Spencer, 1989). It was used on the SNORCLE Line 1 section so that the geometrical aspects of the reflectivity were well displayed (Cook et al., 1999). The denoised section of Fig. 3.10b does not quite show the clarity of reflections observed on the equivalent coherencyfiltered section (Fig. 3.2 in Cook et al., 1999), although all the major aspects of the reflectivity can be seen. However, the curvelet-denoised section preserves the amplitude information of the original data section, whereas all amplitude information is lost through coherency filtering. Thus, curvelet denoising would be particularly useful for enhancing data where further analyses that involve amplitude information are intended. Such analyses might include instantaneous amplitude, trace envelope amplitude, instantaneous frequency, etc. (e.g., Robertson and Nogami, 1984; Robertson and Fisher, 1988; Barnes, 1996). Since the amplitude information on reflection seismic sections are not usually utilized in any quantitative sense, the ability to do so with curvelet-denoised data opens up new opportunities for studying and interpreting crustal reflectivity.  3.6 Conclusions By exploiting the ability of curvelets to detect wavefront-like features and sparsely represent seismic data, a denoising procedure that is a clear improvement over results from existing methods such as adaptive F-X predictive-error filtering has been developed. The latter methods often dim reflection amplitudes and are subject to leaving incoherent noise present. The two-stage procedure applied in curvelet denoising largely eliminates these shortcomings by first detecting a set of curvelet coefficients that are mostly noisefree and contain all the reflection information. During the second iterative stage, as much 71  of the remaining noise as possible is removed while preserving the amplitudes. Application of curvelet denoising to synthetic shot gathers and a stacked section exemplify the efficacy of the method. The technique also was applied to observed shot gathers and a stacked seismic section from a Lithoprobe survey in northwestern Canada. Clear coherent reflectivity to the base of the crust on the shot gather and dipping reflection events and reflections from the Moho that are much more clearly expressed in the post-stack example demonstrate that the method also works well on real data. Curvelet denoising expands the toolbox available for seismic data processing. Although the algorithm was developed for 2D data, it could be extended, using 3D curvelets, to enhance the quality of 3D datasets so that structure in three dimensions can be better exploited. This would require further developments. Another approach for further development would consider choosing a restriction operator that suppresses certain dipping events within the input data. This could lead to a sophisticated dip filtering algorithm using the curvelet transform.  72  Bibliography Abma, R., Claerbout, J., 1995. Lateral prediction for noise attenuation by t-x and f-x techniques. Geophysics 60, p. 1887–1896. Barnes, A.E., 1996. Theory of 2-D complex seismic trace analysis. Geophysics 61, p. 264-272. Candés, E., Demanet, L., Donoho, D., Ying, L., 2006. Fast discrete curvelet transforms. Multiscale Modeling and Simulation, 5, p. 861–899. Candés, E., Demanet, L., Donoho, D., Ying, L., 2004. New tight frames of curvelets and optimal representations of objects with C2 singularities. Communications on Pure and Applied Mathematics 57, p. 219–266. Cook, F.A., van der Velden, A.J., Hall, K.W., Roberts, B.J., 1999. Frozen subduction in Canada’s Northwest Territories: Lithoprobe deep lithospheric reflection profiling of the western Canadian Shield. Tectonics 18, p. 1–24. Daubechies, I., Defrise, M., De Mol, C., 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics 5711, p. 1413–1457. Donoho, D., 1995. De-noising by soft thresholding. IEEE Transactions on Information Theory 41, p. 613–627. Elad, M., Starck, J. L., Querre, P., Donoho, D. L., 2005. Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA). Applied and Computational Harmonic Analysis 19 (3), p. 340–358. Fernandez Viejo, G., Clowes, R.M., 2003. Lithospheric structure beneath the Archaean Slave Province and Proterozoic Wopmay orogen, northwestern Canada, from a LITHOPROBE refraction/wide-angle reflection survey. Geophysical Journal International 153, p. 1-19. Figueiredo, M. Nowak, R., 2002. An EM algorithm for wavelet-based image restoration. IEEE Transactions on Image Processing 12, p. 906-916. Fomel, S., 2002. Applications of plane-wave destruction filters. Geophysics 67, p. 1946–1960. Herrmann, F. J. Hennenfent, G., 2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International 173, p. 233–248.  73  Hennenfent, G., Herrmannn, F. J., 2006. Seismic denoising with nonuniformly sampled curvelets. Computing in Science and Engineering 8 (3), p. 16–25. Herrmann, F. J., Boeniger, U., Verschuur, D. J., 2007. Nonlinear primary-multiple separation with directional curvelet frames. Geophysical Journal International 170, p. 781–799. Jones, I.F., Levy, S., 1987. Signal-to-noise enhancement in multichannel seismic data via the Karhunen-Lo´eve Transform. Geophysical Prospecting 35 (1), p. 12–32. Kong, S. M., Phinney, R. A., Roy-Chowdhury, K., 1985. A nonlinear signal detector for enhancement of noisy seismic record sections. Geophysics 50, p. 539-550. McMechan,G. A., 1983. P-x imaging by localized slant stacks of T-x data. Geophysical Journal of the Royal Astronomical Society 72, p. 213–221. Milkereit, B., 1987. Migration of noisy crustal seismic data. Journal of Geophysical Research 92, p. 7916-7930. Milkereit, B., Spencer, C., 1989. Noise suppression and coherency enhancement of seismic data. In Agterberg, F.P. & Bonham-Carter, G.F.,eds., Statistical Application in the Earth Sciences, Geological Survey of Canada, p. 243–248. Neelamani, R., Baumstein, A. I., Gillard, D. G., Hadidi, M. T., Soroka, W. L., 2008. Coherent and random noise attenuation using the curvelet transform. The Leading Edge 27, p. 240–248. Olhovich,V. A., 1964. The causes of noise in seismic reflection and refraction work. Geophysics 29, p. 1015–1030. Oueity, J., Clowes, R.M., 2010. The nature of the Moho transition in NW Canada from combined near-vertical and wide-angle seismic reflection studies. Lithosphere, 2, p. 377-396. Robertson, J.D., Fisher, D.A. 1988. Complex seismic trace attributes. The Leading Edge 7, p. 2226. Robertson, J.D., Nogami, H.H. 1984. Complex seismic trace analysis of thin beds. Geophysics 49, p. 344-352. Tarantola, A., 1987. Inverse Problem Theory. Elsevier, N.Y. Ulrych, T.J., Sacchi M.D., Graul, J.M., 1999. Signal and noise separation, Art and science. Geophysics 64, p. 1648- 1656.  74  Versteeg, R., 1994. The Marmousi experience: Velocity model determination on a synthetic data set. The Leading Edge 13, p. 927-936 Zhang, R., Ulrych, T. J., 2003. Physical wavelet frame denoising. Geophysics, 68, p. 225–231.  75  Chapter 4  The nature of the Moho transition in NW Canada from combined near-vertical and wide-angle seismic reflection studies3 4.1 Introduction The crust-mantle boundary, or the Moho, is one of the most distinct manifestations of a differentiated Earth. Major changes in petrology, mineralogy, chemistry, seismic wave velocity, density, and rheology occur across it (Jarchow and Thompson, 1989). Originally, the Moho was defined on the basis of seismic head waves (Mohorovičić, 1910). It was described as the depth at which the compressional P-wave velocity increases rapidly or discontinuously to 7.6-8.6 km/s (Steinhart, 1967). This suggested a petrologic interpretation of the Moho as a first-order (i.e., zero thickness) discontinuity in rock composition from predominantly mafic lower crustal rocks to predominantly ultramafic upper mantle rocks. Ophiolite studies indicate that the lithologic sequence of the oceanic Moho includes pelagic sediments, pillow basalts, massive and layered gabbro and residual ultramafic tectonites of the upper mantle (Casey et al., 1981; Karson et al., 1984; Boudier and Nicloas, 1995; Jousselin and Nicolas, 2000; Dilek et al., 2008). These rock cumulates form laterally discontinuous lenses with their thicknesses ranging from 1 cm to several tens of meters and their aspect ratios ranging from 10:1 to 100:1. Prominent reflections from within the oceanic Moho transition zone obtained by seismic reflection methods also support a geologically complex boundary (e.g., Hasselgren and Clowes, 1995; Nedimović et al., 2005). 3  A version of this chapter has been accepted for publication: Oueity, J., and Clowes, R. M., 2010. The  nature of the Moho transition in NW Canada from combined near-vertical and wide-angle seismic reflection studies, Lithosphere, 2, p. 377-396.  76  Unlike abundant ophiolites, which allow direct observation of the oceanic Moho, exposed sections of the continental Moho are rare (Fountain and Salisbury, 1981). Nevertheless, important constraints on the nature of the continental Moho are provided by the small continental dataset of xenoliths (McGetchin and Silver, 1972; DeBari et al., 1987; Kopylova et al. 1998). Xenoliths samples indicate a lower crust composed of mafic granulite, gabbro and amphibolite and an upper mantle composed of several varieties of peridotite with some eclogite. Despite all the geologic studies of the Moho by means of exposed sections (ophiolites) or xenolith samples, seismic data still remain the main source for our understanding of the structure of the Moho and its development. As techniques applied to lithospheric studies have improved considerably in terms of seismic data acquisition technologies (e.g., Meissner, 1986; Mooney, 1987; Mjelde et al., 1997; Pancea et al., 2005) and computational capabilities (e.g., Fuchs and Muller, 1971; Spence et al., 1984; Levander and Holliger, 1992; Zelt and Smith, 1992; Carbonell et al., 2002), it has become apparent that the Moho is typically not a uniform, sharp, first-order discontinuity, but rather a complex and variable transition zone. From the seismic perspective, a frequently cited model for the Moho has been characterized by a variable thickness transition zone comprised of anastomosing layers progressing with depth from mafic to ultramafic rocks (Jarchow and Thompson, 1989; Braile and Chiang, 1986). Alternatively, the Moho transition may represent a change in scale of vertical layering and horizontal extent of inhomogeneities (Tittgemeyer et al., 1999; Hurich, 2003; Carpentier and RoyChowdhury, 2007). Although the Moho is well imaged by near-vertical incidence (NVI) and refraction/wide-angle reflection (R/WAR) data, the exact nature and characteristics of the transition are not well understood. The wavefield (i.e., waveform, duration and amplitude) of a seismic signal (reflected or refracted) from the lower crust and upper mantle may contain significant information about the structure of the Moho transition zone, its fabric and scale of heterogeneities. Only a few studies have specifically addressed this issue. These studies can be grouped into two categories, qualitative and quantitative. Examples of qualitative analyses of reflection patterns near the Moho are 77  those of Gibbs (1986) on the Consortium for Continental Reflection Profiling (COCORP) seismic  reflection  profiles,  Meissner  (2000)  using  Deutsches  Kontinentales  Reflexionsseismisches Programm (DEKORP) data, and Hammer and Clowes (1997), Cook (2002) and Eaton (2006) on LITHOPROBE reflection transects that sample regions of diverse age and tectonic history. From these studies, reflection patterns appear to be divided into three primary classes: 1) a distinct, multicyclic band of reflectors that are laterally continuous or piece-wise continuous over tens of kilometers, 2) fading out of crustal reflectivity with depth, and 3) no clear reflections in the vicinity of the Moho. These studies did not involve wide-angle data. Although the comparisons did not clearly link a specific tectonic process with a distinct reflection pattern, sharp Moho reflections were associated with lower crustal deformational processes including compression, extension, and intrusive features. Quantitative studies use synthetic modeling to construct the seismic response of the continental Moho using seismic velocity data (laboratory, seismic refraction results) together with petrologic information derived from exposures of the crust-mantle boundary when available. In earlier studies, the Moho was modeled as a finite thickness laminated zone of alternating high and low velocity layers (Clowes and Kanasewich, 1970; Hale and Thompson, 1982; Fountain, 1986; Braile and Chiang, 1986). More recently, the Moho has been considered as a laterally and vertically heterogeneous transition where the heterogeneities are represented by randomly distributed ellipsoids (Carbonell et al., 2002; Tittgemeyer et al., 2003). Another representation of the Moho follows a statistical approach, in which the transition zone is characterized by a stochastic heterogeneity distribution (Hurich, 2003; Nielson and Thybo, 2006; Carpentier and RoyChowdhury, 2007). Only a few of these studies included consideration of wide-angle data (Carbonell et al., 2002; Tittgemeyer et al. 2003; Nielson and Thybo, 2006). In this paper we analyze the near-vertical incidence and wide-angle seismic data acquired along Line 1 of Lithoprobe’s Slave-NORthern Cordillera Lithospheric Evolution (SNORCLE) transect in the Paleoproterozoic-Archean domains of the Northwest Territories, Canada (Cook et al. 1999; Fernandez Viejo and Clowes, 2003). We investigate the dynamic characteristics of Moho reflections by calculating the 78  synthetic seismic signature for both near-vertical and wide-angle data of a number of crust-mantle transition models using 1- and 2-D wave propagation algorithms. The models range from a simple first-order discontinuity to more complex, laterally and vertically heterogeneous structures. Unraveling the detailed structure of the crust-mantle boundary and the way in which it forms can help us understand some major tectonic processes such as crustal growth, accretion and delamination (Morozov et al., 2001; Carbonell et al., 2002).  4.2 Tectonic background The Wopmay orogen in northwest Canada (Fig. 4.1) is a Paleoproterozoic assembly of domains that developed during the Calderian orogeny between 1.92 Ga and 1.84 Ga (Hoffman, 1989). The orogen is divided into four major north-south trending domains: the Coronation Supergroup, the Great Bear magmatic arc, the Hottah terrane and the Fort Simpson terrane. Rifting to the west of the Slave province at about 1.90 Ga initiated deposition of the Coronation Supergroup, a west-facing shelf-rise, which shortly thereafter was intruded by a suite of plutons and translated eastward by the Hottah collision to form a foreland thrust-fold belt (Hoffman, 1988). The Hottah terrane formed as a magmatic arc and collided with the western Slave craton ca. 1.89 – 1.88 Ga (Hoffman and Bowring, 1984; Hildebrand and Bowring, 1999). The collision caused compression, shortening and eastward thrusting of the Coronation Supergroup onto the Slave craton (Hoffman and Bowring, 1984). The lack of coeval arc magmas on the western Slave craton indicates that the accretion of the Hottah terrane onto the craton was the result of west-dipping subduction of an ocean basin. The Hottah terrane consists of Paleoproterozoic sedimentary and intermediate volcanic rocks metamorphosed to amphibolite grade and intruded by calc-alkaline granitic plutons during the period 1.940-1.902 Ga (Hildebrand et al., 1987). To the west, another oceanic plate was converging against the western Hottah terrane, generating an east-dipping subduction zone below the Hottah terrane and Coronation Supergroup. It led to the formation of the Great Bear magmatic arc during the period between 1.875 to 1.843 Ga (Hildebrand et al., 1987; Gandhi et al., 2001). 79  Figure 4.1. Tectonic map of the study area and location of SNORCLE near-vertical incidence reflection and refraction/wide-angle reflection Line 1. Stars identify shot locations; some are identified by numbers. Sediments of the Phanerozoic Western Canada Sedimentary Basin overlie the Precambrian domains west of the long dashed line. Short dashed black lines show political boundaries. CS - Coronation Supergroup, SD Sleepy Dragon Complex, GSLsz - Great Slave Lake shear zone, AB - Alberta, BC - British Columbia, NWT - Northwest Territories, YK - Yukon. Inset shows location of map within Canada.  Between the Coronation Supergroup and the Hottah terrane lies the 100 km-wide Great Bear magmatic arc, which consists of volcanosedimentary sequences and plutonic rocks (Hildebrand et al., 1987). Aeromagnetic data indicate that rocks of the Great Bear arc extend from north to south beneath the younger Proterozoic and Paleozoic cover for a total length of about 900 km (Coles et al., 1976). At its southern end, the magmatic arc is about 2.0–4.5 km thick as interpreted from seismic reflection data along SNORCLE Line 1 (Cook et al., 1999). These data show a strong reflection from the top (~0.2 s) of the arc, which has been drilled, and another strong reflection at ~1.5 s from its interpreted base. Whether the thin basin-like feature represents the complete Great Bear section at this latitude or whether there are additional Great Bear rocks below is not known. In either case, development of the arc does not appear to disrupt the older structure beneath it, at least within the resolution of the reflection survey. One possible explanation for the intact 80  structure underneath is that the arc magmas migrated through the crust along relatively narrow preexisting zones of weakness without disrupting the primary crustal structure (Cook et al., 1999). The entire Great Bear - Hottah - Slave assembly was deformed by a set of brittle conjugate faults, which have been interpreted as the result of the terminal collision of the Fort Simpson exotic terrane with the western margin of the Hottah terrane (Hoffman and Bowring, 1984). The collision of the Fort Simpson terrane with the Hottah terrane occurred after the formation of Fort Simpson magmatic rocks, ca. 1.845 Ga (Villeneuve et al., 1991) and pre-1.71 Ga mafic intrusions of sedimentary rocks (Abbott et al., 1997). Formation of the Fort Simpson basin, identified as such on SNORCLE Line 1 by seismic reflection and refraction data (Cook et al., 1999; Welford et al., 2001), is the result of lithospheric extension, which followed the collision of the westernmost terranes with the older Paleoproterozoic domains. From the reflection data, the basin was interpreted as a monocline with a base that dips westwards from a depth of ~3 km at the eastern end to a depth of ~24 km at the western end, revealing at least 20 km of subsurface relief (Cook et al., 1999). The enlarged continental assembly was affected by a number of crustal extension phases that continued to Neoproterozoic time (Thompson et al., 1987). The Proterozoic rocks of the Wopmay orogen along SNORCLE Line 1 are overlain by Phanerozoic sedimentary rocks that thicken from their feather edge near the eastern boundary of the Great Bear magmatic arc to about 1000 m over the Nahanni domain (Fig. 4.1; Cook et al., 1999). The Nahanni domain, west of the Fort Simpson terrane, is the westernmost component of the Wopmay orogen in the study area. It has no surface exposure and is only distinguished by its magnetic signature.  4.3 Seismic experiment and data 4.3.1 SNORCLE Line 1 The main purpose of the SNORCLE transect studies was to investigate the various processes involved in the westward growth of North America from the Archean to the 81  present using a multidisciplinary approach (Cook and Erdmer, 2005). Vibroseis nearvertical incidence reflection data were acquired along SNORCLE seismic line 1, one of four such lines in the transect. Line 1 is 725 km long and crosses the southwest Slave craton and Wopmay orogen (Fig. 4.1). A total of 404 geophones per record with a station spacing of 60 m were used. Shot points were every 90 m with sweep frequencies ranging between 10 and 80 Hz. Complete field acquisition parameters are given in Cook et al. (1999). For such data, the nominal vertical resolution is given by λ/4, where λ is the wavelength, and the horizontal resolution is given by the radius r of the first Fresnel zone, r = λz / 2 , where z is the depth (Yilmaz, 2001). Thus, at a Moho depth of ~32 km where the seismic velocity is ~ 6.8 km/s (Fernandez Viejo and Clowes, 2003) and the dominant frequency determined from Moho reflections is ~21 Hz, the vertical resolution is ~80 m and the horizontal resolution is ~2.2 km. Refraction/wide-angle reflection data were recorded from 13 explosive shots along the same SNORCLE Line 1. The total number of seismographs was 594 with station spacing of 1.2 km. Further details of the field acquisition can be found in Fernandez Viejo and Clowes (2003). Interpretation of the R/WAR data followed standard procedures (Zelt and Smith, 1992) and provided information on the 2-D seismic velocity structure. The estimated uncertainty in velocities based on procedures discussed by Fernandez Viejo and Clowes (2003) are given in their Table 2.3. For the WAR data, traveltime modeling indicated vertical resolution at lower crustal depths on the order of 1 to 2 km and horizontal resolution between 60 and 80 km (Fernandez Viejo and Clowes, 2003). 4.3.2 Processing of the reflection data The reflection data were re-processed following standard procedures with parameters (Table 4.1) that are similar to, but not identical with, those used by Cook et al. (1999). The main difference in our procedure is with handling noise. Incoherent noise present in the data corrupts the quality of the signal and may lead to misinterpretation. Thus, separation of signal and noise is an important step in data processing, particularly in the lower crust where the signal-to-noise (SNR) ratio is low. Traditional seismic processing takes advantage of the redundant (multifold) data to improve SNR. In this study, we 82  adopt an additional and new technique, based on the curvelet transform (Hennenfent and Herrmann, 2006; Neelamani et al., 2008), to suppress incoherent noise and improve the resolution of the seismic data at deeper levels, including the Moho (e.g., Fig. 4.2). Further discussion about the technique and its applications can be found in Chapter (3). The stacked and denoised section of the seismic image beneath the Great Bear magmatic arc shows a distinct Moho transition at ~11.0 s from about 15 to 60 km and a slightly deeper Moho at 11.6 s from 0 to 15 km distance (Fig. 4.2). The transition comprises a multicyclic band of reflections separating a highly reflective crust from a relatively transparent upper mantle. Ideally, we would like to simulate this Moho transition response by calculation of synthetic seismograms from model representations of the Moho. To simulate the stacked section accurately, calculation of a large number of shot gathers along a high- resolution 2-D velocity model is required. This process is computationally very expensive and a high-resolution 2-D velocity model is not available. Parameter Crooked line geometry Static correction Trace balance Bandpass filter Deconvolution Automatic gain control NMO correction Stack Curvelet denoising  Value 30 m x 2000 m bin Average amplitude scaled to be 1.0 10-40 Hz Filter length=150 ms; gap length=2.0 ms 500 ms window V(z) derived from refraction study Nominally 134 fold 5 scales and 64 angles  Table 4.1. Processing parameters for NVI data.  An alternative approach is to identify Moho reflections on individual shot gathers, then attempt to simulate these through synthetic seismogram modeling. Previously, this has not been attempted because individual shot gathers are characterized by low SNR ratios and deeper reflections, particularly at the base of the crust, are usually obscured by random noise. We overcame this problem by utilizing the aforementioned curvelet denoising technique to attenuate incoherent noise on the shot gathers. Our results show remarkable reflections at ~11.0 s that we interpret as the Moho. The reflections are  83  characterized by a narrow band of 5–9 cycles that are laterally variable and piecewise continuous for several kilometers or more (Fig. 4.3).  Figure 4.2. The seismic reflection profile beneath the Great Bear magmatic arc along SNORCLE Line 1. Strong Moho reflections are visible at ~11.0 s between 15 and 60 km and at ~11.5 s between 0 and 15 km. Block arrows bound a region of dipping reflectivity in the middle and lower crust. Curvelet denoising (see text) has been applied to the stacked seismic section; no migration. Inset shows location of the profile (white) with respect to Line 1.  4.3.3 Processing of the refraction data Processing of the refraction data involved editing/removing noisy traces and band-pass filtering (3–15 Hz). To offset distances of ~160 km, the crustal refracted phase, Pg, is the first arrival (Fig. 4.4). Prominent secondary wide-angle reflections from the crust-mantle boundary (PmP) beneath the Great Bear magmatic arc are observed on three shots (1104, 1105 and 1108). They are characterized by a prominent reverberatory pattern with 84  duration of ~0.5 s between 80 and 220 km offset (Fig. 4.4). At offsets greater than 160 km, the refraction phase Pn traveling within the uppermost mantle becomes the first, but rather weak, arrival. Because of the ambiguity of the Pn phase, further amplitude/frequency analysis of this phase was not pursued. A normal-moveout-corrected and stacked section of shot gathers on which the secondary phase PmP was observed showed clear, horizontal reflections at ~11 s that extend laterally over ~200 km below the Hottah terrane and Great Bear magmatic arc (Fernandez Viejo et al., 1999). These horizontal reflections correspond remarkably well with the horizontal Moho reflections observed on the NVI data.  Figure 4.3. Different shot gathers after curvelet denoising (see text). Note the distinct band of reflections at ~11.0 s. Enlargements of Moho reflections are shown in the insets.  85  Figure 4.4. Record section for shots 1104, 1105 and 1108. Data are band-pass filtered from 3-15 Hz and displayed as trace-normalized amplitudes. Pg is a crustal, first arrival, refracted phase. PmP is a prominent, second arrival reflection from the Moho observed a at offsets > 80 km, followed by ~0.5-s- long coda. Pn is a weak, refraction phase in the upper mantle observed at offsets > 160 km.  86  Figure 4.5. (a) Amplitude vs. offset of the observed PmP phase on shots 1104, 1105 and 1108. Amplitude values are normalized to unity. Amplitude behavior is erratic and no common peak is observed among the shots. (b) Normalized amplitude spectra of the observed PmP phase. The spectra are calculated on single traces at 100, 120 and 140 km offsets and over a 0.5 s window.  87  For amplitude analysis, we calculated the PmP amplitudes within a 0.5 s window length after each traveltime pick for this phase (Fig. 4.5a). Only amplitudes associated with pick uncertainties < 0.15 s were considered (Table 4.2). All trace amplitudes were initially scaled such that ‘true relative amplitude’ between data of various instruments types was achieved; see Ellis et al. (2002) for details. The amplitude-versus-offset plots for the three shots show high variability with no obvious distance at which a common amplitude peak is observed. Shot Point  Phase  Average Pick Uncertainty (s)  1104  PmP  0.125  Pn  0.15  PmP  0.125  Pn  0.14  PmP  0.14  Pn  0.15  1105 1108  Table 4.2. Average uncertainties associated with the picks of PmP and Pn phases in shots 1104, 1105 and 11108.  4.4 Synthetic seismic modeling 4.4.1 Methodology Both near-vertical and wide-angle reflection seismic data show very distinct Moho reflections from beneath the Great Bear magmatic arc. The dynamic characteristics of these reflections contain valuable information that can help us shed some light on the structure of the crust-mantle transition (Moho). In order to investigate wavefield characteristics (traveltime, duration, amplitude, frequency, etc.), we calculated the synthetic seismic signature of different crust-mantle transition models using 1- and 2-D wave propagation algorithms (Appendix A). We examined a suite of models ranging from simple to more complex transitions suggested by the many studies mentioned earlier. These include: first-order discontinuities, gradient transitions, layered transitions and laterally heterogeneous transitions. The key parameters of these models, which control the seismic signature of the synthetic seismograms, are the total thickness of the 88  transition zone, the velocity contrast between the transition zone and the surrounding material and the correlation length of the lateral heterogeneities. By modifying these parameters synthetic seismograms that best simulate the real data can be obtained. To examine the seismic signature of the first-order discontinuity, gradient and layered transitions we used a 1-D reflectivity modeling algorithm (Fuchs and Müller, 1971). For the laterally heterogeneous models, a 2-D viscoelastic finite difference algorithm (Robertsson et al., 1994) is required. The 1- and 2-D velocity models that were incorporated into the 1-D reflectivity and 2-D viscoelastic modeling, respectively, are based on the interpretation of the R/WAR study along SNORCLE Line 1 (FernandezViejo and Clowes, 2003). The crustal P-wave velocities increase from 5.9 km/s at the surface to 6.8 km/s at the crust-mantle boundary (~32 km depth). The velocity of the uppermost mantle is 8.1 km/s and linearly increases to 8.4 km/s at 45 km depth. Corresponding S-wave velocities are calculated using the ratios Vp/Vs = 1.78 and Vp/Vs = 1.82 for the crust and upper mantle respectively (Fernandez-Viejo et al., 2005), whereas density values are based on laboratory measurements by Christensen and Mooney (1995) and Christensen (1996). The P- and S-wave quality factors, Qp and Qs respectively, are more difficult to assign since no attenuation studies have been carried out in the region. Values of Qp = 500 and Qs = 200 for the crust and Qp = 1000, Qs = 400 for the upper mantle were used for generating the synthetic seismograms displayed. These are “generic” values typically used within the algorithms by other scientists who have used it for crustal studies. The Qp value for the crust is consistent with values of 400 ± 200 determined for the continental crust in England (Scheirer and Hobbs 1990). In general terms, Qs is about one half the value of Qp (Fowler, 2005). Although Qs must be specified within the algorithm, this study does not include shear-wave arrivals so the value specified is not significant. Values of the quality factor for the uppermost mantle are generally considered higher. The generic value of Qp = 1000 is consistent with that shown by Fowler (2005) for a general earth model. It should be noted that quality factors of 1000 or more are indicative of very limited attenuation. Moreover, changes in the values of Q by factors of one half  89  or two do not materially affect the dynamic characteristics of the seismograms, as established by additional calculations using different Q values. For the wide-angle data, a zero-phase Ricker wavelet is used as the source signal with a bandwidth of 4-16 Hz to simulate an explosive source. For the vibroseis data, which is characterized by higher frequencies, a zero-phase Ricker wavelet with a 20 Hz central frequency of a two octave wavelet (10-40 Hz) is used as the input source. 4.4.2 Synthetic models As represented by WAR and NVI datasets, the characteristic features of the observed wavefield of the Moho below the Great Bear magmatic arc are: (1) prominent PmP phase between 80 and 220 km offsets followed by ~0.5 s long coda, and the absence of this phase at near-vertical-incidence angles (Fig. 4.4); (2) amplitude behavior of the PmP phase with offset is highly erratic and no common amplitude peak is observed among the different shots (Fig. 4.5a); (3) variable dominant PmP frequency, measured on individual traces, at different offsets between 10 and 13 Hz. (However, no clear relation between amplitude spectra and offset is observed; Fig. 4.5B); (4) a very weak amplitude Pn phase observed at offsets > 160 km offset; and (5) a Moho signature on the NVI data is characterized by a very distinct, mutlicyclic  band of reflectors that are laterally variable and piecewise continuous over several kilometres (Fig. 4.3). Amplitude spectra measured over 25 traces show a dominant frequency centered at ~21 Hz. We calculated the synthetic seismograms of 1- and 2-D velocity models to determine the extent to which these models can replicate the wavefield characteristics noted above. 90  Both early studies of the properties of the Moho (e.g., Clowes and Kanasewich, 1970) and more recent ones (e.g., Carbonell et al., 2002) have demonstrated that representation of the crust-mantle transition by a first-order discontinuity, or a linear increase in velocity from lower crustal to upper mantle values over a limited depth interval, does not generate synthetic seismograms that replicate observed wavefield characteristics. Our studies of such Moho models, not discussed herein, further confirmed these results. Accordingly, we initiate our discussion with layered transition models. 4.4.2.1 Layered transition The crust-mantle transition is represented by alternating high and low velocity layers with a velocity contrast of 0.6 km/s (Fig. 4.6). Alternative models with different/variable velocity contrasts (0.2–0.8 km/s) were also examined. Following Carbonell et al. (2002), these layered models are characterized by two parameters, the total thickness of the transition zone (H) and the thickness of the internal layers (h). Initial modeling suggests that (H) controls the length of PmP coda whereas the frequency content is relatively controlled by (h), similar to the results obtained by Carbonell et al. (2002). We kept the total thickness of the transition zone fixed at H = 3 km, which produces an ~ 0.5 s coda length, comparable to the observed data, and we varied the internal thickness (h) between 100 and 500 m. Layered models produce prominent, reverberatory PmP phases beyond 80 km offset but rather weak Pn phases at around 140 km offset (Fig. 4.6). In a test of two models with the same H = 3 km, but one featuring h = 300 m and the other with h = 500 m, the latter generates a phase characterized by a more reverberatory pattern (Fig. 4.7). For very thin internal layers (h ≤ 100 m), the Moho signature is very similar to that generated by a first-order discontinuity model. From models with varying velocity contrasts we note that smaller values produce fewer numbers of wiggles. Amplitude-versus-offset (AVO) curves are relatively simple compared to the observed PmP amplitudes, regardless of the internal thickness of the layers. They are characterized by a clear peak amplitude located at ~120 km offset that corresponds, as expected, to the critical distance (Fig. 4.8). Beyond 120 km there is a nearly steady decrease in amplitudes with offset.  91  Figure 4.6. Synthetic and observed seismograms for a layered transition. (a) Layered velocity model (H = 3 km, h = 500 m) for P-wave and corresponding densities (b) Observed PmP and Pn phases (left) and corresponding synthetic seismogram (right) for SP 1105. All seismograms are displayed in true relative amplitudes and band-pass filtered between 3 - 15 Hz.  The amplitude spectrum of the PmP phase is calculated on single traces and over a 0.5 s window (here and in all subsequent models). The spectrum features a varying dominant frequency between 7 and 12 Hz at different offsets. There is no correlation between the dominant frequency and varying offsets (Fig. 4.9a). Also, the internal thickness of the layers has a limited effect on the shape of the spectrum (Fig. 4.9b). The synthetic NVI seismograms feature a distinct Moho characterized by sharp, mutlicyclic bands of reflectors (4 to 9 bands) depending on the thickness (h) of the internal layers (Fig. 4.10). Higher velocity contrasts between the internal layers result in higher amplitude reflections and vice versa. The dominant frequency of these reflections, measured from averaging the amplitude spectrum over 25 traces and for different layer  92  thicknesses, ranges between 19 and 23 Hz, compared with the 20 Hz central frequency of the source wavelet. This is due to the fine tuning effect of the varying thicknesses of the internal layers (Fig. 4.11). We notice that thicker layers produce more spiked-shape spectra, whereas layers <100 m thick produce a smoother spectrum similar to that of the input source.  Figure 4.7. Two synthetic seismograms resulting from a layered velocity model with the same total thickness, H = 3 km, but a different internal layer thickness h. (a) Internal layer thickness h = 300 m (b) Internal layer thickness h = 500 m. Notice the latter produces more reverberations in the PmP phase.  Previous results and those from our study indicate that the crust-mantle transition is not a simple, uniform discontinuity in which velocities increase suddenly or linearly from crustal to mantle values, but rather a complex and variable transition zone. One possible scenario that more closely replicates the observed data is a transition zone that consists of alternating high and low velocity layers, as described in the preceding based on 1-D modeling. However, a 1-D model does not generate the lateral variability that is observed, particularly on the NVI reflection gathers and stacked sections (Figs. 4.2 and 4.3). In order to account for all the characteristics of the wavefield, a 2-D wave propagation algorithm is needed, one that can handle both vertical and lateral velocity variations.  93  Figure 4.8. Comparison between synthetic (solid line) and observed (dotted line) amplitude vs. offset curves of the PmP phase. Synthetic amplitudes are for a layered transition model with h = 300 m. Amplitude values are normalized to unity.  Figure 4.9. Normalized amplitude spectra of synthetic PmP phase produced by layered discontinuity (h = 300 m). (a) The spectrum is calculated on single traces at 100, 120 and 140 km offsets and over a 0.5 s window. (b) The amplitude spectrum is calculated at 120 km offset on single traces over a 0.5 s window and for different layer thicknesses h = 100, 300 and 500 m.  94  Figure 4.10. (a) Synthetic near-vertical reflection data for a layered transition model with h = 300 m and (b) for a layered transition model with h = 500 m.  Figure 4.11. Normalized amplitude spectra of synthetic near-vertical data produced by layered discontinuity. Each spectrum is averaged over 25 traces. Complexity of amplitude spectra increases with increasing thickness of the layers. Layer thicknesses are 50, 200, 300 and 500 m.  4.4.2.2 Heterogeneous transitions To study the effects of both lateral and vertical heterogeneity at the crust-mantle boundary, we investigate a range of velocity models that represent two end members. In one model we use laterally discontinuous layers with variable thicknesses and velocities. This represents an anastomosing layered structure with juxtaposed rocks from the lower  95  crust and upper mantle (e.g., Jarchow and Thompson, 1989). In the second approach, the velocity model is assumed to be made up of two components (Mereu and Ojo, 1981): (1) a deterministic background component in which the velocity is increasing with depth and (2) a random velocity component comprising a field of small elliptical reflectors embedded in the background velocity field. Random numbers are generated, which in turn are used to distribute the ellipses throughout the transition zone and to assign velocity contrast values to each ellipse (Mereu, 2003) with 2% standard deviation. The ellipses are characterized by their correlation length, which has a horizontal ( a x ) and vertical ( a z ) component. For the WAR data we calculate the synthetic seismograms for 3 shots that simulate shots 1104, 1105 and 1108. In the case of NVI data, a number of synthetic seismograms simulating normal shot gathers are generated. For all cases, the velocity structure on either side of the region of the crust-mantle transition is derived from the refraction analysis of Line 1 data (Fernandez Viejo and Clowes, 2003). Laterally discontinuous layers The Moho transition consists of a complexly layered velocity structure over a 3-km-thick zone. Individual layers have thicknesses varying between 100 and 500 m with their horizontal extent ranging between 5 and 20 km and varying P-wave velocities of 6.9, 7.5 and 8.1 km/s (with 2% standard deviation), characteristic of lower crust and upper mantle values (Fig. 4.12a). Synthetic seismograms generated by the viscoelastic finite-difference algorithm display a prominent PmP phase followed by an ~0.5 s long coda, in agreement with the observed data (Figs. 4.12b and 4.4). Through a variety of tests, layers as short as 3 km long still produce a PmP phase with similar dynamic characteristics. On the other hand, increasing the total thickness of the transition zone to 6 km results in a much longer coda (~1.2 s), which is not observed. In general terms, the synthetic seismograms of Figure 4.12 compare favorably with the observed data of Figure 4.4, including variability of the PmP phase and a weak Pn phase beginning at ~170 km offset. The PmP amplitude analyses show complex behavior with offset in all 3 shots (Fig. 4.13). This resembles the amplitude analysis results for the observed data (Fig. 4.5). Synthetic amplitudes were  96  corrected for geometrical spreading by multiplying all traces by r½ (r = source-receiver distance), assuming a point source in 2D as equivalent to a line source in 3D.  Figure 4.12. Synthetic wide-angle seismograms for laterally discontinuous layers. (a) 2-D velocity model characterized by discontinuous layers of length from 5 to 20 km with thicknesses ranging between 100 and 300 m and alternating velocities of 6.9, 7.5 and 8.1 km/s. (b) Synthetic seismograms simulating SP 1105 and SP 1108. Phases Pg, PmP, and Pn are explained in the text. PiP is a wide-angle phase from a velocity contrast at ~23 km. All seismograms are displayed in true relative amplitudes, and were band-pass filtered between 3-15 Hz. Compare with Fig. 4.4.  97  On the synthetic NVI data, the Moho signature features a strong, mutlicyclic band of reflectors. However, only models comprising layers with shorter horizontal extent (3 to 7 km) display laterally variable and piecewise continuous reflections, similar to the observed data (Fig 4.14).  Figure 4.13. Amplitude vs. offset of the synthetic PmP phase on shots 1104, 1105 and 1108. Amplitude values are normalized to unity. Amplitude shows erratic behavior and no common peak is observed among the different shots.  Lamellae structure The Moho is represented by a 3-km-thick transition zone characterized by randomly distributed ellipses with their correlation lengths a x = 1000 m and a z = 100 m and velocities ranging between 7.1 – 7.8 km/s (Fig. 4.15a). This velocity structure produces a strong and reverberatory PmP phase characterized by ~0.5 s long coda, comparable to the observed data (Fig. 4.15b). Increasing the dimensions of these ellipses ( a x = 3000 m and a z = 300 m) or replacing their varying velocities with one fixed value at 7.2 km/s does  98  not change the dynamic characteristic of the PmP phase substantially. On the other hand, when a x < 700 m, the synthetic seismograms feature a simple PmP phase lacking a visible coda, similar to that generated by a first-order discontinuity. A velocity contrast as small as 0.1 km/s still produces a reverberatory PmP phase. Distributing the ellipses over a 6-km-thick zone produces a distinct PmP phase with a very long coda (~1.5 s), whereas a 1-km-thick transition zone produces a PmP phase without a coda. Table (4.3) summarizes all the different combination of parameters that were tested. All of the above scenarios generate a weak Pn phase observed at offsets >160 km. As in the previous model, the PmP amplitude curves of all the different shots and for different correlation lengths vary considerably with offset (Fig. 4.16).  Figure 4.14. Synthetic near-vertical shot gather for laterally discontinuous layers. Inset shows Moho reflections that are multicyclic, laterally variable and piece-wise continuous over several kilometers.  On the NVI data, the Moho resulting from this lamellae structure is imaged as a narrow (0.4 – 0.7 s TWT) band of 4-8 cycles of reflections (Fig. 17), which is comparable to the observed data (Fig. 3). These reflections are laterally variable and piecewise continuous over ~10 km distance. Our results show that increasing the total thickness of the 99  Figure 4.15. Synthetic wide-angle seismograms for lamellae structure transition zone. (a) 2-D velocity model with randomly distributed ellipses (correlation lengths ax = 1000 m and az= 100 m) and velocities ranging between 7.1 and 7.8 km/s. (b) Synthetic seismograms simulating SP 1104 and SP 1105. Phases Pg, PmP, and Pn are explained in the text. PiP is a wide-angle phase from a velocity contrast at ~23 km. All seismograms are displayed in true relative amplitudes and were band-pass filtered between 3-15 Hz.  100  transition zone from 3 to 6 km would result in a much longer (~0.9 s TWT) band of 1213 cycles of reflections. Different correlation lengths change the number of cycles and even with values as small as a x = 300 m and a z = 30 m, distinct bands of reflections can still be detected, attesting to the higher resolution of the NVI data. However, a minimum velocity contrast of 0.2 km/s is needed between the ellipses and the surrounding material in order to produce visible reflections.  Velocity (km/s)  Correlation length (km)  Thickness (km)  7.1 - 7.8  ax=3, az=0.3  3  7.1 - 7.8  ax=1, az=0.1  3  7.2  ax=3, az=0.3  3  7.2  ax=1, az=0.1  3  8.0  ax=3, az=0.3  3  8.0  ax=1, az=0.1  3  7.1 - 7.8  ax=3, az=0.3  6  7.1 - 7.8  ax=1, az=0.1  6  7.2  ax=3, az=0.3  6  7.2  ax=1, az=0.1  6  7.1 - 7.8  ax=0.5, az=0.1  3  7.2  ax=0.5, az=0.1  3  7.1 - 7.8  ax=3, az=0.3  1  7.2  ax=3, az=0.3  1  PmP characteristics  Strong and reverberatory phase, coda length = ~0.5 s  Strong and reverberatory phase, coda length = ~1.5 s  Strong phase lacking visible coda  Table 4.3. Different combination of all the parameters characterizing the lamellae transition zone.  101  Figure 4.16. Amplitude vs. offset of the synthetic PmP arrivals on shots 1104, 1105 and 1108 produced by the lamellae structure. Amplitude values are normalized to unity. (a) ax = 1000 m, az = 100 m. (b) ax = 3000 m, az = 300 m.  102  Figure 4.17. Synthetic near-vertical seismograms for a lamellae structure transition zone. (a) Correlation length ax = 1000 m and az = 100 m. (b) ax = 3000 and az = 300 m. (c) ax = 300 m and az = 30 m.  103  4.5 Discussion 4.5.1 Synthetic models From the results of synthetic modeling and based on the duration of the PmP coda, a 3km-thick Moho is required to generate 0.5 s long coda. Thicknesses greater than 5 km or less than 2 km produce either a longer coda (~1.0 s) or a PmP phase lacking any visible coda, respectively. Furthermore observed NVI data at Moho levels feature a distinct band of reflections over 0.45 s travel time (0.9 s TWT). Considering average seismic velocity at this level to be ~7.4 km/s, this band translates to ~3.3-km-thick zone, which is in good agreement with the results from the synthetic WAR data. Only a laterally and vertically heterogeneous crust-mantle boundary can simulate all the dynamic features observed on the wide-angle and near-vertical data (Figs 4.12 and 4.15). This transition consists of discontinuous layers/ellipses characterized by variable correlation lengths and velocities. Synthetic models for both WAR and NVI data allow us to place lower and upper limits on the values for such correlation lengths and velocities. A minimum thickness of 100 m and a 0.2 km/s velocity contrast are needed for the layers/ellipses in order to reproduce the reverberatory pattern of the PmP phase. Models with layers/ellipses thicker than 600 m do not simulate the multicyclic band of reflections observed on the NVI data. On the other hand, layers/ellipses that are at least 1 km but no more than 6 to 7 km long are required to generate a distinct PmP coda and to achieve the lateral variability seen on the NVI data. Calculated amplitude versus offset curves for the PmP phase from the layered 1-D velocity model are characterized by one prominent peak amplitude at around 120 km. In contrast, amplitude curves calculated from the heterogeneous models have more complicated shapes with multiple high amplitudes in agreement with the calculated amplitudes from the observed data (Figs. 4.8, 4.13 and 4.16). The deviation between these curves and the ones calculated for simple 1-D velocity models is most probably the result of constructive and destructive interference effects due to the lateral heterogeneity as well as the vertical complexity of the crust-mantle boundary.  104  Based on the spectral analysis, the dominant frequency of the PmP phase varies randomly between 7 and 12 Hz at different offsets. Thus, no relation can be inferred between frequency and offset. In addition, the wide-angle data are less sensitive to the fine structure of the crust-mantle boundary (e.g. internal thickness of the layers) than the near-vertical data (Figs. 4.9b and 4.11). Unfortunately, spectral analysis of the reflection waveforms is insufficient to distinguish between different heterogeneous models. 4.5.2 Multigenetic origin for the continental Moho The Wopmay orogen is an amalgamation of Paleoproterozoic domains that developed and accreted onto the Slave craton between 2.1 Ga and 1.84 Ga by complex geologic processes (Hoffman, 1989). Thus, in this region reflection patterns near the crust-mantle boundary are expected to be similarly complex (e.g., Clowes et al., 1987). However, observations along SNORCLE line 1 reveal simple Moho reflections characterized by subhorizontal surfaces and relatively uniform arrival time (~11.0 s) over the entire profile, even beneath regions with different ages such as the Archean Slave province and the Paleoproterozoic Wopmay orogen (Cook et al., 1999). This uniformity implies that postaccretion, regionally extensive processes have effectively modified rocks and reset the reflection Moho to a uniform depth (~33 km). Therefore, a complete interpretation for the development of the Moho must consider all possible thermal and/or mechanical processes in order to account for its formation and evolution from initial complexities during accretion to increasing coherency over time (Cook, 2002). Several hypotheses have been proposed for the formation and development of the continental Moho, which can be classified under either igneous or non-igneous origins (Eaton, 2006). Igneous origins include thermal front (e.g., partial melt and magmatic intrusion) and magmatic underplating, whereas non-igneous origins include metamorphic or metasomatic front (e.g., mafic granulite to eclogite) and regional décollement (Eaton, 2006; Cook et al., 2010). The spatial correlation between the Moho and crustal reflections provides important constraints on the origin of the crust-mantle boundary. For instance, lower crustal reflections that are listric into the Moho suggest that the crust-mantle boundary is most 105  probably a regional décollement (Cook and Vasudevan, 2003). The transition zone acts as a rheologically weak layer and is underlain by a stronger layer. On the other hand, truncation of lower crustal reflections at Moho levels could be caused by thermal or metamorphic fronts (e.g., Carbonell et al., 2002; Cook et al., 2010). Interestingly, seismic data from the lower crust of the Great Bear magmatic arc, part of SNORCLE line 1, exhibit both types of spatial correlation between the Moho and lower crustal reflections near 10-11 s (Fig. 4.2). Thus, in our synthetic models the lamellae structure with randomly distributed ellipses represents a thermal/metamorphic front, whereas the laterally discontinuous layering represents a regional décollement. 4.5.2.1 Regional décollement Many deep seismic profiles reveal lower-crustal reflections that are listric into the Moho (Cook et al., 1992, 1999; Calvert et al., 1995; Kukkonen et al., 2008). This suggests that the transition zone near the crust-mantle boundary behaved as a mechanically weak layer (more plastic) into which upper crustal detachments soled (Eaton, 2006). This process could be associated with either regional extension and/or contraction (Cook, 2002). Within the lithosphere, the strength of minerals decreases considerably from olivine, amphibole and garnet to feldspar and quartz, the latter being the weakest mineral (Austrheim et al., 1997). The lower crust of the Hottah terrane consists primarily of amphibolite and paragranulite (Fernandez-Viejo et al., 2005). Consequently the rheology is presumably controlled by amphibole and granulite whereas olivine controls the rheology of the mantle. Therefore, at the crust-mantle boundary one expects a distinct strength contrast that may tend to localize detachment surfaces (Cook, 2002; Cook et al., 2010). Based on published geotherms and lithospheric composition, Meissner and Mooney (1998) estimated the depth of low viscosity (weak) zones within the crust and upper mantle along which decoupling may occur. They found these weak zones at three depths: (1) the base of the felsic upper crust; (2) just above the Moho; and (3) some tens of kilometers below the Moho. In the lower crust beneath the Great Bear magmatic arc, a series of discrete, northeastdipping reflections (Fig. 4.2) are listric into the Moho. Moreover, improved seismic 106  images of the same profile (Figure 3d in Cook and Vasudevan, 2003) show distinct upper mantle structures, interpreted as a synform fabric overlain with angular discordance by the flat Moho. Such observations support the interpretation of the Moho as a structural boundary along which lower crustal and upper-mantle rocks were displaced, sheared and blended. However, as Cook and Vasudevan (2003) explained, the Moho here cannot be regarded only as a regional décollement for the following two reasons: (1) the crust is much more reflective than the upper mantle and with a simple detachment the rocks above the detachment surface should be equivalent to, but dislocated from, rocks below; (2) a detachment surface does not explain how rocks below the Moho have significantly higher seismic velocities (~8.1 km/s; Fernandez Viejo and Clowes, 2003) than those above it (~6.8 km/s). From the preceding discussion, consideration of the Moho as a regional décollement represents only one stage of its development. A logical sequence of stages or events is probably necessary to explain the varying characteristics of the crust-mantle boundary observed on both NVI and WAR data. 4.5.2.2 Thermal front (e.g., magmatic intrusion and partial melting) Mafic igneous rocks may be injected into the lower crust and upper mantle in the form of sills during the amalgamation of island arcs at collisional orogens (e.g., Wopmay orogen). These sills are characterized by relatively high density and high seismic velocity (>7.2 km/s). The intermediate composition crust above acts as a density filter for these intrusives and hence they pond within the lower crust and along the crust-mantle boundary, producing seismically laminated lower crust and sharply defined Moho (Furlong and Fountain, 1986; Nelson, 1991). Partial melting can modify the chemical and physical properties of intermediate composition rocks at elevated lower-crustal temperatures and pressures by preferentially removing lower density minerals such as quartz and feldspar. This will leave mafic residuals with higher density and higher seismic velocity (>7.5-7.8 km/s) appropriate for upper mantle (Hynes and Snyder, 1995; Rudnick and Fountain, 1995). The solidus temperature for these restites is considerably elevated and any further partial melting 107  requires substantially higher temperatures. As a result the Moho could be ‘frozen’ unless significantly higher temperatures occur (Cook 2002). Because partial melting might not remove all of the lighter fraction from the lower crust, structural geometry near the Moho can be preserved (Cook et al., 2010). As mentioned previously, termination or truncation of some of the dipping structures in the lower crust of the Great Bear magmatic arc by the Moho, in addition to the discordance between the underlying synform fabric and the Moho (Cook and Vasudevan, 2003), strongly suggest that the Moho postdates the synform and lower crustal structures. Thus, it may have been superimposed on them by thermal processes. The approach and subsequent subduction of oceanic lithosphere beneath the Hottah terrane is probably the source of thermal activity responsible for partial melting and magmatic intrusion in the study area, which in turn is associated with the deposition of the Great Bear magmatic rocks between ~1.88 to 1.84 Ga (Cook et al., 1999; Oueity and Clowes, 2010). 4.5.2.3 Metamorphic front (eclogite phase transition) Since the work of Green and Ringwood (1967) and Ito and Kennedy (1971), it has been recognized that under relatively high P-T conditions, mafic rocks (e.g., mafic granulite, gabbros) might transform into eclogite at depths of ~35 km. Hence the Moho could be interpreted as a metamorphic front associated with a phase change. Eclogitization results in an increase in the density and seismic velocity of the lower crust such that it becomes indiscernible from the uppermost mantle (e.g., Austrheim, 1987; Mengel and Kern, 1992; Fischer, 2002; Kukonnen et al., 2008). Consequently, lower crustal eclogitization has been invoked to explain the observation of a relatively shallow Moho at several places (e.g., the Scottish Caledonides [Hynes and Snyder, 1995]; the Appalachians off northeastern Newfoundland [Chian et al., 1998]). Along the southwestern coast of Norway, eclogites are exposed over a large area in the Western Gneiss Region (WGR). They form pods and lenses ranging in size from decimeters to several 100 meters in length embedded in amphibolites (Austrheim, 1987). This is very similar to the lamellae structure we tested in our synthetic models. Although in the WGR most eclogites are mylonitic, examples of partial eclogitization are abundant. 108  Therefore, in our model the wide range of velocities/densities characterizing the randomly distributed ellipses could be explained by different degrees of eclogitization depending on the rock composition and availability of water.  4.6 Conclusions Observed wide-angle and near-vertical seismic reflection data as well as 1- and 2-D synthetic seismogram modeling for the crust-mantle boundary offer different but complementary information about its nature and formation. A unified interpretation, in which both wide-angle and near-vertical data derive from the same Moho location, can provide a more complete explanation concerning the structure of the Moho, one with higher resolution and perhaps less ambiguity. Our synthetic models lead to the conclusion that beneath the Great Bear magmatic arc, an approximately 3 km thick, laterally and vertically heterogeneous crust-mantle boundary properly simulates the dynamic characteristics of the observed wavefields recorded on the near-vertical and wide-angle reflection data. The heterogeneity can be achieved by either laterally discontinuous layering or a lamellae structure with randomly distributed ellipses. These two models may represent different re-equilibration mechanisms, one characterized by thermal/metamorphic fronts and another that features a regional décollement. Due to insufficient details on the reflection seismic image, it is difficult to establish the exact geometric relationships between the Moho and lower crustal layers and consequently the re-equilibration process. However, assuming that our images of the Moho beneath the Great Bear magmatic arc represent ‘snapshots’ of the crust-mantle boundary ‘taken’ at different stages of its evolution (Keller et al., 2005), we believe that the Moho has formed as a combination of the aforementioned processes. Here we propose a scenario for the evolution of the Moho (Fig. 4.18) that builds on, but differs somewhat from a previous study (Cook and Vasudevan, 2003) as follows: (a) Pre-1.9 Ga the Moho represented the base of the Hottah terrane crust, which was generated as a magmatic arc on cryptic 2.4 to 2.0 Ga crust (Hoffman and Bowring, 1984). The depth to the Moho is unknown but was perhaps at about 35 km. (b) The collision of  109  Figure 4.18. Schematic evolutionary model for the development of the Moho visible beneath the Great Bear arc along SNORCLE Line 1. (A) Formation of the Hottah terrane as a magmatic arc west of the Slave craton. Pre- 1.9 Ga tectonism may have produced the internal structure within the lower crust. (B) Accretion of the Hottah terrane onto the Slave craton as a result of west-dipping subduction. The collision caused folding and faulting within the crust and the Moho acted as a detachment surface. This has probably led to the formation of the listric fabric in the lower crust. (C) Formation and deposition of the Great Bear magmatic arc on both the Hottah terrane and the Coronation Supergroup as a result of west-dipping subduction. The discordance between the arc rocks and the fold structures below entails that the crustal deformation, and consequently Moho detachment, occurred prior to arc magmatism (Cook and Vasudevan, 2003). We suggest that partial melt and magmatic intrusives at the base of the crust, produced by this thermal activity, are responsible for the Moho development as a thermal/metamorphic front.  110  the Hottah terrane onto the Slave craton caused substantial compression and led to the formation of the listric structures that flatten into the Moho. This implies that the Moho at this stage acted as a ductile decoupling zone or a detachment surface. (c) Between 1.88 and 1.84 Ga, the Great Bear magmatic arc was formed as a result of east-dipping subduction of oceanic lithosphere beneath the Hottah and the Slave craton generated by incipient of the Fort Simpson terrane. The thermal activity associated with the previous subduction produced partial melt and magmatic intrusions at the base of the crust. A horizontal Moho was developed as a thermal/metamorphic boundary.  111  Bibliography Abbott, G., Thorkelson, D., Creaser, R., Bevier, M.L., and Mortensen, J., 1997. New correlations among Proterozoic successions and intrusive breccias in the Ogilvie and Wernecke Mountains, Yukon: in SNORCLE Transect and Cordilleran Tectonics Workshop Meeting, University of Calgary, Calgary, Alta., March 7-9, 1997. Compiled by F.A. Cook, and P. Erdmer. LITHOPROBE Secretariat, The University of British Columbia, Vancouver, B.C. Lithoprobe Report No. 56, p. 188–197. Austrheim, H., Erambert, M., and Engvik, A. K., 1997. Processing of crust in the root of the Caledonian continental collision zone: the role of eclogitization. Tectonophysics, 273, p.129– 153. Austrheim, H., 1987. Eclogitization of lower crustal granulites by fluid migration through shear zones. Lithos, 27, p. 29–42. Braile, L.W., and Chiang, C. S., 1986, The continental Mohorovičić discontinuity results from near-vertical and wide-angle seismic reflection studies. In Barazangi, M., Brown, L., eds., Reflection Seismology: A Global Perspective, Washington DC., American Geophysical Union Geodynamic Series, 13, p. 257-272. Boudier, F., and Nicolas, A., 1995. Nature of the Moho transition zone in the Oman Ophiolite. Journal of Petrology, 36, p. 777-796. Calvert, A.J., Sawyer, E.W., Davis, W.J., and Ludden, J.N., 1995. Archean subduction inferred from seismic images of a mantle suture in the Superior Province. Nature, 375, p. 670–674. Candés, E., Demanet, L., Donoho, D., and Ying, L., 2006. Fast discrete curvelet transforms. Multiscale Modeling and Simulation, 5, p. 861–899. Carbonell, R., Gallart, J., and Perez-Estaun, A., 2002. Modeling and imaging the Moho transition: the case of the southern Urals. Geophysical Journal International, 149, p. 134-148. Carpentier, S. and Roy-Chowdhury, K., 2007. Underestimation of scale lengths in stochastic fields and their seismic response: a quantification exercise. Geophysical Journal International, 169, p. 547–562.  112  Casey, J. F., Dewey, J. F., Fox, P. J., Karson, J. A., and Rosencrantz, E., 1981. Heterogeneous nature of the oceanic crust and upper mantle: a perspective from the Bay of Islands Ophiolite Complex, in Emiliani, C., ed., The Sea: v 7. New York John Wiley and Sons, p. 305-338. Chian, D., Marillier, F., Hall, J., and Quinlan, G., 1998. An improved velocity model for the crust and upper mantle along the central mobile belt of the Newfoundland Appalachian orogen and its  offshore  extension.  Canadian  Journal  of  Earth  Sciences,  35,  p.  1238–  1251, doi:10.1139/cjes-35-11-1238. Christensen, N. I., and Mooney, W. D., 1995. Seismic velocity structure and composition of the continental crust: A global view. Journal of Geophysical Research, 100, p. 9761–9788. Christensen, N. I., 1996. Poisson's ratio and crustal seismology. Journal of Geophysical Research, 101, p. 3139–3156. Clowes, R. M., and Kanasewich, E. R., 1970. Seismic attenuation and the nature of reflecting horizons within the crust. Journal of Geophysical Research, 75, p. 6693-6705. Clowes, R.M., Brandon, M.T., Green, A.G., Yorath, C.J., Sutherland-Brown, A., Kanasewich, E.R., and Spencer, C., 1987. Lithoprobe Southern Vancouver Island: Cenozoic subduction complex imaged by deep seismic reflections. Canadian Journal of Earth Sciences, 24, p. 31– 51. Coles, R.L., Haines, G.V., and Hannaford, W., 1976. Large-scale magnetic anomalies over western Canada and the Arctic: a discussion. Canadian Journal of Earth Sciences, 13, p. 790– 802. Cook, F.A., Varsek, J.L., Clowes, R.M., Kanasewich, E.R., Spencer, C.S., Parrish, R.R., Brown, R.L., Carr, S.D., Johnson, B.J., and Price, R.A., 1992. Lithoprobe crustal reflection cross section of the southern Canadian Cordillera, I: Foreland thrust and fold belt to Fraser River fault. Tectonics, 11, p. 12–35. Cook, F.A., van der Velden, A.J., and Hall, K.W., 1999. Frozen subduction in Canada's Northwest Territories: LITHOPROBE deep lithospheric reflection profiling of the western Canadian Shield. Tectonics, 18, p. 1-24. Cook, F.A., 2002. Fine structure of the continental reflection Moho. Geological Society of America Bulletin, 114, p. 64-79.  113  Cook, F. A., and Vasudevan, K., 2003. Are there relict crustal fragments beneath the Moho?. Tectonics, 22, p. 1026, doi:10.1029/2001TC001341. Cook, F. A., and Erdmer, P., 2005. An 1800 km cross section of the lithosphere through the northwestern North American plate: lessons from 4.0 billion years of Earth’s history. Canadian Journal of Earth Sciences, 42, p. 1295–1311, doi: 10.1139/E04-106. Cook, F. A., White, D. J., Jones, A. G., Eaton, D. W., Hall, J., and Clowes, R. M., 2010. How the crust meets the mantle: Lithoprobe perspectives on the Mohorovičić discontinuity and crustmantle transition. Canadian Journal of Earth Sciences, 47, p. 315-351. DeBari, S., Kay, S. M., Kay, R. W., 1987. Ultramafic xenoliths from Adagdak Volcano, Adak, Aleutian Islands, Alaska: deformed igneous cumulates from the Moho of an island arc. Journal of Geology, 95, p. 329-341. Dilek, Y., Furnes, H., and Shallo, M., 2008. Geochemistry of the Jurassic Mirdita Ophiolite (Albania) and the MORB to SSZ evolution of a marginal basin oceanic crust. Lithos, 100, p. 174-209, doi: 10.1016/j.lithos.2007.06.026. Eaton, D. W., 2006, Multi-genetic origin of the continental Moho: insights from LITHOPROBE: Terra Nova, 18, p. 34-43, doi: 10.1111/j.1365-3121.2005.00657. Elad, M. J., Starck, L., Querre, P., and Donoho, D. L., 2005. Simultaneous cartoon and texture image in painting using morphological component analysis (MCA). Applied and Computational Harmonic Analysis, 19, p. 340–358. Ellis, R., Clowes, R.M., Spence, G.D., Asudeh, I., Hajnal, Z., and Amor, J.R., 2002. SlaveNorthern Cordillera Refraction Experiment (SNORE’97). Field Acquisition and Preliminary Data Processing Report: University of British Columbia, LITHOPROBE Report No. 83, 14 pp, 10 apps. Fernandez Viejo, G., Clowes, R.M., and Amor, J.R., 1999. Imaging the lithospheric mantle in northwestern Canada with seismic wide-angle reflections. Geophysical Research Letters, 26, p. 2809-2812. Fernandez Viejo, G., and Clowes, R.M., 2003. Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay orogen, northwestern Canada, from a LITHOPROBE refraction/wide-angle reflection survey. Geophysical Journal International, 153, p. 1-19.  114  Fernandez Viejo, G., Clowes, R.M., and Welford, J. K., 2005. Constraints on the composition of the crust and uppermost mantle in northwestern Canada: Vp/Vs variations along Lithoprobe's SNorCLE transect. Canadian Journal of Earth Sciences, 42, p. 1205-1222, doi:10.1139/e05028. Fischer, K.M., 2002. Waning buoyancy in the crustal roots of old mountains. Nature, 417, p. 933–935. Fountain, D. M., and Salisbury, M. H., 1981. Exposed cross-sections through the continental crust: implications for crustal structure, petrology and evolution. Earth and Planetary Sciences Letters, 56, p. 263-277. Fountain, D. M., 1986. Implications of deep crustal evolution for seismic reflection interpretations, in Barazangi, M., Brown, L., eds., Reflection Seismology: A global Perspective: Washington DC, American Geophysical Union Geodynamic Series, 14, p. 95106. Fowler, C. M. R., 2005. The Solid Earth: An Introduction to Global Geophysics. Cambridge, UK, Cambridge University Press, 685 p. Fuchs, K ., and Müller, G., 1971. Computation of synthetic seismograms with the reflectivity method and comparison with observations. Geophysical Journal of Royal Astronomical Society, 23, p. 417-433. Furlong, K.P., and Fountain, D.M., 1986. Continental crust underplating: thermal considerations and seismic-petrologic consequences. Journal of Geophysical Research, 91, p. 8285–8294. Gandhi, S.S., Mortensen, J.K., Prasad, N., and van Breemen, O., 2001. Magmatic evolution of the southern Great Bear continental arc, northwestern Canadian Shield: geochronological constraints. Canadian Journal of Earth Sciences, 38, p. 767-785. Gibbs, A. K., 1986. Seismic reflection profiles of Precambrian crust: a qualitative assessment, in Barazangi, M., Brown, L., eds., Reflection Seismology: A global Perspective: Washington DC, American Geophysical Union Geodynamic Series, 14, p. 95-106, Green, D.H., and Ringwood, A.E., 1967. An experimental investigation of the gabbro to eclogite transition and its petrological applications. Geochimica et Cosmochimica Acta, 31, p. 767833.  115  Hale, L. D., and Thompson, G. A., 1982. The seismic reflection character of the continental Mohorovičić discontinuity. Journal of Geophysical Research, 87, p. 4625-4635. Hammer, P. T. C., and Clowes, R.M., 1997. Moho reflectivity patterns - a comparison of Canadian LITHOPROBE transects. Tectonophysics, 269, p. 179-198. Hasselgren, E. O., and Clowes, R. M., 1995. Crustal structure of northern Juan de Fuca plate from multichannel reflection data. Journal of Geophysical Research, 100, p. 6469-6486. Hennenfent, G. and Herrmann, F.J., 2006. Seismic denoising with nonuniformly sampled curvelets. Computing in Science and Engineering, 8, p. 16-25. Herrmann, F. J., and Hennenfent, G., 2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International, 173, p. 233–248. Hildebrand, R. S., Hoffman, P.F., and Bowring, S.A., 1987. Tectono-magmatic evolution of the 1.9-Ga Great Bear magmatic zone, Wopmay orogen, northwestern Canada. Journal of Volcanology and Geothermal Research, 32, p. 99-118. Hildebrand, R. S., and Bowring, S. A., 1999. Crustal recycling by slab failure. Geology, 27, p. 11-14. Hoffman, P. F., and Bowring, S. A., 1984. Short-lived 1.9 Ga continental margin and its Destruction. Wopmay orogen, NW Canada. Geology, 12, p. 68-72. Hoffman, P.F., 1988. United plates of America, the birth of a craton: Early Proterozoic assembly and growth of Laurentia. Annual Reviews of Earth and Planetary Sciences, 16, p. 543-603. Hoffman, P.F., 1989. Precambrian geology and tectonic history of North America. In Bally, A. W., and  Palmer, A. R., eds., The Geology of North America: An Overview, Geological  Society of America, Boulder , Colorado v. A, p. 447-512. Hurich, C., 2003. The Nature of Crustal Seismic Heterogeneity: A Case Study From the Grenville Province. In Goff, J. A. and Holliger, K., eds., Heterogeneity in the crust and upper mantle: nature, scaling, and seismic properties, Kluwer Academic Publishers, New York, p. 299–320. Hynes, A., and Snyder, D.B., 1995. Deep-crustal assemblages and potential for crustal rocks below the Moho in the Scottish Caledonides. Geophysical Journal International, 123, p. 323– 339.  116  Ito, K., and Kennedy, G.C., 1971. An experimental study of the basalt-garnet granulite-eclogite transition. In Heacock, J.G., ed., The Earth beneath the continents: American Geophysical Union Geophysical Monograph 14, p. 303–314. Jarchow, C.M., and Thompson, G.A., 1989. The nature of the Mohorovicic discontinuity. Annual Review of Earth and Planetary Sciences, 17, p. 475-506. Jousselin, D., and Nicolas, A., 2000. The Moho transition zone in the Oman Ophiolite; relation with wehrlites in the crust and dunites in the mantle. Marine Geophysical Researches, 21, p. 229-241. Karson, J. A., Collins, J. A., and Casey, J. F., 1984. Geologic and seismic velocity structure of the crust/mantle transition in the Bay of Islands Ophiolite Complex. Journal of Geophysical Research, 89, p. 6126-6138. Keller, G. R., Karlstrom, K. E., Williams, M. L., Miller, K. C., Andronicos, C., Levander, A.R., Snelson, C. M., Prodehl, C., 2005. The dynamic nature of the continental crust-mantle boundary: Crustal evolution in the southern rocky mountain as an example. In Karlstrom, K. E., and Keller, G. R., eds., The Rocky Mountain Region: An Evolving Lithosphere, Washington, American Geophysical Union Geophysical Monograph Series 154, , p. 403-420. Kopylova, M. G., Russell, J. K., and Cookenboo, H., 1998. Upper-mantle stratigraphy of the Slave craton, Canada: Insights into a new kimberlites province. Geology, 26, p. 315-318. Kukkonen, I.T., Kuusisto, M., Lehtonen, M., Peltonen, P., 2008. Delamination of eclogitized lower crust: Control on the crust-mantle boundary in the central Fennoscandian shield. Tectonophysics, 457, p. 111-127, doi:10.1016/j.tecto.2008.04.029. Kumar, V., 2009. Incoherent noise suppression and deconvolution using curvelet-domain sparsity [Master thesis]: Vancouver, University of British Columbia, 58 p. Levander, A. R., and Holliger, K., 1992. Small scale heterogeneity and large scale velocity structure of the continental crust. Journal of Geophysical Research, 97, p. 8792-8804. McGetchin, R. R., and Silver, L. T., 1972. A crustal-upper mantle model for the Colorado Plateau based on observations of crystalline rock fragments in the Moses Rock dike. Journal of Geophysical Research, 77, p. 7022-7037.  117  Meissner, R. 1986. The Continental Crust – A Geophysical Approach. Academic Press, International Geophysics Series, 34, 426 p. Meissner, R., 2000. The mosaic of terranes in central Europe as seen by deep reflection studies. In Mohriak, W. and Talwani, M., eds., Atlantic Rifts and Continental Margins, Washington. American Geophysical Union Geophysical Monograph Series 115, p. 33-55. Meissner, R., and Mooney, W., 1998. Weakness of the lower continental crust: a condition for delamination, uplift, and escape. Tectonophysics, 296, p. 47-60. Mengel, K. & Kern, H., 1992. Evolution of the petrological and seismic Moho-implications for the continental crust-mantle boundary. Terra Nova, 4, p. 109-1 16. Mereu, R. F., 2003. The heterogeneity of the crust and its effect on seismic wide-angle reflection fields. In Goff. J. A., and Holliger, K., eds., Heterogeneity in the crust and upper mantle: nature, scaling, and seismic properties, Kluwer Academic Publishers, New York, p. 299–320. Mereu, R. F., and Ojo, S. B., 1981. The scattering of seismic waves through a crust and upper mantle with random lateral and vertical inhomogeneities. Physics of the Earth and Planetary Interiors, 26, p. 233-240. Mjelde, R., Fjellanger, J. P., Digranes, P., Kodaira, S., Shimamura, H., and Shiobara, H., 1997. Application of the single-bubble airgun technique for OBS-data acquisition across the Jan Mayen Ridge. North Atlantic: Marine Geophysical Researches, 19, p. 81-96. Mohorovičić, A., 1910. Das Beben vom 8. X. 1909, Jahrbuch des meteorologischen Observatoriums in Zagreb, 9/4, p. 1-63. Mooney, W. D., 1987. Seismology of the continental crust and upper mantle. Reviews of Geophysics, 25, p. 1168-1176. Morozov, I. B., Smithson, S. B., Chen, J., and Hollister, L. S., 2001. Generation of new continental crst and terrane accretion in southeastern Alaska and western British Columbia: constraints from P- and S-wave wide-angle seismic data (ACCRETE). Tectonophysics, 341, p. 49-67. Nedimović, M. R., Carbotte, S. M., Harding, A. J., Detrick, R. S., Canales, J. P., Diebold, J. B., Kent, G. M., Tischer, M., and Babcock, J. M., 2005. Frozen magma lenses below the oceanic crust. Nature, 436, p. 1149-1152, doi:10.1038/nature03944.  118  Neelamani, R., Baumstein, A.I., Gillard, D.G., Hadidi, M.T., and Soroka, W.L. , 2008. Coherent and random noise attenuation using the curvelet transform. The Leading Edge, 27, p. 240-248. Nelson, K.D., 1991. A unified view of craton evolution motivated by recent deep seismic reflection and refraction results. Geophysical Journal International, 105, p. 25–35. Nielsen, L., and Thybo, H., 2006. Identification of crustal and upper mantle heterogeneity by modeling of controlled-source seismic data. Tectonophysics, 416, p. 209–228. Oueity, J., and Clowes, R. M., 2010. Paleoproterozoic subduction in northwestern Canada from near-vertical and wide-angle seismic reflection data. Canadian Journal of Earth Sciences, 47, p. 35–52, doi:10.1139/E09-073. Pancea, I., Stephenson, R., Knapp, C., Mocanu, V., Drijkoningen, G., Matenco, L., Knapp, J., Prodehl, K., 2005. Near-vertical seismic reflection image using a novel acquisition technique across the Vrancea Zone and Foscani Basin, south-eastern Carpathians (Romania). Tectonophysics, 410, p. 293-309. Robertsson, J., Balnch, J., and Symes, W., 1994. Viscoelastic finite-difference modeling. Geophysics, 59, p. 1444 1456. Rudnick, R.L., and Fountain, D.M., 1995. Nature and composition of the continental crust: A lower crustal perspective. Reviews of Geophysics, 33, p. 267–310. Scheirer, D. S., and Hobbs, R. W., 1990. Seismic attenuation in the continental crust SW of England. Geophysical Journal International, 103, p. 553-540. Spence, G.D., Whittall, K.P., and Clowes, R.M., 1984. Practical synthetic seismograms for laterally varying media calculated by asymptotic ray theory. Bulletin of Seismological Society of America, 74, p. 1209-1223. Steinhart, J. S., 1967. Mohorovičić discontinuity, In Runcorn, S. K., ed., International Dictionary of Geophysics, Volume: 2, New York, Pergamon Press, p. 991-994. Thompson, R.I., Mercier, E., and Roots, C.F., 1987. Extension and its influence on Canadian Cordilleran passive margin evolution. In Coward, M. P., Dewey, J. F., and Hancock, P. L., eds., Continental Extension Tectonics, Geological Society of London, Special Publication, p. 409-417 Tittgemeyer, M., Wenzel, F., Ryberg, T., and Fuchs, K., 1999. Scales of heterogeneities in the continental crust and upper mantle. Pure and Applied Geophysics, 156, p. 29-52. 119  Tittgemeyer, M., Ryberg, T., Wenzel, F., and Fuchs, K., 2003. Heterogeneities of the uppermost mantle inferred from controlled-source seismology. In Goff, J. A., and Holliger, K., eds., Heterogeneity in the crust and upper mantle: nature, scaling, and seismic properties, Kluwer Academic Publishers, New York, p. 299–320. Villieneuve, M.E., Theriault, R.J., and Ross, G.M., 1991. U-Pb ages and Sm-Nd signature of two subsurface granites from the Fort Simpson magnetic high, northwest Canada. Canadian Journal of Earth Sciences, 28, p. 1003-1008. Welford, J. K., Clowes, R. M., Ellis, R. M., Spence, G. D., Asudeh, I., and Hajnal, Z., 2001. Lithospheric structure across the craton-Cordilleran transition of northeastern British Columbia. Canadian Journal of Earth Sciences, 38, p. 1169–1189, doi:10.1139/cjes-38-81169. Yilmaz, O., 2001. Seismic data analysis: Tulsa, Oklahoma,Society of Exploration Geophysicists, 2024 p. Zelt, C.A., and Smith, R.B., 1992. Seismic travel-time inversion for 2-D crustal velocity structure. Geophysical Journal International, 108, p. 16-34.  120  Chapter 5  Heterogeneity of the lower crust, Moho and upper mantle from statistical analysis of near-vertical and wide-angle seismic reflection data4 5.1 Introduction Traditionally, the assumption that the Earth’s physical structure is layered and that fluctuations of the physical parameters within individual layers are smooth in structure and small in magnitude was used to facilitate the interpretation of seismic and other geophysical data. However, such an oversimplified geophysical model has resulted in some source of disagreement among the various disciplines of Earth sciences. From a petrological perspective, the Moho represents either a phase transition from mafic granulites (lower crust) to eclogite (upper mantle) (Ito and Kennedy, 1971) or a compositional transition from mafic granulites to peridotite (Green and Ringwood, 1972). On the other hand, the seismic refraction Moho is interpreted as the depth at which the Pwave velocity exceeds 7.6 km/s (Steinhart, 1967) which may not be spatially coincident with the petrologic Moho. Exposures of crystalline rocks that probably represent the crust-mantle transition, such as the Ivrea Zone in northern Italy (Fountain, 1976), exhibit a high degree of irregularity in lithology and spatial complexity, part of which reflects the structures at depth, and part of which results from the exhumation process itself. The spatial distribution of this 4  A version of this chapter will be submitted for publication.: Oueity, J., Carpentier, S. F. A., and Clowes,  R. M. Heterogeneity of the lower crust, Moho and upper mantle from statistical analysis of near-vertical and wide-angle seismic reflection data.  121  heterogeneity is a function of the magmatic and tectonic processes that progressively redistribute the various lithologic components of the crust (Fountain and Salisbury, 1981). Unlike sedimentary structures, which are often repetitive in depth with specific periods and scales, crustal heterogeneity exists over a wide range of scales and has a fractal-like geometry due to the characteristic ratio among the various scales over which the heterogeneity occurs (Mandelbrot, 1983). Although multiscale fractals can have a systematic distribution, most of the crustal heterogeneities display random distribution. Thus, one can describe crustal heterogeneity as a stochastic distribution that unites both terms, ‘fractal’ and ‘random’ (Carpentier, 2007). Stochastic heterogeneity can also be defined as the shorter-scale fluctuations in the physical properties of the rocks which in turn are superimposed on longer-scaled deterministic variations in average properties. Digitization and statistical analysis of fine-scaled geologic maps of exhumed crystalline rocks reveal that the stochastic distributions of heterogeneity exhibit spatially self-similar or power-law scaling, often described by the von Karman autocorrelation function (Goff and Jordan, 1988; Holliger and Levander, 1992; Holliger, 1996; Dolan et al., 1998; Goff and Jennings, 1999). Such analyses have been done on a number of exposed lower crustal sections including the Ivrea Zone in northern Italy (Holliger and Levander, 1992), the Lewisian gneissic complex in Scotland (Levander et al., 1994), the Franciscan complex in the Diablo Range of California (Goff and Levander, 1996) and the Hafafit complex in eastern Egypt (Goff and Levander, 1996). The von Karman function is parameterized by 2 fractal parameters: correlation length and Hurst number. Nyquist wavelength, dependent on the sampling interval of an autocorrelated sequence, and correlation length represent the lower and upper limits of the von Karman power spectrum respectively. Correlation length is an upper limit for the scale invariance in the heterogeneity, where scale-invariance is defined as the absence of typical scales within a range of observations. Hence, a scale-invariant morphology appears statistically identical at all scales (Marsan and Bean, 2003). The Hurst number is a special decay exponent from fractal mathematics that controls the degree of scale invariance in heterogeneity below the correlation length scale. It is related to the Haussdorf-Besicovitch fractal dimension D = E+1 –v, where E is the Euclidian dimension (2 in this study) and v is the 122  Hurst number (Haussdorf, 1918). In general terms one can think of the Hurst number as a measure of the ‘roughness’ of lithological interfaces between the heterogeneities. A series of studies have provided empirical evidence that the von Karman type distribution of causative scattering heterogeneity is conserved in both synthetic and real near-vertical seismic reflection data (Gibson, 1991; Holliger et al, 1992; Hurich, 1996; Pullammanappallil et al, 1997; Hurich and Kocurko, 2000; Poppeliers, 2007; Carpentier et al., 2009). Because of this direct relation between the reflection data and the heterogeneity in causative geological structure/velocity models, extracted stochastic parameters can be used as constraints on crustal/upper mantle structures. A fundamental underestimation in correlation length and in Hurst number, which was explained by Carpentier and Roy-Chowdhury (2007), was observed in the method. Although the explanation was founded both by analytical and numerical analyses, a unique upscale factor for the underestimated correlation lengths and Hurst number could not be established. Empirical values were derived as working values, but with relatively large uncertainties. Thus far, the stochastic characterization studies have been directed into two separate approaches, the forward problem and the inverse problem. In the former, a stochastic velocity model is constructed based on global statistical distribution of lithologies in exhumed rocks and synthetic seismograms are then calculated and compared to field data (e.g. Holliger and Levander, 1992; Holliger et al., 1993; Tittgemeyer et al., 1999; Ryberg et al., 2000; Nielsen and Thybo, 2006). In the latter, the stochastic power law distribution with the associated fractal parameters are estimated directly from the reflection seismic data and then objective maps of the crustal/upper mantle fabrics are made (Holliger et al., 1994; Pullammanappallil et al., 1997; Hurich 2003; Carpentier et al., 2009; Carpentier et al., 2010). In this study we combine both approaches, the forward and inverse problems, by taking advantage of the coincident, high quality near-vertical and wide-angle reflection data recorded along Line 1 of LITHOPROBE’s Slave-NORthern Cordillera Lithospheric Evolution (SNORCLE) (Cook et al., 1999; Fernandez Viejo and Clowes, 2003). This is  123  achieved in three steps: 1) estimates for correlation length, including their associated maximum-likelihood uncertainties, are taken directly from the near-vertical reflection dataset to yield a range of likely power law parameters for the transition of lower crust into uppermost mantle; 2) a synthesis of these correlation length estimates, Hurst numbers from the literature (Holliger and Levander, 1992) and velocities from wideangle studies (Fernandez Viejo and Clowes, 2003) leads to candidate velocity models for the lower crust and uppermost mantle; and 3) reflection shot gathers and wide-angle reflection data are simulated from these models through visco-elastic full wavefield solutions and compared to observed shot gathers and wide-angle seismograms. For the wide-angle data, we go beyond visual comparison and directly compare attributes of the synthetic recordings (complex trace envelope, envelope energy decay and frequencyspace panels) amongst each other and with the observed data. . Our goal is two-fold: (1) to improve the calibration of the upscaling values for the fractal parameters by using two independent datasets and (2) to provide stochastic signatures for the Moho transition zone, the lower crust and uppermost mantle from the best possible estimates of the fractal parameters . In turn, this may provide additional understanding of the geodynamical processes associated with the formation of that part of the lithosphere (e.g., Oueity and Clowes, 2010).  5.2 The Wopmay orogen The Wopmay orogen is a Paleoproterozoic amalgamation of four principal north-south trending domains (Fig. 5.1) that developed west of the Slave craton between 1.92 Ga and 1.84 Ga (Hoffman, 1989). From east to west, the tectonic domains are: (1) the Coronation Supergroup (2) the Great Bear magmatic arc (3) the Hottah terrane and (4) the Fort Simpson terrane. The Coronation Supergroup represents a westward-facing passive margin, which developed between 1.92 Ga and 1.90 Ga (Hoffman and Bowring, 1984). The Hottah terrane formed as a magmatic arc and was accreted onto the Slave craton as a result of west-dipping subduction of an ocean basin around 1.89–1.88 Ga (Hildebrand and Bowring, 1999). Following the Hottah-Slave collision, the initiation of an eastdipping subduction zone led to the generation of the Great Bear magmatic arc during the 124  period between 1.875 to 1.843 Ga (Hildebrand et al., 1987; Gandhi et al., 2001). The Great Bear arc comprises a 100- by 1,000-km zone of calcalkaline, mainly intermediate, volcanic and plutonic rocks developed primarily on top of the Hottah terrain (Hoffman, 1989). Cessation of the arc magmatism may coincide with the collision of the Fort Simpson exotic terrane with the western margin of the Hottah terrane some time between 1.844 Ga and 1.710 Ga (Hoffman and Bowring, 1984). Along SNORCLE Line 1, Phanerozoic sedimentary rocks overlie the Proterozoic rocks of the Wopmay orogen with increasing thickness from the feather edge at the eastern boundary of the Great Bear magmatic arc to about 1000 m over the Nahanni domain (Cook et al., 1999). The Nahanni domain, west of the Fort Simpson terrane, is an enigmatic feature because it does not outcrop and is only defined by its magnetic signature (Villeneuve et al., 1991).  Figure 5.1. Tectonic map of the study area and location of SNORCLE near-vertical incidence reflection and refraction/wide-angle reflection Line 1 (white segment is the focus of the study). Stars identify shot locations; some are identified by numbers. Sediments of the Phanerozoic Western Canada Sedimentary Basin overly the Precambrian domains west of the long dashed line. Short dashed black lines show political boundaries. CS - Coronation Supergroup, SD - Sleepy Dragon complex, GSLsz - Great Slave Lake shear zone, AB - Alberta, BC - British Columbia, NWT - Northwest Territories, YK - Yukon. Inset shows location of map within Canada.  125  5.3 Seismic data Two types of active-source seismic datasets were recorded along SNORCLE Line 1, near-vertical and refraction/wide-angle reflection data. Line 1 is 725 km long and crosses the southwest Slave craton and the Wopmay orogen (Fig. 5.1). The reflection experiment comprised a total of 404 geophones per shot with a 60-m receiver spacing, resulting in a stacked trace every 30 m. Shot points were occupied every 90 m providing a nominal stacking fold of 134 with sweep frequencies ranging between 10 and 80 Hz (Cook et al., 1999). For complete field acquisition parameters the reader may refer to Chapter 1; also see Cook et al. (1999). In the case of the wide-angle reflection experiment, 13 explosive shots ranging in size between 400 and 10,000 kg were used (Fig. 5.1). The total number of seismographs was 594 with geophone spacing of 1.2 km. Further details of the field acquisition can be found in Fernandez Viejo and Clowes (2003). Parameter Crooked line geometry Refraction statics correction Bandpass filter Deconvolution Automatic gain control NMO correction Stack Trace balance  Value 30 m x 2000 m bin From original processing 10-40 Hz Filter length=150 ms; gap length=2.0 ms 500 ms window V(z) derived from refraction study Nominally 134 fold Average amplitude scaled to be 1.0  Table 5.1. Basic processing parameters for NVI data prior to stochastic analyses.  The near-vertical data were processed using standard procedures (Table 5.1) similar to those found in Cook et al. (1999). Except for trace balancing, no further processing was applied to the stacked image. This is due to that fact that processes such as migration, coherency filtering and curvelet denoising increase the lateral coherency and smoothness of the data. In turn, this can lead to erroneous results when estimating the stochastic parameters. On the segment of the stacked profile that includes the Great Bear magmatic arc (Fig. 5.2), the arc is about 2.0–4.5 km thick and is delineated by a strong reflection from the top (~0.2 s), which has been drilled, and another strong reflection at ~1.5 s from its interpreted base (Cook et al., 1999). On the southwestern end of the profile, the lower crust contains a series of discrete, northeast-dipping reflections within a zone that is  126  around 3.0 s thick and extends to the Moho. The Moho is relatively flat over the width of the profile and well imaged at ~11.0 s except for the first 15 km where it is about 0.8 s deeper. The Moho transition zone is characterized by a sharp, narrow band of reflections (~300 ms), which separates the reflective crust from the transparent upper mantle.  Figure 5.2. Stacked, unmigrated seismic reflection profile beneath the Great Bear magmatic arc along SNORCLE Line 1. Moho reflections (dashed line) are visible at ~11.0 s between ~20 and 60 km and at ~11.5 s between 0 and 20 km. Block arrows bound a region of dipping reflectivity in the middle and lower crust. Inset shows location of the profile (white) with respect to Line 1.  Processing of the wide-angle data included editing/removing noisy traces and bandpass filtering between 3–15 Hz (Fig. 5.3). The crustal refracted phase Pg is the first arrival to offsets of ~160 km. Strong wide-angle reflections from the Moho transition zone along the same segment of Line 1 as the reflection data are observed on three different shots (1104, 1105 and 1108) between 80 and 220 km offset. This PmP phase is followed by a ~0.5 s-long coda. An intracrustal PcP phase reflects from the discontinuity between the upper-middle crust and the lower crust. At offsets greater than ~160 km, a  127  weak Pn phase traveling within the uppermost mantle becomes the first arrival. A 2-D crust and upper mantle seismic velocity model for all of Line 1 (Fig. 5.1) was derived from interpretation of the wide-angle dataset from all shots (Fernandez-Viejo and Clowes, 2003). In this study, we use the velocity model for the segment of Line 1 between shots 1104 and 1108 (Fig. 5.1) to provide the deterministic velocities used with the derived stochastic parameters.  Figure 5.3. Example of record section for shot 1105. Data are bandpass filtered from 3-15 Hz and displayed as trace-normalized amplitudes. Pg is a crustal, first arrival, refracted phase; PmP is a prominent, second arrival reflection from the Moho observed at offsets > 80 km, followed by ~0.5 long coda; Pn is a weak, refraction phase in the upper mantle observed at offsets > 160 km; PcP is an intracrustal phase reflected from an interface between the upper-middle crust and the lower crust.  5.4 Estimation of power law parameters Assuming that the continental lower crustal structure can be described by a 2-D von Karman autocorrelation function (ax, az, v) and this structure is inherited by the recorded seismic wavefield, we can back-estimate the lateral von Karman function (ax,v) from the reflection response. ax and az are the horizontal and vertical correlation lengths respectively that control the characteristic structural scale of the heterogeneous fabric; below this scale the heterogeneity manifests itself according to a power law with exponent v, which is effectively a measure of roughness. For the estimation we follow an  128  approach similar to that used by Holliger et al. (1994), Hurich and Kocurko (2000) and Carpentier and Roy-Chowdhury (2009). We compute normalized 2-D autocorrelation functions in windows of seismic data from which 1-D averaged lateral autocorrelation for the entire window is derived (Appendix B). The 1-D lateral autocorrelation curve is then fitted with a 1-D von Karman function as an estimator for the average horizontal von Karman content of the window, in terms of horizontal correlation length and power law exponent. The fitting is done in a least squares sense with ax and v supplied by a gridded model-parameter-space search. Accompanying the best fitting parameters estimates are associated maximum likelihood uncertainties, which are derived from probability density functions. To achieve maximum spatial resolution when mapping the variations in stochastic parameters, the 2-D autocorrelation windows must be as small as possible. However, to ensure stability and robustness for the lateral autocorrelation curve, the windows must be large enough for a proper sampling of the correlation lengths in space and the reflections in time. Thus, a compromise is needed regarding the window size. We have found from the literature and experience that around 2π correlation lengths should fit in the window spatially and 4 dominant periods of the source wavelet should fit in the time window.  5.5 Two-dimensional finite-difference modeling Another major element of our analysis is the calculation of synthetic seismograms for heterogeneous lithospheric models using a 2-D visco-elastic finite-difference wave propagation algorithm (Robertsson et al., 1994). The velocity models are defined on regular grids of cells for which the velocity values are derived from the earlier interpretation of the wide-angle data (Fernandez-Viejo and Clowes, 2003). The P-wave velocity increases from ~5.9 km/s at the surface to ~6.9 km/s at the crust-mantle boundary (~32 km depth). The velocity of the uppermost mantle is 8.1 km/s and linearly increases to 8.3 km/s at 43 km depth. Corresponding S-wave velocities are calculated using the ratios Vp/Vs = 1.78 and 1.82 for the crust and upper mantle respectively (Fernandez-Viejo et al., 2005). Crustal density values are ρ = 2700-2900 kg/m3 and the uppermost mantle density is ρ = 3300 kg/m3 (Christensen and Mooney, 1995; 129  Christensen, 1996). Intrinsic attenuation is also considered by using quality factors Qp = 500 and Qs = 200 for the crust and Qp = 1000, Qs = 400 for the upper mantle (Scheirer and Hobbs 1990; Fowler, 2005). The grid cell size and time step size of the calculations are chosen based on the stability criterion, which depends on the velocities and the maximum frequency of the signal propagating through the model (Robertsson et al., 1994). In order to maintain stability for the finite-difference algorithm up to a maximum frequency of 16 Hz, we use a grid spacing of 64 m and time steps of 2.34 ms in the case of the wide-angle simulations and values of 32 m and 1.17 ms in the case of the nearvertical data. This results in a grid with more than 2 million grid points for the 260 km long and 43 km deep wide-angle model (65 km long and 44 km deep for the near-vertical simulations). Three sources replicating shots 1104, 1105 and 1108 in addition to 160 receivers are placed along the surface of each model with 1.5 km receiver spacing. For the wide-angle synthetic seismograms, a zero-phase Ricker wavelet is used as the source signal with a center frequency of 8 Hz and a maximum frequency of 16 Hz, whereas a center frequency of 18 Hz and a maximum frequency of 32 Hz is used for the nearvertical seismograms. The velocity models are padded along the sides with a 10-km-wide damping zone characterized by very low Q value. This is done to ensure that the undesired edge effects are effectively suppressed (Nielson and Thybo, 2006).  5.6 Calibration approach The main purpose of the SNORCLE transect studies was to investigate the various processes involved in the westward growth of North America from the Archean to the present using a multidisciplinary approach, including the deep seismic reflection profile Line 1 (Cook and Erdmer, 2005). In this study we only examine that section of the line beneath the southern end of the Great Bear magmatic arc. Figure 5.2 shows the stacked unmigrated section, a 60 km long by 15 s (~48 km) deep profile. The processing parameters and procedures are given in Table (5.1) and a detailed interpretation was published by Cook et al. (1999). Our focus is to calibrate the lateral correlation length and the aforementioned upscaling factor of the deeper structures including the lower crust, the Moho transition zone and the uppermost mantle. This is achieved in the  130  following three steps. 1) estimation of the horizontal correlation length ax and its uncertainties directly from the near-vertical data; 2) construction of a suite of 2-D stochastic velocity models based on the measured and known fractal parameters; and 3) 2-D visco-elastic finite-difference forward modeling of near-vertical shot gathers and wide-angle record sections, for comparison with the observed data. 5.6.1 Estimation of von Karman parameters from SNORCLE Line 1 NVI data 5.6.1.1 Window size The windows used for the 2-D autocorrelations are 64 traces wide and 64 time samples deep (1920 m by 256 ms). This should ensure robust statistics for correlation lengths ax up to 306 m (1920/2π m) based on the criterion of Frenje and Juhlin (2000). Also, the dominant frequency of the near-vertical SNORCLE data is ~26 Hz, which corresponds to a period of 1000/26 Hz = 38.5 ms. Thus, the time-width of our window (256 ms) includes more than 6 cycles of the dominant period, which will ensure a robust average over the time direction. The sliding window shift is 1/4 the window length for both space and time leading to 28182 (122 by 231) overlapping windows in total. 5.6.1.2 Correlation lengths and their uncertainties Following the procedures outlined in section 5.4, the estimated correlation lengths were generated from the NVI data section (Fig. 5.4a). The vast majority of the estimated values of ax lie beneath 306 m, except for some high values (> 500 m) mainly in the shallowest part of the section at ~0.5 s. The map shows a clear differentiation between the crust and the uppermost mantle where the latter is characterized by much smaller ax values (about a factor of 3). The Moho is seen as a distinct feature separating the crust from the upper mantle and is continuously horizontal across the entire profile, except at ~18 km where it is offset by about 0.8 s. The outline for the base of the Great Bear magmatic arc (at < 0.5) is well delineated by very high correlation length with respect to the surrounding crust. Narrow undulating stripes of low ax values are visible in the profile at ~4.0 s between 40 and 60 km and just above the Moho transition between 25 km and 45 km. Other patches of low ax values are visible in the middle crust at ~20 km.  131  Figure 5.4. (a) Map of the estimated correlation lengths (ax) from the reflection data (Fig. 5.2). White rectangle indicate the area over which the average ax value is measured for the Moho transition zone. (b) Associated maximum likelihood uncertainties in meters.  132  Figure 5.4b displays the maximum-likelihood uncertainties in ax in the form of standard deviation in meters. The overall low to moderate values of the uncertainties, especially at the Moho transition and the uppermost mantle, give us some confidence in the estimates for ax and for further interpretation. There is no obvious relation between the ax uncertainties and regions with high reflective power (high SNR) or with low reflective power (low SNR). 5.6.1.3 Power law exponent estimates Estimates of the power law exponent v (Hurst number) show a similar pattern to that of ax but they are associated with relatively larger uncertainties (Fig. 5.5a and b). However, there is a general trade-off between estimated ax and v. Zones of high ax tend to correlate with zones of low v and vise versa. Such a trade-off has been observed in earlier studies by Hurich and Kocurko (2000), Carpentier and Roy-Chowdhury (2007), Carpentier et al. (2009) and Carpentier et al. (2010). Some estimates exceed the nominal maximum value of 1.0; however, the majority (except within the uppermost mantle) fall into the strict von Karman parameter regime (0-1). An increase in the values of v with depth, indicating increasing smoothness in the data with depth is partially due to the intrinsic attenuation of the elastic wavefield. From the preceding discussion and other studies (e.g., Pullammanappallil et al, 1997; Carpentier et al., 2009), the Hurst number is not a very reliable parameter to use for constraining lithospheric fabric. 5.6.1.4 Goodness of fit The average least squares residuals (Fig. 5.6) is another quantity for assessing the quality of the von Karman content of a seismic window (Carpentier et al., 2009). These residuals are a direct measure of how well the measured correlation curves are fit by the best fitting von Karman curve in a window. The least squares value map displays a relatively uniform distribution of misfit over the seismic profile, which is a good sign in terms of lacking bias. However, the residuals tend to be slightly higher on the left side of the map than on its right side. This is probably due to more diffractions crossing the data on the left. Another interesting trend is the decrease of residuals at the uppermost mantle level. Normally, one expects a general decrease in the quality of seismic data due to attenuation 133  Figure 5.5. (a) Map of estimated power law v exponent from the reflection data (Fig. 5.2). (b) Associated maximum likelihood uncertainties in the power law exponents.  134  Figure 5.6. Map of average minimum least squares residuals for best fit. Note that the misfit distribution over the profile is relatively flat. This means that there is little bias in the data regarding the estimates.  and lower signal-to-noise ratios. In this case however, it seems that the high power of the seismic sources, the acquisition setup and excellent recording conditions have prevented such usual data deterioration at these late arrival times. This is consistent with the recording of reflections from within the upper mantle as discussed in Chapter 2. 5.6.2 Synthesis of velocity models The stochastic velocity fields used for the finite-difference modeling have a structure according to the radial von Karman power spectrum P(k), which is given by P(k) = 4πνaxaz/(1+k2)ν+1 where  (5.1)  is the weighted radial wavenumber, kx = 2π/λx and kz = 2π/λz  are the horizontal and vertical wave numbers respectively, ax is the horizontal correlation length, az is the vertical correlation length and ν is the Hurst number. For creating the discrete 2-D stochastic field, we follow the procedures described by Goff et al. (1994). First we define the 2-D von Karman power spectrum as a discrete signal. Then we apply a uniform random phase φ to the square root of the signal and finally an inverse Fourier transform yields a 2-D stochastic field. The velocities are then  135  mapped to a tri-modal distribution while preserving the von Karman power spectrum of the structure (Goff and Holliger, 1994; Carpentier and Roy-Chowdhury, 2007). Here we test two substantially different models that describe the crust-mantle boundary. In Model 1, the Moho represents a 3-km-thick, stochastic transition zone with velocities characteristic of both the lower crust and upper mantle (Fig. 5.7a). The reason for choosing such a model is based on our results from Chapter 4, which showed that a 3km-thick zone was needed to produce the observed 0.5-s-long coda. This stochastic transition zone is then embedded within a deterministic velocity model described earlier by Fernandez-Viejo and Clowes (2003). In Model 2, the Moho is an interface/change in the correlation length between a stochastic uppermost mantle (35-43 km depth interval) and a stochastic lower crust (22-35 km depth interval) overlain by a deterministic upper crust; i.e, there is no Moho transition zone (Fig. 5.7b). From the correlation length map and its associated uncertainties we calculate the average ax value over the Moho transition (Fig. 5.4a), which gives a value of 86±17 m. However, because of the systematic underestimation of ax (Hurich, 1996; Hurich and Kocurko, 2000; Poppeliers and Levander, 2004; Carpentier and Roy-Chowdhury, 2007), we upscale the horizontal correlation length and its uncertainty for one model by a factor of 6 and for another model by a factor of 9 (Hurich, 2003; Carpentier and RoyChowdhury, 2007 and Carpentier and Roy-Chowdhury, 2009). To calculate the values of az, we use the aspect ratios (ax/az) of 4 and 10 as proposed by others (Holliger and Levander, 1992; Holliger, 1996). For the Hurst number ν we use a value of 0.3 as extracted from outcrops in the Ivrea Zone by Holliger and Levander (1992); as noted previously, Hurst numbers derived from the data are not reliable. These fractal parameters lead to 8 different combinations of a tri-modal velocity distribution (Table 5.2). The associated P-wave velocities of the tri-modal distribution are 6900, 7500 and 8100 m/s (3864, 4125 and 4455 m/s S-velocities). Figure 5.7a shows two examples of the discrete von Karman velocity fields representative of the Moho transition zone. The upper model in (a) has the smallest horizontal and vertical correlation lengths (ax= 409 m and ax/az =10), hence we call it the absolute minimum model. Whereas the lower model in (a) has the largest horizontal and vertical correlation lengths (ax=936 m and ax/az =4), 136  Figure 5.7. Tri-modal velocity fields with two distinct Moho models embedded within a deterministic background. (a) Model 1 includes a 3-km-thick Moho transition, (top) absolute min. correlation length; (middle) absolute max. correlation length, see Table 5.2. (b) Model 2 contains a Moho represented as an interface between stochastic lower crust: ax=732 m, az=183 m, and stochastic upper mantle: ax=261 m, az=62 m. Stars indicate locations of the three shot points for which wide-angle data are computed.  137  hence we call it the absolute maximum model. Note in these models the lack of periodicity in the heterogeneities and their self-similar character. For Model 2 we follow the same procedures described above to assign ax, az and ν values for the lower crust (between 7.0 and 11.0 s) and uppermost mantle (below 11.0 s), which are both characterized by a tri-modal velocity distribution (Table 5.2). The associated P-wave velocities for the lower crust are 6700, 6800 and 6900 m/s and for the upper mantle 8000, 8100 and 8200 m/s (Fig. 5.7).  Velocity models Model 1  az (m)  ax (m)  az (m)  (86-18) x 6 = 409  102  409  40.9  (86+18) x 6 = 614  154  614  61.4  (86-18) x 9 = 624  156  624  62.4  (86+18) x 9 = 936  234  936  93.6  462  116  462  46.2  abs. max  732  183  732  73.2  abs. min abs. max  135 261  34 62  135 261  13.5 26.1  Abs. min. correlation length Min. correlation length Max. correlation length  Upper mantle  ax/az = 10  ax (m)  Abs. max. correlation length Lower crust abs. min Model 2  ax/az = 4  Table 5.2. Types of velocity models derived from stochastic parameters.  5.6.3 Simulation of wide-angle and near-vertical data The calculated synthetic wide-angle seismograms for the different combinations of the transition zone (Model 1) display a prominent PmP phase between ~80 km and 210 km offset trailed by ~0.5 s coda, consistent with the observed data, and a weak but visible Pn phase at offsets greater than 170 km (Fig. 5.8 c-f). The changes in the appearance of the PmP coda and its duration as a result of the different correlation lengths are subtle. However, models with larger correlation lengths tend to produce longer and more pronounced reverberatory coda. Synthetic seismograms resulting from Model 2 show a strong PmP phase but very weak Pn phase (Fig. 5.8 g-h). In addition, the intracrustal PcP phase, which can be attributed to the interface between the deterministic upper crust and the stochastic lower crust, is clearly visible. The major difference between Model 1 and 2 is the absence of a complex coda following the PmP phase in Model 2, regardless of the size of the correlation lengths and their aspect ratios. 138  Figure 5.8. Observed (a-b) and synthetic (c-h) wide-angle seismograms for shots 1105 and 1108. Different seismic phases are explained in the text. (c-d) Synthetic seismograms from Model 1 with absolute max. correlation length. (e-f) Synthetic seismograms from Model 1 with absolute min. correlation length. (g-h) Synthetic seismograms from Model 2.  139  Figure 5.9. Observed (a) and synthetic (b-d) near-vertical seismograms. (b) Synthetic seismogram from Model 1 with absolute min. correlation length and ax/az=10. (c) Synthetic seismogram from Model 1 with absolute max. correlation length and ax/az=4. (d) Synthetic seismogram from Model 2 with ax/az=10. Enlargements of Moho reflections are shown in the insets.  140  The synthetic near-vertical seismograms from Model 1 image the Moho transition as a multicyclic band of reflectors (Fig. 5.9b and c) that are laterally variable and piecewise continuous, in a good agreement with the observed data (Fig. 5.9a). On the other hand, the Moho transition in Model 2 consists of a single strong reflector with no lateral variability (Fig. 5.9 d). However, this model successfully explains the reflective response of the lower crust and the lack of reflectivity in the upper mantle.  5.7 Attributes analysis Unlike Model 1, Model 2 fails to simulate the observed Moho signature in both datasets; hence, further investigation with this model is not pursued. Although Model 1 does not currently produce the PcP phase, it produces the PmP phase, needed for investigation. Therefore, we will use Model 1 for the specific goal of faithfully modeling the Moho reflection response. However, visual inspection is insufficient to favor specific correlation lengths within Model 1 due to the subtle differences among the wide-angle seismograms. Thus, we further quantify the PmP phase by means of careful examination of complex trace envelope, envelope energy decay and frequency-space panels. The nearvertical data were not used for attribute analyses for two reasons: (1) limited offsets (< 11 km) are not suitable for amplitude vs. offset (AVO) and frequency-space analysis and (2) multicyclic reflections at Moho level do not exhibit any visible energy decay. 5.7.1 Complex trace envelope Before comparing the trace envelopes between the synthetic and real data, a number of corrections were first needed. The various types of geophones in the field were initially calibrated and appropriate scaling factors were applied to the wide-angle trace amplitudes to achieve ‘true relative amplitude’ (Ellis et al., 2002). To correct for geometrical spreading in the case of the synthetic seismograms, all trace amplitudes were multiplied by a factor of r½ (r = source-receiver distance), assuming that a point source in 2-D is equivalent to a line source in 3-D. We calculate the trace envelopes for the PmP phase, also known as instantaneous amplitude, as a function of offset (Fig. 5.10). To compute the envelope curves, we extract  141  Figure 5.10. Computed PmP envelope curves vs. offset for observed and synthetic seismograms of the three different shots. Observed data (solid line), absolute maximum correlation length model with ax/az =4 (dotted line) and absolute minimum correlation length model with ax/az=10 (short dashed line). Trend lines are calculated based on the least-squares fit and are shown in black.  142  the PmP amplitudes on each trace over a 0.5 s long window. Then we add the Hilbert transform of that window as the imaginary part and take the absolute value of this complex signal. The curves measured from the real data display a complex and erratic behavior with offset and lack any distinct peak. For shot 1104, the envelopes were not calculated beyond 175 km offset distance due to the overlap between the Pg and the PmP phases. The envelope curves measured from the various synthetic seismograms, also exhibit complex and erratic behavior but with high peaks at ~120 km offset (Fig. 5.10), which can be explained best as the critical distance. All the curves share the same general trend of decreasing energy with offset as in the real data but any difference caused by the extremes in correlation length is too subtle to notice. 5.7.2 Envelope energy decay In order to reveal the subtle variations in trace envelopes for the different synthetic seismograms, we measure the envelope energy decay per trace. This is accomplished by fitting an exponential function with a certain exponent to the decay curve for a trace using the least-squares method (Fig. 5.11). The decay exponents are then plotted as a function of offsets of the traces for the different seismograms (Fig. 5.12). The decay curves calculated from the observed data (Fig. 5.12 solid line) show smaller exponent values on all three shots suggesting a slower decay rate. For the synthetic data, there is a small but pervasive trend of higher decay exponents for the absolute minimum von Karman media (Fig. 5.12 dashed line). This means that the coda dissipates faster with smaller correlation lengths. Also, there is a small increase in the exponent values at farther offsets, which can be attributed to attenuation. As a result, the absolute maximum correlation lengths can better explain the observed data on this point. 5.7.3 F-X analysis We calculate the trace-balanced, continuous frequency spectra by taking the Fourier transform of each trace and plotting all traces side by side. The frequency content of the observed PmP phase fluctuates with offset between 5 and 10 Hz and the average frequency is around 8 Hz (Figs. 5.13 and 5.14). On the other hand, average frequency of the synthetic PmP phase is slightly lower than that of the observed data by about 1.5 Hz 143  Figure 5.11. Example of best-fitting exponential curve (dashed) to the envelope decay (solid) of one trace at a specific offset.  Figure 5.12. Envelope decay curves for observed and synthetic seismograms of the three different shots. (ac) Comparison of decay curves between observed (solid line) maximum (dotted line) and minimum (dashed line) correlation lengths for models with ax/az=4. (d-f) Same as before but for models with ax/az=10.  144  Figure 5.13. Trace-balanced F-X spectra maps for observed (a-c) and synthetic (d-i) seismograms. (d-f) Maps for the absolute max. correlation length model with ax/az=4. (g-i) Maps for the absolute min. correlation length model with ax/az=10. The trace in the right panel of the maps represents the average F-X spectra of all traces.  145  Figure 5.14. Trace-balanced peak frequency curves for observed and synthetic seismograms of three different shots. Observed data (solid line), absolute max. correlation lengths model with ax/az=4 (dotted line) and absolute min. correlation lengths model with ax/az=10 (dashed line).  (~6.5 Hz). However, between the same shots there is a persistent trend of narrower range of frequencies for the absolute maximum correlation length (right panel on F-X maps in Fig. 5.13). In other words, there is a more pronounced peak frequency on each shot, resembling the sharper peak frequencies of the observed data. The peak frequency curves of the absolute maximum correlation lengths in Figure 5.14 are slightly closer to the ones of the observed data as well.  5.8 Discussion 5.8.1 Fabric in the crust and uppermost mantle The heterogeneity information derived from mapping directly the von Karman parameters of the seismic reflection wavefield is closely related to the spatial distribution of acoustic impedance fluctuations or macroscale petrofabrics. The petrofabrics are in turn the result of a combination of microscale processes such as fracturing and partial melting. Thus, variations in stochastic parameters may indicate a transition from one lithospheric fabric to another or from one mode of geodynamical process to another. On the other hand, clustering of parameters with similar values may imply domains of similar rheological behavior. A significant difference between the lateral correlation functions and reflection amplitudes lies in the fact that when the former are computed from the seismic data, they are normalized by amplitude and in turn insensitive to reflection strength. As a result, regions in the seismic section with low amplitudes can still offer important internal delineations based on their different correlation lengths (Carpentier et al., 2009). 146  Comparison of the correlation length map with an earlier geologic interpretation (Fig. 5.15) based on the same reflection data (Cook et al., 1999) shows good agreement between the two sections. One important observation is that zones of higher correlation lengths do not generally coincide with higher amplitudes in the reflection data and vice versa. The Great Bear arc rocks probably do not possess a von Karman structure due to their bulk homogeneity; hence they display artificially high correlation lengths (Fig. 5.15a). Beneath the Great Bear magmatic arc, Cook et al., (1999) divided the Hottah terrane into upper (HC1), middle (HC2) and lower (HC3) crust. To a certain extent, the boundaries between these three zones can be associated with the undulating low correlation length zones (short dashed lines). Associated low uncertainties and small residuals (Figs. 5.4 and 5.6) support the validity of these low correlation length values. HC1 and HC2 zones were interpreted as being uplifted and folded in response to the accretion of the Hottah terrane onto the Slave craton. Within HC3, the main feature is a series of discrete, east-dipping layered structures that are also characterized by low correlation lengths (dotted lines). Thus, there appears to be a consistent relation between these detachment faulting/contractional folding surfaces and decreasing correlation lengths. In contrast, the Moho transition zone is characterized by relatively higher correlation lengths (long-dashed lines). In previous studies (Carpentier et al., 2009; Carpentier et al., 2010), horizontal mafic intrusions (sills) have been suggested to explain the increasing correlation lengths at Moho levels. This is consistent with our findings in Chapter 4 (Oueity and Clowes, 2010). Mafic sills intruded from the upper mantle into the crust can cause elongated reflections with large lateral correlation lengths (Lyngsie et al., 2007). These intrusions are probably part of regionally extensive, postaccretion processes that have effectively modified rocks and re-equilibrated the reflection Moho to a uniform depth (~33 km) (Oueity and Clowes, 2010; Thybo and Nielsen, 2010) . In Chapter 4, we proposed that the source of thermal activity responsible for magmatic intrusions in the study area is most likely associated with the deposition of the Great Bear magmatic rocks between ~1.88 to 1.84 Ga.  147  Figure 5.15. (a) Interpretation for the correlation length map. Short-dashed lines delineate boundaries between the crustal regions; dotted-lines follow east-dipping surfaces; long-dashed lines delineate the Moho transition zone; white line indicate a possible fault. (b) Earlier interpretation based on reflection data (from Cook et al., 1999). HC1, HC2 and HC3 are upper, middle and lower Hottah crust respectively. F is a thrust fault.  148  A west-dipping surface (F; white line), previously interpreted as a thrust fault that projects to about 6.0 s, is probably responsible for the juxtaposition between the short dashed lines of low correlation lengths to the southwest and higher values to the northeast. This could indicate slightly different fabrics for the rocks on either side of the interpreted fault, but the results are too subtle to be definitive. The uppermost mantle is characterized by low correlation lengths. Associated good fits and low uncertainties give us confidence that the correlation length values are valid. The high quality seismic data is partly responsible for this. A reason for the lower correlation lengths is probably that the origin of the stochastic fabric in the upper mantle is fundamentally different from that of the lower crust. When considering that the upper mantle is the departing point of the intrusions into the lower crust, one can imagine that, to a first order, the migration paths of the partial melts are vertical. Vertically migrating magmas will not only have a lesser expression in horizontal width, but they can also interrupt existing horizontal lineations in the upper mantle. Cook and Vasudevan (2003) have applied a combination of image enhancement techniques on SNORCLE Line 1 reflection data in order to analyze the uppermost mantle structure beneath the Great Bear magmatic arc. The improved image shows steeply west dipping reflections just below the Moho. One explanation they proposed for such reflections is that they could be magmatic intrusions related to the Great Bear arc. This is in agreement with the interpretation of our results. Other evidence for low uppermost mantle correlation lengths comes from the work of Thybo (2006) and Nielsen and Thybo (2006). Based on their interpretation for the lack of pronounced Pn coda in many wide-angle reflection datasets, they argue against the possibility of well developed horizontally elongated heterogeneity in the uppermost mantle between 40-100 km depth. 5.8.2 Constraining the correlation length and the upscaling factor The diffusive reflection response of the Great Bear lower crust, Moho and uppermost mantle regions can be represented through stochastic models of heterogeneity, for which the stochastic parameters were chosen from our analysis and other studies (Holliger and Levander, 1992; Carpentier and Roy-Chowdhury, 2009). This is one of the goals of our  149  study, but the second goal relates to one of the stochastic parameters under consideration. Von Karman correlation length for the Moho region is still an uncertain parameter, even though it can be directly estimated from the associated near-vertical reflection data. The inverted correlation lengths have a maximum-likelihood uncertainty and, in addition, the upscaling factor to correct for the underestimation ranges between 6 and 9. Our second goal in this study is to constrain the correlation length and its upscaling factor with the wide-angle reflection data. We achieve this by finite-difference forward modeling and comparison of attributes between synthetic and observed data. In the comparison of trace energy envelope, envelope decay and F-X spectra between the observed and synthetic data, we noticed that the values of these attributes are comparable. Figures 5.12-14 demonstrate that following the same processing steps for the observed and synthetic data, values for envelope decay and peak frequency match relatively well. This indicates that our choice for the current model of heterogeneity and its variations is reasonably appropriate. Unfortunately, the sensitivity of the compared attributes to the correlation lengths under investigation is not high. Apparently, the energy of the backscattered wavefield does not really depend on correlation length in heterogeneity, given the same recording geometry (Fig. 5.10). However, the distribution of backscattered energy with time does depend on correlation length as the envelope decay in time suggests (Fig. 5.12). The envelope decay rates for larger correlation lengths are lower, enabling a richer coda as seen in the observed data as well. F-X spectra point to the same conclusion; the wavefield reflected from a model with absolute maximum correlation length exhibits peaks in frequency that are relatively sharper (right panels on Fig. 5.13d-f). This is in better agreement with the more pronounced peak frequencies seen on the observed data (right panels on Fig. 5.13a-c). Although the evidence is not compelling and the effects are subtle, the comparisons between synthetic and observed data appear to favor the models with absolute maximum correlation lengths and an upscaling factor of 9 for the Moho heterogeneity.  150  5.8.3 Character and attributes of wide angle reflection data From the previous section we determined that the Moho model with absolute maximum correlation lengths best explained the observed wide angle reflection data. In this section we provide physical explanations for the observed changes in attributes as a function of correlation length. One important realistic property of the reflection response of the model with absolute maximum correlation length was its lower decay in coda envelope. We propose that this lesser decay with larger correlation length can be explained by scattering attenuation. As studies of Frenje and Juhlin (2000) and Hong and Kennett (2003) show, forward and multiple scattering in complex heterogeneous media will cause destructive interference in the transmitted and reflected wavefield. This type of destructive interference is often called scattering attenuation and is a complex function of many medium and wavefield parameters. Incidence and reflection angle, elastic parameter ratio, scale length of heterogeneity and wavelength are all highly influential parameters. Both studies mentioned above have done a careful analytical and numerical analysis of the strength of scattering attenuation as a function of these parameters. In particular, 1/Qs (the amount of scattering attenuation) as a function of the product of k (wavenumber of incident wavefield) and a (correlation length) is a well established relation. If we inspect the graphs of 1/Qs versus ka in both studies, we observe that depending on incidence angle and P/S velocity ratio, there is always a maximum of 1/Qs at a certain ka value. We can determine our values of ka for the absolute minimum and absolute maximum correlation length models: taking k = 2π/λ, λ = Vave / fcenter, Vave = 6900 m/s, fcenter = 8 Hz and considering the two extreme values for a, 409 m and 936 m, we get two values for ka: 3.02 and 6.92. Taking into account that we are dealing with large angles of incidence (θ > 30°) and a P/S ratio (γ) of about 1.78, we make an interesting observation. In both studies, the scattering attenuation for our parameters ka, θ and γ is seen to decrease with larger correlation length a. This provides a good explanation of the less attenuated coda that we observe for the Moho model with absolute maximum correlation length.  151  Another realistic attribute characteristic of the absolute maximum correlation length model was the peak frequency versus offset function, which features a sharper peak frequency (Fig. 5.13) and was overall closer to the observed data than the other model (Fig. 5.14). An explanation for the sharper peak frequency cannot be found in frequency tuning by dominant thickness in layering, as there are no typical scale lengths or layering in von Karman media. This fact has been demonstrated by Holliger et al. (1994), who presented a relevant theoretical background and synthetic example. A valid explanation could come again from scattering attenuation. In the two aforementioned studies of Frenje and Juhlin (2000) and Hong and Kennett (2003), no frequency dependency of scattering attenuation was examined and discussed. It could however be possible that a frequency dependency in scattering attenuation is quite prominent and that we are witnessing the effect of this in our modeled data. A further examination of the scattering attenuation theory is necessary. We can only say that the absolute maximum correlation length model has more realistic peak frequencies with respect to the observed data. 5.8.4 The Moho transition zone We have calculated the two-dimensional visco-elastic seismic response of two end members for the crust-mantle boundary (Fig. 5.7). Model 1 includes a 3-km thick heterogeneous transition zone sandwiched between homogeneous crust and upper mantle, both with positive velocity gradients. In Model 2, the Moho is an interface between heterogeneous lower crust and upper mantle. The velocity fluctuations within the heterogeneous zones follow a von Karman distribution function characterized by different correlation lengths. Both models cannot fully explain all the phases and scattered features (i.e., coda) observed on the wide-angle and near-vertical data. Thus, our preferred representation of the Moho transition zone is a combination of Models 1 and 2 (Fig. 5.16a). This model includes all three stochastic regions, each with a Hurst number of 0.3: (1) the lower crust (22-33 km depth interval) with horizontal and vertical correlation lengths of 732 and 183 m respectively; (2) the Moho transition (33-36 km depth interval) with horizontal and vertical correlation lengths of 936 and 234 m respectively and (3) the uppermost mantle (36-43 km depth interval) with horizontal and vertical correlation lengths of 261 and 62 m respectively. The synthetic seismograms 152  Figure 5.16. (a) Composite velocity model including heterogeneity in the lower crust (22-33 km), the Moho transition (33-36 km) and the uppermost mantle (below 36 km). (b) Observed wide-angle record section. (c) Observed near-vertical shot gather. (e) Synthetic wide-angle seismogram. Different seismic phases are explained in the text. (d) Synthetic near-vertical seismogram. The distinct reflector at ~5.5 s is due to the velocity interface at around 15 km depth.  153  (Fig. 16 d and e) provide a good match to all the main features of the observed data (Fig. 16 b and c). The near-vertical seismogram exhibits reflective lower crust, a narrow band of multicyclic reflections at the Moho level and relatively transparent uppermost mantle. On the wide-angle section, we notice a distinct coda following the PmP phase, a clear intracrustal PcP phase, and a very weak Pn phase.  5.9 Conclusions Estimation of the von Karman stochastic parameters (lateral correlation length ax and power law exponent v) for the Wopmay orogen beneath the Great Bear magmatic arc were carried out using LITHOPROBE’s SNORCLE Line 1 seismic reflection data. Simulated wide-angle record sections and near-vertical shot gathers from representative stochastic velocity fields were compared to observed data and used to constrain the stochastic nature of the Moho transition zone and its associated lateral correlation length. Our results show that unlike the lateral correlation length ax, the power law exponent v is not a reliable parameter to use for constraining the lithospheric fabric. Thus, for the purpose of this study we use a value of v = 0.3 taken from the literature. The Great Bear arc rocks are characterized by artificially high correlation lengths implying a relatively homogeneous fabric. Narrow bands of low correlation lengths can be linked to intracrustal boundaries and east-dipping layered structures previously interpreted as faulting/folding surfaces. The Moho transition zone is associated with higher correlation lengths possibly due to upper mantle mafic intrusions causing elongated reflections. On the other hand, the uppermost mantle is characterized by low correlation lengths, which we interpret as the result of vertically migrating magmas that may interrupt existing horizontal lineations in the upper mantle. Comparisons between the simulated synthetic seismograms and observed data indicate that the lower crust, the Moho transition zone and the uppermost mantle all possess a von Karman-type heterogeneity distribution, each characterized by different lateral correlation lengths. The Moho transition zone is approximately 3-km-thick between 33 and 36 km depth interval. This transition zone is responsible for the generation of the pronounced coda trailing the PmP phase seen on the wide-angle data and the multicyclic 154  band of reflections seen at Moho level on the near-vertical data. The results from the seismic attribute analysis favor an upscaling factor of 9 resulting in lateral correlation length values of 732 m, 936 m and 261 m for the lower crust, the Moho transition zone and the uppermost mantle, respectively.  155  Bibliography Carpentier, S. F. A., 2007. On the estimation of stochastic parameters from deep seismic reflection data and its use in delineating lower crustal structure [Ph. D. Thesis]: Utrecht, University of Utrecht, 147 p. Carpentier, S. and Roy-Chowdhury, K., 2007. Underestimation of scale lengths in stochastic fields and their seismic response: a quantification exercise. Geophysical Journal International, 169, p. 547–562. Carpentier, S. and Roy-Chowdhury, K., 2009. Conservation of lateral stochastic structure of a medium in its simulated seismic response. Journal of Geophysical Research, 114, doi:10.1029/2008JB006123. Carpentier, S. F. A., and Roy-Chowdhury, K., Stephenson, R. A., and Stovba, S., 2009. Delineating tectonic units beneath the Donbas Fold Belt using scale lengths estimated from DOBRE 2000/2001 deep reflection data. Journal of Geophysical Research, 114, B10315, doi: 10.1029/2008JB006124. Carpentier, S. F. A., Roy-Chowdhury, K., and Hurich, C.A., 2010. Mapping correlation lengths of lower crustal heterogeneities together with their maximum-likelihood uncertainties. Tectonophysics, doi:10.1016/j.tecto.2010.07.008. Christensen, N. I., and Mooney, W. D., 1995. Seismic velocity structure and composition of the continental crust. A global view: Journal of Geophysical Research, 100, p. 9761–9788. Christensen, N. I., 1996. Poisson's ratio and crustal seismology. Journal of Geophysical Research, 101, p. 3139–3156. Cook, F.A., van der Velden, A.J., and Hall, K.W., 1999. Frozen subduction in Canada's Northwest Territories. LITHOPROBE deep lithospheric reflection profiling of the western Canadian Shield: Tectonics, 18, p. 1-24. Cook, F. A., and Vasudevan, K., 2003. Are there relict crustal fragments beneath the Moho?. Tectonics, 22, p. 1026, doi:10.1029/2001TC001341. Cook, F. A., and Erdmer, P., 2005. An 1800 km cross section of the lithosphere through the northwestern North American plate: lessons from 4.0 billion years of Earth’s history. Canadian Journal of Earth Sciences, 42, p. 1295–1311, doi: 10.1139/E04-106.  156  Dolan, S., Bean, C., and Riollet, B., 1998. The broadband fractal nature of heterogeneity in the upper crust from petrophysical logs. Geophysical Journal International, 132, p.489–507. Ellis, R., Clowes, R.M., Spence, G.D., Asudeh, I., Hajnal, Z., and Amor, J.R., 2002. SlaveNorthern Cordillera Refraction Experiment (SNORE’97): Field Acquisition and Preliminary Data Processing Report. University of British Columbia, LITHOPROBE Report No. 83, 14 pp, 10 apps. Fernandez Viejo, G., and Clowes, R.M., 2003. Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay orogen, northwestern Canada, from a LITHOPROBE refraction/wide-angle reflection survey. Geophysical Journal International, 153, p. 1-19. Fernandez Viejo, G., Clowes, R.M., Welford, J. K., 2005. Constraints on the composition of the crust and uppermost mantle in northwestern Canada: Vp/Vs variations along Lithoprobe's SNorCLE transect. Canadian Journal of Earth Sciences, 42, p. 1205-1222, doi:10.1139/e05028. Fountain, D. M., 1976. The Ivrea-Verbano and Strona-Ceneri zones, northern Italy: A cross section of the continental crust-New evidence from seismic velocities of rock samples. Tectonophysics, 33, p. 145-165. Fountain, D. M., and Salisbury, M. H., 1981. Exposed cross-section through the continental crust. Implications for crustal structure, perology and evolution, Earth and Planetary Science Letters, 56, p. 263-277. Fowler, C. M. R., 2005. The Solid Earth: An Introduction to Global Geophysics. Cambridge, UK, Cambridge University Press, 685 p. Frenje, L. and Juhlin, C., 2000. Scattering attenuation: 2-d and 3-d finite difference simulations vs. theory. Journal of Applied Geophysics, 44, p. 33–46. Gandhi, S.S., Mortensen, J.K., Prasad, N., and van Breemen, O., 2001. Magmatic evolution of the southern Great Bear continental arc, northwestern Canadian Shield: geochronological constraints. Canadian Journal of Earth Sciences, 38, p. 767-785. Gibson, B., 1991. Analysis of lateral coherency in wide-angle seismic images of heterogeneous targets. Journal of Geophysical Research, 96, p. 10,261–10,273. Goff, J. and Jordan, T., 1988. Stochastic modeling of seafloor morphology: inversion of sea beam data for second-order statistics. Journal of Geophysical Research, 93, p.13589–13608. 157  Goff, J., Holliger, K., and Levander, A., 1994. Modal fields: A new method for characterization of random seismic velocity heterogeneity. Geophysical Research Letters, 21, p. 493–496. Goff, J. and Levander, A., 1996. Incorporating ”sinuous connectivity” into stochastic models of crustal heterogeneity: Examples from the Lewisian gneiss complex, Scotland, the Fransiscan formation, California, and the Hafafit gneiss complex, Egypt. Journal of Geophysical Research, 101, p. 489–501. Goff, J. and Jennings, J. J., 1999. Improvement of fourier-based unconditional and conditional simulations for band limited fractal (von karman) statistical models. Mathematical Geology, 31, p. 627–649. Green, D. H., and Ringwood, A. E., 1972. A comparison of recent experimental data on the gabbro-garnet granulite-eclogite phase transition. Journal of Geology, 80, p. 277-288. Haussdorf, F., 1918. Dimension und äußeres Maß, Mathematische Annalen, 79, p. 157-179. doi: 10.1007/BF01457179. Hildebrand, R. S., Hoffman, P.F., and Bowring, S.A., 1987. Tectono-magmatic evolution of the 1.9-Ga Great Bear magmatic zone, Wopmay orogen, northwestern Canada. Journal of Volcanology and Geothermal Research, 32, p. 99-118. Hildebrand, R. S., and Bowring, S. A., 1999. Crustal recycling by slab failure. Geology, 27, p. 11-14. Hoffman, P.F., 1989. Precambrian geology and tectonic history of North America. In Bally, A. W., and  Palmer, A. R., eds., The Geology of North America: An Overview, v. A, p. 447-  512, Geological Society of America, Boulder , Colorado. Hoffman, P. F., and Bowring, S. A., 1984. Short-lived 1.9 Ga continental margin and its Destruction. Wopmay orogen, NW Canada: Geology, 12, p. 68-72. Holliger, K., 1996. Upper-crustal seismic velocity heterogeneity as derived from a variety of pwave sonic logs. Geophysical Journal International, 125, p. 813–829. Holliger, K. and Levander, A.,1992. A stochastic view of lower crustal fabric based on evidence from the Ivrea Zone. Geophysical Research Letters, 19, p. 1153–1156. Holliger, K. and Levander, A., 1994a. Seismic structure of gneissic/granitic upper crust: geological and petrophysical evidence from the Strona-Ceneri Zone (northern italy) and  158  implications for crustal seismic exploration. Geophysical Journal International, 119, p. 497– 510. Holliger, K., Carbonell, R., and Levander, A., 1992. Sensitivity of the lateral correlation function in deep seimic reflection data. Geophysical Research Letters, 19, p. 2263–2266. Holliger, K., Levander, A., and Goff, J., 1993. Stochastic modeling of the reflective lower crust: petrophysical and geological evidence from the Ivrea Zone (northern Italy). Journal of Geophysical Research, 98, p. 11967–11980. Holliger, K., Levander, A., Carbonell, R., and Hobbs, R., 1994. Some attributes of wavefields scattered from ivrea-type lower crust. Tectonophysics, 232, p. 267–279. Hong, T. K. and Kennett, B. L. N., 2003. Scattering attenuation of 2D elastic waves: theory and numerical modeling using a wavelet-based method. Bulleting of the Seismological Society of America, 93, p. 922-938. Hurich, C., 1996. Statistical description of seismic reflection wavefields: a step towards quantitative interpretation of deep seismic reflection profiles. Geophysical Journal International, 125, p. 719–728. Hurich, C., 2003. The nature of crustal seismic heterogeneity: A case study from the Grenville Province. In Goff, J. A. and Holliger, K., eds., Heterogeneity in the crust and upper mantle: nature, scaling, and seismic properties, Kluwer Academic Publishers, New York, p. 299–320. Hurich, C. and Kocurko, A., 2000. Statistical approaches to interpretation of seismic reflection data. Tectonophysics, 329, p. 251–267. Ito, K., and Kennedy, G.C., 1971. An experimental study of the basalt-garnet granulite-eclogite transition. In Heacock, J.G., ed., The Earth beneath the continents: American Geophysical Union Geophysical Monograph 14, p. 303–314. Levander, A., Hobbs, R., Smith, S., England, R., Snyder, D., and Holliger, K., 1994. The crust as a heterogeneous ”optical” medium, or ”crocodiles in the mist”. Tectonophysics, 232, p. 281– 297. Lyngsie, S. B., Thybo, H., and Lang, R., 2007. Rifting and lower crustal reflectivity: A case study of the intracratonic Dniepr-Donets rift zone. Ukraine. Journal of Geophysical Research, 112, B12402, doi:10.1029/2006JB004795.  159  Mandelbrot, B., 1983. The fractal geometry of nature. W.H.Freeman and company, New York, 468 p. Marsan, D., and Bean, C. J., 1999. Multiscaling nature of sonic velocities and lithology in the upper crystalline crust: Evidence from the KTB Main Borehole. Geophysical Research Letters, 26, p. 275-278. Nielsen, L., and Thybo, H., 2006. Identification of crustal and upper mantle heterogeneity by modeling of controlled-source seismic data. Tectonophysics, 416, p. 209–228. Oueity, J., and Clowes, R. M., 2010. The nature of the Moho transition in NW Canada from combined near-vertical and wide-angle seismic reflection studies. Lithosphere, 2, p. 377-396. Poppeliers, C., 2007. Estimating vertical stochastic scale parameters from seismic reflection data: deconvolution with non-white reflectivity. Geophysical Journal International, 168, p. 769– 778. Poppeliers, C. and Levander, A., 2004. Estimation of vertical stochastic scale parameters in the earth’s crystalline crust from seismic reflection data. Geophysical Research Letters, 31, p. 19538–19541. Pullammanappallil, S., Levander, A., and Larkin, S., 1997. Estimation of crustal stochastic parameters from seismic exploration data. Journal of Geophysical Research, 102, p. 15,269– 15,286. Robertsson, J., Balnch, J., and Symes, W., 1994. Viscoelastic finite-difference modeling. Geophysics, 59, p. 1444 1456. Ryberg, T., Tittgemeyer, M., and Wenzel, F., 2000. Finite difference modeling of p-wave scattering in the upper mantle. Geophysical Journal International, 141, p. 787–800. Scheirer, D. S., and Hobbs, R. W., 1990. Seismic attenuation in the continental crust SW of England. Geophysical Journal International, 103, p. 553-540. Thybo, H., 2006. The heterogeneous upper mantle low velocity zone. Tectonophysics, 416, p. 5379. Thybo, H., and Nielsen, C. A., 2009. Magma-compensated crustal thinning in continental rift zones. Nature, 457, p. 873-876, doi:10.1038/nature07688.  160  Tittgemeyer, M., Wenzel, F., Ryberg, T., and Fuchs, K., 1999. Scales of heterogeneities in the continental crust and upper mantle. Pure and Applied Geophysics, 156, p. 29-52. Villeneuve, M.E., Theriault, R.J., and Ross, G.M., 1991. U-Pb ages and Sm-Nd signature of two subsurface granites from the Fort Simpson magnetic high, northwest Canada. Canadian Journal of Earth Sciences, 28, p. 1003-1008.  161  Chapter 6  Conclusions In this Chapter we summarize key contributions of this thesis and discuss potential limitations. We also outline future research opportunities to further develop ideas presented in the thesis.  6.1 Main contributions The two main topics of this thesis are the characterization of the structure and evolution of the Moho transition zone and the geodynamical processes associated with its formation, as well as the nature and extent of the Paleoproterozoic relict subducted slab in the Wopmay orogen of Canada’s Northwest Territories. In our research we make use of the high quality, coincident refraction/wide-angle reflection and near-vertical incidence seismic reflection data recorded along corridor 1 of LITHOPROBE’S SlaveNorthern Cordillera Lithospheric Evolution (SNORCLE) transect. These two datasets represent a unique opportunity, at least in North America, to tackle the qualitative and quantitative aspects of the Moho transition zone and the subcrustal reflections in the region. Our approach includes 1- and 2-D forward and inverse traveltime modeling, 2-D finite-difference wave propagation algorithm, seismic attribute analysis and stochastic modeling based on the von Karman autocorrelation function. The reported results show that the methods used are well suited for our purposes and their applications can provide us with a better understanding about the complex nature of the lower crust, the Moho transition zone and the uppermost mantle. 6.1.1 Shallow relict subduction zone We use 2-D forward ray tracing and 3-D forward and inverse traveltime inversion methods to identify the extent, geometry and nature of a previously interpreted relict subducted slab beneath the Wopmay orogen in the Northwest Territories (Cook et al., 1999). The Paleoproterozoic relict slab is the result of shallowly east-dipping subduction 162  of oceanic lithosphere, which once separated the Fort Simpson and Hottah terranes, beneath the Hottah terrane. Although the angle of dip is very low, the slab is not regarded as a flat subduction. We demonstrate that the main characteristics found in modern flat subduction zones are missing. The slab extends over a lateral distance of ~300 km and between the 35 to 100 km depth interval where its seismic signature disappears. This can be explained either as the result of slab delamination in which any remainder of the slab has sunk into the mantle, or by reaching a level at which metamorphic minerals were converted to higher grade and became seismically indiscernible from mantle rocks. The subduction zone was active between 1.88-1.84 Ga and is probably responsible for the formation of the Great Bear magmatic arc. The deep upper mantle reflections imaged beneath the slab at around 100 km depth may represent a petrological transition from spinel-garnet peridotite to garnet-only peridotite (Kopylova et al., 1998) or a subcrustal imbricate structure (Helmstaedt and Schulze 1989). We also highlight the implications of projecting a crooked 3-D seismic line onto a 2-D cross section and the resulting artifacts arising in reflection images, a procedure that is sometimes overlooked and may lead to misinterpretation. 6.1.2 Application of curvelet denoising Here we advocate the novel use of curvelet-based denoising technique in removing incoherent noise present in deep crustal seismic reflection data where the signal-to-noise ratio (SNR) is low. This technique was previously used for enhancing industry-type reflection data but never for crustal data (Neelamani et al., 2008; Herrmann and Hennenfant, 2008). We apply the iterative curvelet denoising method to both pre-stack shot gathers and post-stack data recorded along SNORCLE Line 1 in the Northwest Territories. Our results demonstrate the effectiveness of this procedure in removing incoherent noise in both datasets as well as some of the coherent noise such as ground roll. Unlike the F-X deconvolution and the slant stack methods (Abma and Claerbout, 1995; Milkereit, 1987) that both work in a windowed sense, assume linear events and predict single dips at a time, the curvelet denoising can handle multiple dips and more complex and curved structures. Comparing the results for the stacked data with those from F-X deconvolution, curvelet denoising outperforms the latter by attenuating 163  incoherent noise with minimal harm to the signal. The enhanced clarity of deep crustal reflectivity due to the application of curvelet denoising allows for improved seismic interpretation. Hence, this technique can be regarded as a substantial new addition to the processing toolbox of crustal seismic reflection data. 6.1.3 The nature of the Moho transition Although the year 2009 marked one century since the discovery of the Moho, its exact nature remains one of the major uncertainties in lithospheric studies. In the past years a number of studies in Canada have addressed the issue of the Moho transition zone; however, they all followed a qualitative approach (e.g., Hammer and Clowes, 1997; Cook, 2002; Eaton, 2006; Cook et al., 2010). To quantitatively characterize the detailed structure of the Moho transition, we simulate the scattered wavefield (both reflected and refracted) by the calculation of both near-vertical and wide-angle synthetic seismograms for a number of velocity models using 1- and 2-D wave propagation algorithms. Only a combined interpretation of the two datasets can provide a more complete explanation regarding the structure of the Moho. The results indicate that only synthetic seismograms resulting from an approximately 3 km thick, laterally and vertically heterogeneous transition zone best fit the observed wavefield recorded on both data sets. The heterogeneity is represented by either laterally discontinuous layers or a lamellae structure. The discontinuous layers have velocities characteristic of both lower crust and upper mantle. On the other hand, the lamellae structure is described by randomly distributed ellipses embedded within a deterministic background velocity. The ellipses are characterized by their lateral (1000-3000 m) and vertical (100-300 m) scale lengths. However, these models are proposed to explain the seismic data observed beneath the Great Bear magmatic arc and may not be applicable in other parts of the world. 6.1.4 The stochastic lower crust, the Moho and the uppermost mantle Following a relatively new statistical analysis approach, we estimate the stochastic power law parameters directly from the near-vertical seismic data. This was possible due to the direct relation existing between the reflection data and the heterogeneity present in causative geological structure (Gibson, 1991; Holliger et al, 1992; Hurich, 1996; 164  Pullammanappallil et al, 1997; Hurich and Kocurko, 2000; Poppeliers, 2007). In turn, this allows us to construct a range of stochastic velocity models representative of the lower crust, the Moho transition zone and the uppermost mantle. To constrain the stochastic parameters, we take a step further than previous studies and simulate the wide-angle and near-vertical data from the stochastic velocity models using the finite-difference wave propagation algorithm. Our results strongly support the hypothesis that the lower crust, the Moho transition zone and the uppermost mantle all have a von Karman heterogeneity distribution. The Moho transition zone is characterized by the highest lateral correlation length with a value of 936 m, which we associate with mafic intrusions from the upper mantle. The transition zone separates the lower crust with a lateral correlation length of 732 m from the uppermost mantle with the lowest lateral correlation length of 261 m. We attribute the low correlation length in the mantle to the vertical migration paths of magma that may interrupted the horizontal lineation. 6.1.5 The origin and tectonic evolution of the Moho Based on our results from synthetic seismic modeling and statistical analyses, a combination of two geodynamical processes can explain the origin of Moho formation in the study area: (1) a regional décollement and (2) magmatic intrusions and partial melting. We offer a two-stage scenario describing the tectonic evolution of the Moho. In the first stage the Moho acted as a detachment surface to accommodate the compressional movement due to the accretion of the Hottah terrain onto the Slave craton ca. 1.9 Ga. In the second stage, the Moho developed as a thermal boundary due to the generation of partial melts and the placement of mafic intrusions at the base of the crust that were associated with the deposition of the Great Bear magmatic arc between 1.88 and 1.84 Ga.  6.2 Limitations Most of the limitations encountered throughout the course of this thesis can be attributed to two main sources; (1) the coverage and quality of the datasets themselves and (2) the computational restrictions. For example in Chapter 2 and based on the 3-D finitedifference traveltime inversion, our preferred model is a continuously dipping slab from  165  around 35 km (Moho level) to ~100 km where its seismic signature disappears. The accuracy of the model is affected by the coverage extent of depth points (midpoints between the shots and receivers) that sample the slab. In our case, these points are limited to the crooked 2-D geometry of SNORCLE Line 1. We interpret the apparent flattening of the slab as an artifact resulting from the projection of a 3-D structure onto the 2-D cross section. From the forward traveltime calculations, two models appear to produce synthetic reflections that fit the observed data well. In Model 1, the slab is represented by a dipping plane that flattens for several kilometers before dipping back again. On the other hand, Model 2 represents a continuously dipping slab. Without more complete 3-D seismic coverage, the possibility of a slab that flattens out cannot be dismissed. In Chapters 4 and 5, examination of the detailed structure of the Moho transition zone requires two important components; (1) coincident and (2) high quality near-vertical and wide-angle seismic reflection data. Successful and conclusive comparisons between synthetic and observed data, as well as the estimation of the stochastic parameters from the near-vertical data, which describe the von Karman heterogeneity distribution, may not be possible in the absence of these components. The SNORCLE Line 1 data are exceptionally high quality. Therefore, the applicability of our forward and inverse approach to further constrain the nature of the Moho transition zone, and to a certain extent the lower crust and uppermost mantle, has not been confirmed for the case where more normal, poorer quality data are involved. Moreover, some of the results we present in this thesis are subject to our own interpretation. For instance, in Chapter 4 our choice for representing the Moho transition zone by either laterally discontinuous layers or randomly distributed ellipses within a zone 3 km thick is subjective and other models that can still explain the observed data are possible as well. This follows from the forward modeling approach that is used (Chapter 4); it is a non-unique procedure with a subjective interpretation. Accordingly, we apply a more objective inverse approach involving statistical analyses of the near-verticalincidence reflection data (Chapter 5). Given the data set that we are using, the results are positive and reasonable. However, what type of modality (Gaussian, bi-modal, tri-modal,  166  multi-modal) that is best suited to map the von Karman velocity fields is still a matter of debate. Due to the strong and well-defined Moho reflections beneath the Great Bear magmatic arc, all the models we proposed for the transition zone in Chapters 4 and 5 are limited to that specific region. However, the Moho signature along SNORCLE Line 1, as well as other lines in the SNORCLE program (Cook et al., 2003), varies significantly from one place to another. For example, it changes from a multicyclic band of reflections beneath the Great Bear arc to a complete absence of reflections beneath the Fort Simpson-Hottah terrane. On the other lines, which are located within the Canadian Cordillera, the Moho signature varies from a distinguishable reflector through a prominent decrease in reflectivity to diffuse reflectivity or little reflectivity at all. Thus, our results probably are not representative of the general nature of the continental Moho. Many more similar studies would have to be undertaken to illuminate the varying complexity of this important lithospheric boundary. In terms of verifying interpretations of crustal/upper mantle reflection and refraction data, a severe limitation is the fact that no “ground-truthing” is possible. Sampling and analysis of the continental deep crust, Moho and uppermost mantle by drilling and logging is not feasible; temperatures are simply too hot. This contrasts with the use of the seismic reflection method in petroleum exploration or near-surface studies for which drilling and logging results are available to complement and check the seismic interpretations. Finally, the structural features and inherent heterogeneity that we interpret for the lower crust, the Moho transition zone and the uppermost mantle in Chapters 4 and 5 are three-dimensional in reality. We try to accommodate this fact with various approximations. However, we are restricted to two-dimensional full wavefield calculations due to computational limitations. Unfortunately, calculation of the wavefield in large-scale three-dimensional heterogeneous media, such as the lower crust and upper mantle, is not practically possible on present-day computers that are available to the scientific community.  167  6.3 Future work 6.3.1 Application of curvelet denoising We have clearly demonstrated in Chapter 3 the efficacy of the curvelet denoising algorithm for removing unwanted random noise from pre- and post-stack reflection data. However, the technique has only been used on a small segment (60 km) of SNORCLE Line 1 and was limited to crustal depths (~34 km). The application of curvelet denoising to the entire 725-km-long seismic line and down to upper mantle levels (~90 km) may result in a revised and improved interpretation, especially for the abundant upper mantle reflections. For example, the denoising technique may help identify any weak upper mantle reflectors beneath the Hottah-Slave boundary. This would assist investigation of any geometric relation between the relict subducted slab and the anisotropic layered structure imaged by the teleseismic data beneath the Slave craton (Bostock, 1998). Another set of east-dipping upper mantle reflections east of the subducted slab (Cook et al., 1999), observed as a series of discontinuous events, can be traced from Moho levels to about 70 km depth. The interpretation of these reflectors remains speculative because of their ambiguous geometric relation with the Moho due to conflicting dips and poor image quality.  This again represents an opportunity for the application of curvelet  denoising. More generally, the approach may be useful on both legacy and new reflection data worldwide to provide higher quality images that may enable better interpretations. 6.3.2 Application of quantitative analysis By means of combining the results of synthetic wide-angle and near-vertical seismograms along with the estimation of stochastic parameters directly from nearvertical data, we have established a solid approach to quantitatively characterize the nature of deep crustal and uppermost mantle reflections, particularly in the Moho transition zone and surrounding regions. Along SNORCLE Line 1, the apparent reflectivity of the Moho transition imaged by the near-vertical data varies considerably. Reflectivity patterns range from a strong, multicyclic band of reflectors that are piecewise continuous over many kilometers, to a complete absence of reflectors. In addition,  168  strong Moho reflections on the near-vertical data are not necessarily accompanied by distinct wide-angle PmP phase and vise versa. For example, beneath the Hottah-Fort Simpson transition zone the near-vertical data do not show any Moho reflections, whereas the corresponding wide-angle data exhibit a distinct reverberatory PmP phase. The reason for such observations is not clear. One explanation is that fluids released by the dehydration of the subducted slab could cause serpentinization, which in turn would lower the acoustic impedance contrast between lower crust and upper mantle (Fernandez Viejo and Clowes, 2003). Another explanation is the presence of pyroxenite, which has a similar effect on the acoustic impedance (Fernandez Viejo et al., 2005). The application of quantitative analysis would probably help us better understand the causes of the variable reflectivity patterns and the discrepancy observed occasionally between nearvertical and wide-angle data. Our approach cannot only be applied to SNORCLE line 1, but also to other SNORCLE lines and LITHOPROBE transects where high quality nearvertical and wide-angle datasets are available, and similar studies elsewhere. For LITHOPROBE, the seismic data archive includes an extensive collection of such data; more information can be found at http://gdr.nrcan.gc.ca/seismtlitho/archive/index_e.php. 6.3.3 More and better data leading to 3-D surveys During acquisition of refraction/wide-angle reflection data along SNORCLE Line 1 in the Northwest Territories in 1997, nearly 600 individual seismographs were deployed at an average spacing of 1.2 km along the 725-km profile (Chapter 3; Fernandez Viejo and Clowes, 2003). Thirteen shot points with typical spacing of 50 to 60 km provided the energy source. During acquisition of a similar survey in central British Columbia in 2009, an experiment in which this writer participated, more than 2000 individual seismographs were deployed at an average spacing of 0.2 km along a 400-km profile. Sixteen shots were detonated in fourteen locations to generate the seismic energy, giving an average shot spacing of about 30 km. The increase in density of seismographs and shot points is important because this enables much higher resolution in the interpreted models. Indeed, this experiment was approaching a merger of R/WAR techniques and coarse NVI methods. As seismographs become more sophisticated and less expensive, a fully 3-D lithospheric survey in which 10,000 or more seismographs are deployed looms as a 169  distinct possibility. Unfortunately, cost and permitting issues for the necessary explosive energy sources for R/WAR surveys will continue to limit a similar increase in the number of shots and decrease in shot spacing. Nevertheless, a focused 3-D survey is possible. The SNORCLE Line 1 NVI vibroseis reflection survey, carried out in 1996, recorded 404 channels with a receiver spacing of 60 m and a shot spacing of 90 m (Chapter 3; Cook et al., 1999). A similar vibroseis reflection survey, carried out in central British Columbia in 2008, recorded 960 channels with a receiver spacing of 20 m and a shot spacing of 40 m. Again, the enhanced parameters allow for a higher resolution image to be generated. However, through new technological developments, industry can now acquire 10,000 or even 20,000 channels of data for their 3-D surveys (information obtained at GeoCanada 2010, May 2010 in Calgary). If adequate funds were available, similar experiments could be carried out for lithospheric studies using a larger receiver and shot spacing than is typical in industry. Given the effort and expense in carrying out such enhanced field experiments for lithospheric studies, all possible analysis techniques would need to be applied to the acquired dataset to maximize the scientific output. This thesis describes some of those analysis techniques, albeit mainly in a 2-D implementation. The 3-D versions would require further development of the algorithms.  170  Bibliography Abbott, D. 1991. The case for accretion of the tectosphere by buoyant subduction. Geophysical Research Letters, 18, p. 585– 588. Abma, R., and Claerbout, J., 1995. Lateral prediction for noise attenuation by t-x and f-x techniques. Geophysics 60, p. 1887–1896. Bostock, M.G. 1998. Mantle stratigraphy and evolution of the Slave province. Journal of Geophysical Research, 103, p. 183-200. Cook, F.A., 2002. Fine structure of the continental reflection Moho. Geological Society of America Bulletin, 114, p. 64-79. Cook, F.A., van der Velden, A.J., and Hall, K.W. 1999. Frozen subduction in Canada's Northwest Territories: Lithoprobe deep lithospheric reflection profiling of the western Canadian Shield. Tectonics, 18, p. 1-24. Cook, F. A., White, D. J., Jones, A. G., Eaton, D. W., Hall, J., and Clowes, R. M., 2010. How the crust meets the mantle: Lithoprobe perspectives on the Mohorovičić discontinuity and crustmantle transition. Canadian Journal of Earth Sciences, 47, p. 315-351. Eaton, D. W., 2006. Multi-genetic origin of the continental Moho: insights from Lithoprobe. Terra Nova, 18, p. 34-43, doi: 10.1111/j.1365-3121.2005.00657. Fernandez Viejo, G., and Clowes, R.M., 2003. Lithospheric structure beneath the Archean Slave Province and Proterozoic Wopmay orogen, northwestern Canada, from a Lithoprobe refraction/wide-angle reflection survey. Geophysical Journal International, 153, p. 1-19. Fernandez Viejo, G., Clowes, R.M., Welford, J. K., 2005. Constraints on the composition of the crust and uppermost mantle in northwestern Canada: Vp/Vs variations along Lithoprobe's SNorCLE transect. Canadian Journal of Earth Sciences, 42, p. 1205-1222, doi:10.1139/e05028. Gibson, B., 1991. Analysis of lateral coherency in wide-angle seismic images of heterogeneous targets. Journal of Geophysical Research, 96, p. 10,261–10,273. Hammer, P. T. C., and Clowes, R.M., 1997. Moho reflectivity patterns - a comparison of Canadian LITHOPROBE transects. Tectonophysics, 269, p. 179-198.  171  Helmstaedt, H., and Schulze, D. J. 1989. Southern African kimberlites and their mantle sample: Implications for Archean tectonics and lithosphere evolution. In Kimberlites and Related Rocks, Blackwell , Cambridge, Mass., p. 358– 368. Herrmann, F. J., and Hennenfent, G., 2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International, 173, p. 233–248. Holliger, K. and Levander, A.,1992. A stochastic view of lower crustal fabric based on evidence from the Ivrea Zone. Geophysical Research Letters, 19, p. 1153–1156. Hurich, C., 1996. Statistical description of seismic reflection wavefields: a step towards quantitative interpretation of deep seismic reflection profiles. Geophysical Journal International, 125, p. 719–728. Hurich, C. and Kocurko, A., 2000. Statistical approaches to interpretation of seismic reflection data. Tectonophysics, 329, p. 251–267. Kopylova, M. G., Russell, J. K., and Cookenboo, H., 1998. Upper-mantle stratigraphy of the Slave craton, Canada: Insights into a new kimberlites province. Geology, 26, p. 315-318. Milkereit, B., 1987. Migration of noisy crustal seismic data. Journal of Geophysical Research 92, p. 7916-7930. Neelamani, R., Baumstein, A.I., Gillard, D.G., Hadidi, M.T., and Soroka, W.L. , 2008. Coherent and random noise attenuation using the curvelet transform. The Leading Edge, 27, p. 240-248. Poppeliers, C., 2007. Estimating vertical stochastic scale parameters from seismic reflection data: deconvolution with non-white reflectivity. Geophysical Journal International, 168, p. 769-778.  172  Appendix A  Wave propagation algorithms A.1 One dimensional reflectivity method The reflectivity method was developed by Fuchs and Müller (1971). In this approach the numerical integration of the reflectivity (or plane wave reflection coefficient) of a layered medium is carried out in the horizontal wavenumber or angle of incidence domain. The earth model for which the reflectivity method is developed consists of n homogeneous and isotropic layers on top of a half-space. The Fourier transform of the compressional potential of the wave from the explosive point source is given by:  Φ (r , z, w ) = F (w )∫  ∞ 0  k iν  j 0 ( kr ) exp( − i ν z ) dk  (A.1.1)  where F (w) is the Fourier transform of the excitation function f (t ) and j0 (kr ) is the Bessel function of the first kind and zero order, i the imaginary number, k the horizontal wavenumber and ν is the vertical wave number. The first step is to transmit the P wave through the layers until it reaches the reflective interface, reflects and then is transmitted back across the interfaces to the surface. The resulting Fourier transformed potential has the same form as (A.1) but also includes transmission and reflection coefficients: ∞  ~ k j 0 (kr ) Pd ( w, k ) R pp ( w, k ) Pu ( w, k ) × iv1 0  Φ ( r , z , w) = F ( w) ∫  (A.1.2)  m ⎡ ⎤ exp ⎢− i (2∑ hi vi + zv1 ) ⎥ dk i =1 ⎣ ⎦  where Pd ( w, k ) and Pu ( w, k ) are the product of the transmission coefficients for an ~  downgoing and upgoing wave respectively, and R pp ( w, k ) are the product of reflection coefficients. 173  In the next step, the vertical and horizontal displacement at the free surface is determined from the potential of the reflected P and S waves, where the relations between the vertical/horizontal displacement and potential are given by: u z =  ∂Φ ∂Φ and u r = ∂z ∂r  ,  respectively. The final step is changing the variable of integration k into a new variable defined as the angle γ. As long as γ is real, it can be interpreted as the angle of incidence at the top of the reflecting zone. Inverse Fourier transformation results in the synthetic seismogram for the displacement component.  A.2 Two dimensional visco-elastic finite difference method The finite difference scheme simulates wave propagation in viscoelastic media (Robertsson et al., 1994). The basic hypothesis in viscoelastic theory is that the current value of the stress tensor σ depends on the history of the strain tensor ε that can be described as: (A.2.1) where * denotes time convolution, the dot denotes derivation with respect to time and G is a fourth-order tensor called the relaxation function, which collapses into two independent functions for an isotropic and homogeneous material. For the 1-D case in an isotropic homogeneous material, equation (A.2.1) reduces to (A.2.2) The relaxation function G determines the behavior of the material and is given by:  (A.2.3) where MR is the relaxation modulus of the medium, θ(t) is the Heaviside function, and the elements τel and τσl are the stress and strain relaxation times of the lth mechanism. For simplicity, the equations describing wave propagation in viscoelastic media are derived here for the 1-D case, which can be generalized to higher dimensions.  174  Taking the time derivative of equation (A.2.3) and using  g  where p is  the pressure and v is the velocity, gives: (A.2.4) The convolution terms can be eliminated by introducing memory variables defined as rl. Equation (A.2.4) thereby reduces to:  (A.2.5) where  (A.2.6)  A set of first-order linear differential equations can be obtained as follows. First, by taking the time derivative of equation (A.2.6), we obtain: (A.2.7)  Newton’s second law completes the full description of wave propagation in a viscoelastic medium; (A.2.8)  Equations (A.2.5), (A.2.7), and (A.2.8) comprise the system of first-order linear differential equations of 1-D viscoelastic wave propagation.  175  Bibliography Fuchs, K ., and Müller, G., 1971, Computation of synthetic seismograms with the reflectivity method and comparison with observations, Geophysical Journal of Royal Astronomical Society, 23, p. 417-433. Robertsson, J., Balnch, J., and Symes, W., 1994, Viscoelastic finite-difference modeling: Geophysics, 59, p. 1444 1456.  176  Appendix B  Supplementary information for Chapter 5 B.1 Estimation of power law parameters First, the normalized 2-D autocorrelation in an N x M window of seismic data is computed (Carpentier, 2007):  (B.1) where p(xi, zj) is a discrete 2-D seismic signal, l is the autocorrelation lag in the x-direction and τ is the autocorrelation lag in the z-direction (or time direction). From this 2-D matrix, we select the first row  as the average 1-D lateral (horizontal) correlation curve for the  entire N x M window. The procedure is illustrated by Fig (B.1).  Figure B. 1. Illustration of the computation of the 2-D autocorrelation for a window of data p(x,z). The red, blue, and yellow squares represent its horizontal, vertical and diagonal components, respectively. We select the horizontal component (first row is First column is  , red arrow) as an average 1-D autocorrelation for the whole window.  , blue arrow (from Carpentier, 2007).  177  Next, the observed  curve is fitted with a 1-D von Karman autocorrelation function  as an estimator for the average horizontal von Karman content of the window, in terms of parameters ax (horizontal correlation length) and v (Hurst power). For this purpose, the righthand-side of the expression (B.2) has to be minimized (Press et al., 1992) :  , where  (B.2)  is the least-squares misfit as a function of ax and v, are the observed datapoints of  and  are the predicted points of the 1-D von Karman function with parameters ax and v. Parameters ax and v that best predict the observed data are supplied by a gridded modelparameter-space search. A typical fitting result is shown in Fig (B.2).  Figure B. 2. Illustration of typical lateral correlation function (  , square dots plus solid line) and best fit  von Karman curve (dashed line), for the input window in Figure B.1 (from Carpentier, 2007).  178  Bibliography Carpentier, S. F. A., 2007, On the estimation of stochastic parameters from deep seismic reflection data and its use in delineating lower crustal structure. Ph. D. Thesis. University of Utrecht. Press, W., Teukolsky, S. A., Vetterling, W. T., and Flannery, B., 1992, Numerical recipes in C. The art of scientific computing. Second edition. Cambridge Univ. Press.  179  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0052966/manifest

Comment

Related Items