Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Lithospheric structure across the Canadian cordillera of Northeastern British Columbia from seismic refraction… Welford, Joanna Kim 1999

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

Item Metadata

Download

Media
831-ubc_2000-0309.pdf [ 42.22MB ]
Metadata
JSON: 831-1.0089503.json
JSON-LD: 831-1.0089503-ld.json
RDF/XML (Pretty): 831-1.0089503-rdf.xml
RDF/JSON: 831-1.0089503-rdf.json
Turtle: 831-1.0089503-turtle.txt
N-Triples: 831-1.0089503-rdf-ntriples.txt
Original Record: 831-1.0089503-source.json
Full Text
831-1.0089503-fulltext.txt
Citation
831-1.0089503.ris

Full Text

L I T H O S P H E R I C S T R U C T U R E ACROSS T H E C A N A D I A N C O R D I L L E R A O F N O R T H E A S T E R N BRITISH C O L U M B I A F R O M SEISMIC R E F R A C T I O N A N D P O T E N T I A L F I E L D D A T A By Joanna K i m Welford B. Sc. (Geophysics) McGil l University A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E STUDIES E A R T H A N D O C E A N SCIENCES We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH C O L U M B I A December 1999 © Joanna Kim Welford, 1999 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract The region of northeastern British Columbia includes the transition from cratonic North America to the deformed Canadian Cordillera. While geological mapping has progressed for over 120 years in the northern Cordillera, the subsurface has remained terra incognita as few geophysical surveys have been conducted. The SNoRE '97 seismic refraction/wide-angle reflection ( R / W A R ) survey was carried out to determine the velocity structure of the lithosphere in northwestern Canada. It sampled four separate seismic profiles along the three corridors of the L I T H O P R O B E SNorCLE transect. This thesis focuses on the results for Line 21. From east to west, the line begins in the Proterozoic Fort Simpson magmatic arc terrane, crosses the accreted Nahanni terrane, transects the Foreland Belt and samples the edge of the Omineca Belt just west of the Tintina Fault. The refrac-tion data set shows good signal-to-noise ratios; primary phases (Pg, PmP and Pn) are clearly observed for most shots. A velocity structural model is achieved first through r — p analysis, then forward ray trace modeling and 2-d inversion of traveltimes and for-ward amplitude modeling using R A Y I N V R , a widely applied, ray-theory-based algorithm developed for R / W A R studies. The resolved velocity structural model demonstrates a westward thickening package of low upper crustal velocities (6.2 km/s and less) extending to almost 20 km depth along the eastern portion of the profile below the W C S B . Beneath the westward thickening low velocities, a west-facing ramp of higher velocities, resolved in the middle and lower crust, is interpreted as the rifted margin of cratonic North America, which was overlain by low velocity, passive margin sediments during the mid-to-late Proterozoic and early Paleozoic. The eastern part of the Foreland Belt shows high velocities (6.4 km/s) in n the upper five kilometers, indicating rocks upthrust from deeper in the crust. At the western end of the line, a narrow trench of low velocities at the surface of the model is inferred to represent sediments and brecciated surface rocks filling in the scar of the crustal-scale Tintina Fault. The location of the fault is resolved ~ 20 km to the west of its current inferred location. The depth to Moho is ~ 34 km in the Omineca Belt, deepens to ~ 38 km below the eastern Foreland Belt, and shallows to 34 km at the Cordilleran deformation front. The upper mantle is laterally heterogeneous with anomalously high velocities (up to 8.3 km/s) beneath the Foreland Belt flanked by regions of low velocities (7.7-7.8 km/s). Results from 3-D Bouguer gravity inversion, 2.5-D Bouguer gravity forward mode-ling and aeromagnetic comparisons performed in this thesis are combined with results from heat flow analyses and previous geophysical studies to check the consistency of the resolved velocity model and to constrain the interpretation of the tectonic evolution of northeastern British Columbia. m Table of Contents Abstract ii List of Tables viii List of Figures ix Acknowledgement xiv Dedication xv 1 I N T R O D U C T I O N 1 1.1 Thesis Outline 2 1.2 Tectonic and Geological Background 2 1.2.1 The Canadian Cordillera and ancestral North America 2 1.2.2 Regional Tectonic Setting for Line 21 8 1.2.3 Geology of Line 21 12 1.3 Previous work in northeastern British Columbia 16 1.3.1 Geological Mapping 16 1.3.2 Geophysics in the Northern Canadian Cordillera 16 1.4 Experimental Objectives 21 2 D A T A ACQUISITION A N D P R O C E S S I N G 22 2.1 The 1997 Slave-Northern Cordillera Refraction Experiment . . 22 2.1.1 Location 23 2.2 Data Acquisition . . . 23 iv 2.2.1 Drilling and Loading of Source Locations 24 2.2.2 Surveying Shots and Receivers 25 2.2.3 Instrument Deployment 26 2.2.4 Recording Instruments 27 2.2.5 Timing of Shots and Receivers 30 2.3 Preparation of Data 31 2.3.1 Processing Routines 31 2.3.2 Gain 31 2.3.3 Updating of SEG-Y Headers 32 2.3.4 Merging Data and Creating Final Data Sets 32 2.4 Data Processing 33 2.4.1 True Relative Amplitude Scaling 33 2.4.2 Filtering . . 33 3 D A T A D E S C R I P T I O N A N D ANALYSIS 34 3.1 Data Quality 34 3.1.1 Observations of Pg Phase 34 3.1.2 Observations of Crustal Reflectors 37 3.1.3 Observations of PmP Phase 38 3.1.4 Observations of Pn Phase 39 3.1.5 Summary of phase observations 41 3.2 Data Interpretation and Picking 42 3.3 Modeling Algorithms and Techniques 46 3.3.1 Tau-p Analysis 47 3.3.2 R A Y I N V R Modeling 49 3.3.3 Amplitude Modeling 55 v 4 D E V E L O P M E N T O F T H E V E L O C I T Y M O D E L 57 4.1 Modeling the Refraction Data 57 4.1.1 Detailed modeling results 60 4.1.2 Spatial Resolution and Absolute Parameter Uncertainty 72 4.2 The Velocity Structural Model for Line 21 76 4.2.1 Main Features of the Velocity Model 79 4.2.2 Velocity Anomaly Calculation 82 5 P O T E N T I A L F I E L D D A T A M O D E L I N G 84 5.1 Gravity Data and Modeling 84 5.1.1 3-d Inversion Modeling using GRAV3D 86 5.1.2 2.5-d Forward Modeling using S A K I 96 5.2 Isostatic Balance Considerations . 100 5.3 Aeromagnetic Data Interpretation 105 5.3.1 Reduction to pole 106 5.3.2 Interpretation of Magnetic Anomalies 107 6 DISCUSSION 110 6.1 Comparison with other Geophysical Studies in the Region 110 6.1.1 Comparison with Refraction Studies in the Northern Cordillera . 110 6.1.2 Comparison with Refraction Studies in the Southern Cordillera . 119 6.1.3 Comparison with Reflection Studies in the Northern Cordillera . . 122 6.1.4 Comparison with Heat Flow Studies 124 6.2 Final Velocity Model within the Context of the Regional Tectonics . . . . 126 6.2.1 Tectonic Implications of Upper and Middle Crustal Results . . . . 126 6.2.2 Tectonic Implications of Lower Crustal and Upper Mantle Results 128 v i 6.2.3 Interpretation of Variations along-strike of the Cordillera-Craton Transition Zone 129 6.2.4 Interpreted Lithospheric Structure for Northeastern British Columbia 132 6.3 Conclusions 134 References 137 Appendices 148 A Data and Related Results 148 B Inverse Problem Mathematics 173 vii Lis t of Tables 2.1 Shot point information for Line 21. 24 2.2 Seismic Instrument Information 29 2.3 Geophone Information 29 2.4 Gain applied to each instrument-geophone combination 32 3.1 Source-receiver offsets of identified phases 44 3.2 Number of observations for each shot point with corresponding average traveltime uncertainties 45 4.1 Results of spatial resolution and absolute parameter uncertainty tests. . . 75 4.2 Model fitting statistics for each phase 77 4.3 Model fitting statistics for each shot point 79 vm List of Figures 1.1 SNorCLE transect of L I T H O P R O B E 3 1.2 The morphogeological belts of the Canadian Cordillera 6 1.3 Terrane map of SNorCLE transect surrounding Line 21 9 1.4 Diagram of Proterozoic Wopmay Orogen along the margin of cratonic North America 10 1.5 Geological map for the region surrounding SNoRE '97 Line 21 13 1.6 Structural paleogeographic feature map for Line 21 14 1.7 Cross-section of surface sediments 15 1.8 West-facing ramp of cratonic North America's edge 18 3.1 Trace normalized plot for shot 2108 showing slowed Pg phase 36 3.2 Example data plots where Pg amplitudes appear muted 37 3.3 Trace normalized plot for shot 2107 showing examples of the PmP phase. 39 3.4 Trace normalized plot for shot 2102 showing an example of the Pn phase. 40 3.5 True relative amplitude plot for shot 2108 41 3.6 Initial r - p Velocity Model for Line 21 49 3.7 Schematic diagram of R A Y I N V R parameterization 51 4.1 Flow chart for modeling of refraction data 58 4.2 Nodal parameterization of final velocity model for Line 21 61 4.3 Ray coverage for Pg phase 63 4.4 Ray coverage for R phases 66 4.5 Ray coverage for PmP phase 67 ix 4.6 Ray coverage for R l phase . 68 4.7 Ray coverage for Pn phase 69 4.8 Comparison of observed and calculated traveltimes from shot 2105 . . . . 70 4.9 Trace normalized plot for shot 2101 showing PmP phase masked by earlier R l arrivals. . 71 4.10 Final velocity structural model for Line 21 78 4.11 Ray coverage of velocity model for Line 21 80 4.12 Velocity anomaly model for Line 21 83 5.1 Bouguer gravity map of region surrounding Line 21 . . . 85 5.2 Averaged observed Bouguer gravity map with regional removed 88 5.3 Topographical information for the Line 21 region 89 5.4 Initial predicted Bouguer gravity map 90 5.5 Slice through initial density anomaly model at 300 km northing 91 5.6 Reference density anomaly model used in the inversion 91 5.7 Comparison of observed and predicted Bouguer gravity maps from the GRAV3D inversion 92 5.8 Depth Slices through the density anomaly model generated by G R A V 3 D 94 5.9 Slice through GRAV3D density anomaly model along modeled refraction line 95' 5.10 Comparison of gravity trends along parallel lines. 98 5.11 2.5-d Density model derived for Line 21 99 5.12 Model for isostatic balance calculations. • • • 102 5.13 Results from calculations of isostatic compensation 104 5.14 Aeromagnetic map of region surrounding Line 21 106 5.15 Comparison of aeromagnetic data before and after reduction to pole . . . 108 x 5.16 Close-up of reduced to pole aeromagnetic map at the Yukon/BC border. 109 6.1 Comparison of velocity models from Line 11 and Line 21 114 6.2 Combination of velocity models for Line 21 and Line 22 117 6.3 Velocity models for lines 21 and 22 projected perpendicular to geological strike 118 6.4 Comparison between interpreted reflection results from Line 11 and veloc-ity model for Line 21 123 6.5 Interpreted Foreland-craton transitions for the southern and northern Cor-dilleras 130 6.6 Interpreted evolution of Foreland-craton transition in the northern Cordilleral31 6.7 Interpreted cross-section for the Cordilleran-North American craton tran-sition 133 A.1 Shot 2101. a) Trace normalized plot. Expected location of PmP arrivals are indicated in the coda of the R l arrivals (see text for discussion); on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed and calculated traveltimes 149 A.2 True relative amplitude plot for shot 2101 150 A.3 Shot 2101. a) Synthetic seismogram record section; b) A V O plot for ob-served and synthetic plots 151 A.4 Shot 2102. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes . • . 152 A.5 True relative amplitude plot for shot 2102 . 153 A.6 Shot 2102. a) Synthetic seismogram record section; b) A V O plot for ob-served and synthetic plots 154 xi A.7 Shot 2104. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes 155 A.8 True relative amplitude plot for shot 2104 - . . . . 156 A.9 Shot 2104. a) Synthetic seismogram record section; b) AVO plot for ob-served and synthetic plots 157 A. 10 Shot 2105. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes 158 A.11 True relative amplitude plot for shot 2105 159 A. 12 Shot 2105. a) Synthetic seismogram record section; b) A V O plot for ob-served and synthetic plots 160 A.13 Shot 2106. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes 161 A. 14 True relative amplitude plot for shot 2106 162 A.15 Shot 2106. a) Synthetic seismogram record section; b) AVO plot for ob-served and synthetic plots 163 A. 16 Shot 2107. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes 164 A.17 True relative amplitude plot for shot 2107 165 A. 18 Shot 2107. a) Synthetic seismogram record section; b) AVO plot for ob-served and synthetic plots 166 A. 19 Shot 2108. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes 167 A.20 True relative amplitude plot for shot 2108 168 A.21 Shot 2108. a) Synthetic seismogram record section; b) AVO plot for ob-served and synthetic plots 169 xi i A.22 Shot 2109. a) Trace normalized plot; b) Comparison of observed and calculated traveltimes . 170 A.23 True relative amplitude plot for shot 2109 171 A.24 Shot 2109. a) Synthetic seismogram record section; b) A V O plot for ob-served and synthetic plots 172 xm Acknowledgement I wish to thank my supervisor Dr. Ron M . Clowes for his support and guidance over the past two years. Despite his hectic schedule as Lithoprobe Director, his many hours of dedication to my cause, his unwavering faith in my abilities and his ability to always be "cool" have made this degree tremendously enjoyable. Through his insightful criticism of my work and his attention to detail, he has taught me how to be a better scientist. Thank you also to Dr. Robert M . Ellis and Dr. James K . Mortensen for critically reviewing this thesis and to Dr. Tadeusz J . Ulrych for ensuring that my defense ran smoothly and humorously. To Phil Hammer, Andrew Gorman, Gabriela Fernandez Viejo, Baishali Roy and Barry Zelt, I would like to offer my thanks for putting up with my incessant seismic questions and for acting as my surrogate advisors when Ron was out of town. A huge thank you to John Amor for dealing with all my computational and technical "issues" and for bailing me out of many minor catastrophes. To everyone in the Geophysics Branch of the Department of Earth and Ocean Sciences, I would like to offer my thanks for making this such a great place to spend way too much of my time. I'm truly blessed to have had the opportunity to work with such brilliant minds. A huge thank you to Colin Farquharson for listening to all my "rants" over the past two years, for putting up with all my inverse theory questions over dinner and for always knowing how to make me laugh. Lastly, I would like to thank my parents for providing me with the drive and self-confidence to pursue and achieve my goals. xiv To my parents and to Colin, home away from home xv Chapter 1 I N T R O D U C T I O N The dramatic topographic transition from the low plains of the Western Canada Se-dimentary Basin to the mountainous vistas of the Canadian Cordillera has fascinated geoscientists for over a century. This fascination has led to extensive geological mapping and geophysical surveying of the accessible southern Canadian Cordillera with the goal of better understanding the evolution of the North American Cordillera and the interaction between ancestral North America and the Cordillera. These studies have provided de-tailed information about the geology, fossil history, velocity structure, rheological fabric, potential field signature, thermal regime and tectonic evolution of the crust and upper mantle in the southern Canadian Cordillera. However, due to monetary constraints, rugged terrain and remoteness, research in the northern Canadian Cordillera has until recently only been pursued by a small number of geoscientists. Through L I T H O P R O B E , Canada's national geoscience multidisciplinary research project, attention recently has been focused on this region through the SNorCLE (Slave-Northern Cordillera Lithosphere Evolution) transect. The project has enabled a large variety of surveys from many geoscientific disciplines to be collected in this relatively remote portion of Canada. This thesis focuses on the processing, modeling and inter-pretation of Line 21 of the extensive refraction/wide-angle reflection data set collected during the SNorCLE Refraction Experiment, SNoRE '97, and the integration of these results with other geophysical and geological data. 1 Chapter 1. INTRODUCTION 2 1.1 Thesis Outline Chapter 1 describes the tectonic setting and summarizes previous geophysical studies of the Canadian Cordillera and cratonic North America focusing, in particular, on those aspects pertaining to Line 21 of SNoRE '97. Chapter 2 discusses the SNoRE '97 refraction experiment and presents the specific experimental details for Line 21. Chapter 3 discusses the data associated with the experiment and the detailed methodology used to process and model the data collected along Line 21. Chapter 4 discusses the modeling results and presents the final velocity model. To provide further constraints on the interpretation of the final velocity model, potential field data analyses were carried out; Chapter 5 describes the results. Chapter 6 presents a discussion of the results of the thesis along with documentation supporting these results and offers some final conclusions. Appendix A contains all relevant data and analysis plots. 1.2 Tectonic and Geological Background 1.2.1 The Canadian Cordillera and ancestral North America The study area for this thesis lies in northeastern British Columbia and spans, from east to west, the western edge of cratonic North America and the eastern portion of the Canadian Cordillera (Fig. 1.1). The diverse terranes making up the region span geological history from Proterozoic to Phanerozoic times and provide the opportunity to examine how the Canadian Cordillera was accreted to and deformed the western edge of cratonic North America. The western margin of ancestral North America, Laurentia, developed into a rifted margin in Proterozoic time and was overlain by a passive margin miogeoclinal sequence that spanned from the Middle Proterozoic (1500 Ma) to the Middle Devonian (380 Ma) F i g u r e 1.1: S N o r C L E transect of L I T H O P R O B E . The numbers 11, 21, 22 and 31 correspond to the four refraction lines collected during the S N o R E '97 Exper iment . A l l 36 shots of the experiment are shown as stars. The dashed box outlines the project ion used for a l l subsequent regional maps. Inset map shows locat ion of the S N o r C L E Transect w i t h respect to N o r t h A m e r i c a . (Monger and Price, 1979). Following the rifting event, the miogeocline was affected by almost continuous subsidence during the Cambrian and Ordovician (Cecile et al., 1997). The cratonic margin demonstrates a continental shelf to slope succession which was affected by at least 3 phases of crustal extension since its formation (Thompson et al., 1987). Evidence of these phases of extension exists in two of the main stratigraphic sequences of the Canadian Cordillera (sequences A and C; Hoffman (1989), Aitken and McMechan (1991), Gabrielse and Campbell (1991)) . Chapter 1. INTRODUCTION 4 Sequence A , deposited between 1.7 and 1.2 Ga, is primarily associated with the Belt-Purcell Supergroup which is exposed in southeastern British Columbia, northeastern Washington State and northwestern Montana. It is composed of a thick accumulation (11 to 20 km) of clastic sedimentary rocks, sandstones and carbonates interlayered with minor mafic volcanic sills which provide age constraints. Its base is exposed in Montana where it is found to overlie Precambrian basement rocks. Due to the fan-shaped depositional distribution of the Belt-Purcell Supergroup, it has been interpreted as an intracratonic basin that developed in the southern Cordillera at the western edge of the preserved North American craton (Hoffman, 1989; Aitken and McMechan, 1991). The isotopic signature of the sediments and paleocurrent directions suggest that the sediments originated from the west from another continental land mass that lay adjacent to western North America at the time. This continental land mass is widely interpreted to have been southern Australia (Ross et al., 1992). , In the northern Cordillera, Sequence A is exposed in two regions. Firstly, in north-eastern British Columbia to the northeast of the Northern Rocky Mountain Trench, it is exposed as the 6 km thick Muskwa assemblage which unconformably underlies Sequence C rocks. Secondly, in the Yukon Territory, Sequence A exists as the Wernecke Mountain Supergroup which is roughly stratigraphically equivalent to the Belt-Purcell Supergroup. However, recent isotopic dating of the mafic intrusions of the Wernecke Supergroup re-veal that it is at least 150 million years older than the Belt-Purcell indicating that it may actually be Early Proterozoic in age (Mortensen, 1998; personal communication). Sequence C, also known as the Windermere Supergroup, is the only sequence to be present along the entire length of the Cordillera. The Windermere Supergroup is characterized by the presence of glaciogenic gravel-sand-mud mixtures (diamictites) and mafic volcanic rocks near its base. It is generally agreed that Sequence C, which was deposited from 0.78 to 0.57 Ga, represents a rifted margin sequence that developed along Chapter 1. INTRODUCTION 5 the western margin of cratonic North America. Isotopic dating of the mafic intrusions at the base of the Windermere indicate two separate rifting events (780-730 M a and 555 Ma), each of which contributed to the opening of the proto-Pacific Ocean. The evidence for these extensional events is overprinted by shortening, folding and overthrusting of the supracrustal miogeoclinal sediments over the undeformed cratonic basement. These miogeoclinal sediments have been found to extend as far west as the Cordilleran Omineca Crystalline Belt (Monger and Price, 1979). To the west of the Cordilleran Deformation Front, the Canadian Cordillera can be divided into five morphogeological belts (Fig. 1.2) of distinct lithologies, structural styles and morphologies which were accreted to and deformed the western margin of ances-tral North America. From east to west, the five belts are the Foreland Belt, Omineca Belt, Intermontane Belt, Coast Belt and the westernmost Insular Belt. In the southern Cordillera outboard of the Insular Belt lies an accretionary wedge resulting from the eastward subduction of the Juan de Fuca plate beneath the North American plate off the west side of Vancouver Island. In the northern Cordillera, the plate boundary exists as the N-S trending transform Queen Charlotte-Fairweather Fault system separating the Pacific and North American plates. The Cordillera's morphogeological belts can be subdivided into two groups: those made up of deformed North American rocks and those made up of pericratonic or exotic rocks. The former group consists of the Foreland and Omineca belts while the latter group is made up of the Intermontane, Coast and Insular belts. This thesis focuses on the belts of North American origin that abut against the edge of cratonic North America, specifically the Foreland and Omineca belts. During the compressional events that resulted from the accretion of the western Cor-dilleran belts during the Middle Jurassic (ca. 170 Ma), the passive margin sequences along the rifted edge of cratonic North America were folded, detached from crystalline Chapter 1. INTRODUCTION 6 Figure 1.2: The morphogeological belts of the Canadian Cordillera (redrawn from Monger 1993). The Pacific, Explorer and Juan de Fuca tectonic plates are included in the diagram outboard of the Insular Belt. basement and upthrust onto the edge of ancestral North America to form the moun-tainous Foreland Belt. The rocks of the North American miogeocline that make up the Foreland Belt have been thrust over western Laurentia as far as 50 km in the northern Cordillera and 200 km in the southern Cordillera (Gabrielse et al., 1991). To the west, the Omineca Belt, which consists of the westernmost section of cratonic North America's miogeocline, was uplifted and metamorphosed during the compressional events that formed the Foreland Belt. The belt has been subjected to significant tecto-nism and deformation as evidenced by widespread folding and an imbricated succession of folded thrust sheets of Mesozoic age with eastward directed thrusting (Gabrielse et al., 1991). The Omineca Belt's western boundary corresponds to the 8 7 S r / 8 6 S r line. This line separates rocks with 8 7 S r / 8 6 S r ratios less than 0.706, which are associated with oceanic Chapter 1. INTRODUCTION 7 crust, from rocks with greater 8 7 S r / 8 6 S r ratios, which are associated with continental crust. 8 7 S r / 8 6 S r ratios act as indicators of basement composition since magma will take on the geochemical signature of the crust through which it travels. These isotopic results indicate that the Omineca Belt and all regions east of it in the Cordillera are primarily of North American origin. Thus, the belt is said to straddle the boundary between the accreted terranes to the west and the North American miogeocline to the east. During the Eocene (100 to 40 Ma), extension was the primary tectonic process af-fecting the eastern portion of the southern Canadian Cordillera. This extensional period resulted in transtension which followed the earlier formation of large scale right-lateral transcurrent faults with north-westward trends such as the Northern Rocky Mountain Trench in British Columbia and the Tintina Fault in the Yukon Territories. The Rocky Mountain Trench and the Tintina Fault together form a topographic lineament that extends 2600 km from Alaska to Montana (Gabrielse, 1991). While the Tintina Fault and the Northern Rocky Mountain Trench are characterized by major dextral transcurrent displacements, the Southern Rocky Mountain Trench which consists of three arc-like segments, is not associated with any transcurrent displacements and appears to have undergone significant Paleocene compression directed at right angles to it (Gabrielse, 1985). It is possible that the Southern and Northern Rocky Mountain Trenches once formed a continuous transcurrent fault that was disrupted during the Eocene extensional period. The Tintina Fault, which is a Late Cretaceous to Oligocene brittle structure exten-ding up to 2000 km in length, defines the boundary between the Foreland Belt and the Omineca Belt in the southern Yukon. The fault has been interpreted as a crustal-scale intracontinental transform fault due to its unusual continuity (Gabrielse, 1967) and its juxtaposition of crust with distinct physical properties (Lowe et al., 1994). For these Chapter 1. INTRODUCTION 8 same reasons, it is believed to be the northern extension of the Northern Rocky Moun-tain Trench. The Tintina Fault is responsible for 500 km (Tempelman-Kluit, 1977) to 750 km (Gabrielse, 1985) of right-lateral displacement of the cratonic edge of ancestral North America. 1.2.2 Regional Tectonic Setting for Line 21 While the morphogeological belts of the Canadian Cordillera provide a large scale frame-work for interpretation of the tectonic processes that have affected the Cordillera and the edge of cratonic North America, terrane analysis, introduced by Coney et al. (1980), al-lows the division of the Cordillera into distinct terranes for more detailed interpretation. A terrane is defined as "a body of rocks with an internally coherent stratigraphy and a tectonic history that differs from that of adjacent rocks in ways that cannot be easily explained in terms of facies changes" (Coney et al., 1980). The regional tectonic setting for Line 21 (Fig. 1.3) will be presented in terms of both terranes and morphogeological belts. Line 21 of the refraction experiment lies within Corridor 2 of the SNorCLE transect. It extends from the northwest portion of the Western Canada Sedimentary Basin (WCSB), across the deformation front of the northern Cordillera into deformed North American rocks and across the Tintina Fault, the eastern boundary of accreted terranes. With respect to current lines of longitude, the belts transected by Line 21 were all in place by Cretaceous time (Monger, 1993). The easternmost coverage of Line 21 extends onto the western margin of ancestral North America which experienced significant tectonic activity during the Proterozoic (see Fig. 1.4) evidenced by exposed basement rocks and seismic reflection results from Line 11 of the SNorCLE Transect (Cook et al., 1999). Chapter 1. INTRODUCTION 9 -130° -125° -120° Insular Superterrane E£H Alexander Coast Belt Undivided metamorphics and intrusives Accreted Terranes N Taku | H Stikinia I CC I Cache Creek I DY I Dorsey Quesnellia Slide Mountain Pericratonic Terranes I KO j Kootenay Displace Continental Margin I CA I Cassiar Miogeocline 3 Ancestral North America * , Morphogeological belts ""•t. Faults Deformation Front *fV_ SNoRE Line 21 Shotpoints and Seismographs (inline) ""N^  Seismographs along Line 22 (broadside) v ^ Political Boundaries -130° -125° -120° 100km 200 Figure 1.3: Terrane map of SNorCLE transect surrounding Line 21. Both Line 21 and the northern portion of Line 22 along which broadside receivers were deployed are shown. Shot locations along Line 21 are numbered and shown as stars. Morphogeological belt boundaries are superimposed on the map with thick dashed lines. Both the Tintina Fault and the Great Slave Lake Shear Zone (GSLSZ) are indicated as thinner dashed lines. Dashed line with west-pointing triangles delineates the Cordilleran deformation front. From 1.92 to 1.90 Ga, the western margin of the Archean Slave Craton existed as a rifted margin overlain by a passive margin sequence known as the Coronation Supergroup. This rifted margin was subsequently affected by the Wopmay Orogen (Fig. 1.4) which involved: (1) the Slave's collision with the Hottah Arc between 1.90-1.88 Ga due to westward directed subduction, (2) the development of the Great Bear Magmatic Arc between the Hottah and the Coronation Supergroup from 1.88 to 1.84 Ga due to eastward subduction and (3) the final collision of the Fort Simpson and Nahanni terranes at around 1.84 Ga. Recent seismic reflection results (Cook et al., 1999) have revealed that the Fort Simpson collision resulted in delamination of the Fort Simpson lower crust beneath the Chapter 1. INTRODUCTION 10 Figure 1.4: Diagram of Proterozoic Wopmay Orogen along the margin of cratonic North America based on description from Cook et al. (1999). Hottah Terrane. The easternmost extent of Line 21 begins in the Fort Simpson terrane and extends westward into the enigmatic Nahanni terrane, the westernmost margin of ancestral North America. While both the Nahanni and Fort Simpson terranes are overlain by significant thicknesses of sediments from the Western Canada Sedimentary Basin, their classification into two separate terranes is possible based on their distinctive aeromagnetic characteri-stics. Two basement drill cores collected for industry in the Fort Simpson terrane have revealed undeformed biotite granitic basement with an isotopic crystallization age of 1845 Chapter 1. INTRODUCTION 11 M a (Villeneuve et al., 1991). These isotopic results have linked the Precambrian basement of the Fort Simpson terrane with all the Cordilleran basement gneisses delineating a ribbon of Early Proterozoic rocks that may have controlled the rifting geometry that shaped the cratonic edge of North America (Ross, 1991). The Fort Simpson terrane is believed to have been accreted during the Wopmay orogen due to its age overlap with the end stages of magmatism in the Great Bear Magmatic Arc, 200 km to the east. Due to its high aeromagnetic signature, which may be due to uplifted highly magnetic igneous basement rocks (Cook et al., 1991; Cook et al., 1992), the Fort Simpson terrane is interpreted as a continental arc, that developed on the eastern edge of the Nahanni terrane. ''' Further west, the nature of the Nahanni terrane is not well understood as it lies beneath a considerable thickness of at least 10 km of passive margin sediments. To the north of Line 21, the Nahanni terrane is pinched out by the sharp eastward bulging of the Cordilleran deformation front. To the south of Line 21, the Nahanni terrane is truncated by the Great Slave Lake Shear Zone. As one moves west along Line 21 into the Cordillera, the first Cordilleran belt that is encountered is the mountainous Foreland Belt which forms the eastern mountain ranges and foothills of the Canadian Cordillera (Gabrielse et al., 1991). This classic fold and thrust belt is composed of folded miogeoclinal and clastic wedge assemblages of Mid-Proterozoic to Jurassic age that were deposited on and along the flank of the western rifted margin of ancestral North America. The eastern limit of the Foreland Belt corresponds with the Cordilleran deformation front which trends roughly north-west from south-central Alberta to the Yukon/BC bor-der where it sharply swings east at the Liard Line (explained in next subsection) and continues north to the Arctic coast. The belt's western limit in northern British Co-lumbia corresponds with the Northern Rocky Mountain Trench/Tintina Fault. In the Chapter 1. INTRODUCTION 12 Yukon, the western limit of the Foreland Belt diverges eastward from the Tintina Fault and follows the same trend as the Cordilleran deformation front. The westernmost ex-tent of Line 21 crosses the Tintina Fault into the Omineca Belt, a mountainous region of uplift underlain by metamorphic and granitic rocks. The belt is primarily composed of terranes of North American origin that have been displaced north-westward along the Northern Rocky Mountain Trench, its eastern limit. 1.2.3 Geology of Line 21 The surface geology along Line 21 (Fig. 1.5) is dominated by sedimentary rocks that mask the basement structure. The easternmost portion of the line samples the Phanerozoic cover of the Western Canada Sedimentary Basin, with Mesozoic sediments dominating the line east of the Cordilleran deformation front. To the west of the deformation front along Line 21, the surface sediments are predominantly Paleozoic with some basins of Mesoproterozoic sediments present in the foreland of the northern Rocky Mountains. Studies of the early Paleozoic tectonic framework of ancestral North America's mio-geocline characterize the easternmost portion of Line 21 as sampling the MacDonald Block which is sandwiched between the northeast trending Great Slave Lake Shear Zone to the south and the northeast trending Liard Line to the north (Fig. 1.6; Cecile et al., 1997). The Liard Line, interpreted as an ancestral transfer fault that coincides approxi-mately with the Yukon-BC border, represents a first-order structural discontinuity in the cratonic miogeocline. To the north, the Paleozoic miogeoclinal sediments are preserved and extend north-eastward 400 km onto the Canadian Shield beyond the Cordilleran de-formation front's eastward bulge; to its south, the time-equivalent miogeoclinal sediments are restricted to the deformed Cordilleran belts. The sedimentary discrepancy across the Liard Line is attributed to periods of Early Paleozoic erosion in the south due to elevated topography (Cecile et al., 1997). Chapter 1. INTRODUCTION 13 Geologic Assemblages • I I : I Cenozoic sedimentary rocks Cenozoic volcanic rocks Mesozoic-Cenozoic sedimentary rocks Mesozoic sedimentary and volcanic rocks Mesozoic sedimentary rocks Mesozoic volcanic rocks Paleozoic-Mesozoic sedimentary rocks Paleozoic-Mesozoic sedimentary and volcanic rocks Paleozoic-Mesozoic volcanic rocks Paleozoic sedimentary rocks Paleozoic metamorphic rocks Paleozoic sedimentary and volcanic rocks Paleozoic volcanic rocks Mesoproterozoic sedimentary rocks Mesoproterozoic-Neoproterozoic sedimentary rocks Neoproterozoic sedimentary rocks Neoproterozoic-Mesozoic sedimentary rocks Neoproterozoic-Paleozoic sedimentary rocks Precambrian metamorphic rocks metamorphic rocks of unknown age Plutonic and Ultramafic Rocks • Cenozoic Intrusive rocks • Mesozoic-Cenozoic Intrusive rocks • Mesozoic Intrusive rocks • Paleozoic Intrusive rocks Morphogeological Belts v> Faults Deformation Front SNORE Line 21 Shotpoints and Seismographs (line) Political Boundaries 100km 200 Figure 1.5: Geological map for the region surrounding SNoRE '97 Line 21. Line 21 is shown with its shot locations numbered and shown as stars. Morphogeological belt boundaries are superimposed on the map with thick dashed lines. Both the Tintina Fault and the Great Slave Lake Shear Zone (GSLSZ) are indicated as thinner dashed lines. Dashed line with west-pointing triangles delineates the Cordilleran deformation front. Chapter 1. INTRODUCTION 14 -130° -125° -120° Cratonic North America Figure 1.6: Structural paleogeographic feature map for Line 21 adapted from Cecile et al., 1997. Line 21 is shown with its shot locations numbered and shown as stars. Dashed line with west-pointing triangles delineates the Cordilleran deformation front. Drill well information is readily available for the easternmost 150 km of the line that transects the Western Canada Sedimentary Basin. Using information obtained from the Geological Atlas of the Western Canada Sedimentary Basin (Mossop and Shetsen, 1994), a cross-section for the upper 5 kilometers of the easternmost portion of Line 21 was constructed (Fig. 1.7). The geology shown along the surface of the cross-section was based on information included in Fig. 1.5. The prominent vertical offset in the Precambrian basement structure is an inferred feature based on the location of the Bovie Lake Fault Zone just west of shot location 2102. This northeast-trending normal fault is considered to have offset the basement by 2.5-5 km (Mossop and Shetsen, 1994). The cross-section was constructed using the minimum vertical offset. Sediment layer Chapter 1. INTRODUCTION 15 thicknesses below shot 2104 (i.e., the thicknesses west of the Cordilleran deformation front) provided in the Atlas (Mossop and Shetsen, 1994) had been determined using palinspastic reconstructions. Their positions were then assigned relative to a 2.5 km vertical offset in the Precambrian basement. 2104, 1000 -1000 -2000 -3000 -4000 -5000 Cordilleran Deformation Front Level .5000 Figure 1.7: Cross-section of surface sediments for easternmost portion of Line 21. Question marks indicate depths for which no information is available from the Atlas. Relevant shot points are indicated along the simplified topography. From the fold and thrust belt westward to the Tintina Fault, the region is overlain by drift cover sediments of the Liard Plain. The fault's surface expression along the Yukon/BC border is masked by these same sediments (Gabrielse, 1985). However, as the Tintina is believed to be the northern extension of the Northern Rocky Mountain Trench, it is most often mapped by extrapolating the Tintina Fault's known location Chapter 1. INTRODUCTION 16 southward and the Northern Rocky Mountain Trench's known location northward until they intersect. As such, the Tintina Fault's exact location at the Yukon/BC border is uncertain. 1.3 Previous work in northeastern British Columbia 1.3.1 Geological Mapping The Geological Survey of Canada has been collecting geological data from the Cana-dian Cordillera for a little over 120 years (reviewed by Monger, 1993). Their efforts were initially concentrated in the southern Cordillera but were soon expanded into the northern Cordillera such that by 1889, some data had been acquired for all regions of the Cordillera. Full regional geological coverage was completed by the early 1980s. The simplified geological map of Fig. 1.5 derives from the many geological studies. 1.3.2 Geophysics in the Northern Canadian Cordillera In contrast to geological surveying, geophysical studies of the Canadian Cordillera only began in the 1950s (reviewed by Berry et al., 1971). These have consisted of gravity, heat flow, magnetotelluric and geomagnetic depth sounding, and paleomagnetic, aero-magnetic and seismic surveys which have been concentrated primarily in the southern Cordillera with some studies near the Yukon/Alaska border and in the Northwest Territo-ries (Sweeney et al., 1991). Geophysical surveying of the region surrounding Line 21 also began in the 1950s; however the number of surveys have been few. Recent geophysical studies undertaken as part of the L I T H O P R O B E project are improving the geophysical database of the region. Chapter 1. INTRODUCTION 17 Seismic Surveys Several seismic surveys have been conducted in the northern Cordillera in areas sur-rounding the Line 21 area. The Peace River Arch Seismic Experiment (PRASE) was conducted to the southeast of Line 21 (Zelt and Ellis, 1989) although its contribution to the understanding of Line 21 is questionable due to the experiment's location south-east of the Great Slave Lake Shear Zone. Nonetheless, Line B of the experiment which samples an E-W line to the south of Line 21, resolves a westward thickening of low velocity sediments at the surface and a diffuse crust-mantle boundary. Two deep reflection lines recorded further north in the Northwest Territories (Cook et al., 1991), a reflection line along the Alaska Highway to the south of Line 21 (Cook and Van der Velden, 1993), and the reflection survey recorded along Line 11 of the SNorCLE transect (Cook et.al., 1999) across the western margin of ancestral North America have revealed a regionally extensive west-facing ramp at depth which has been confirmed by borehole data (Fig. 1.8). The ramp has been interpreted as the rifted Proterozoic margin of ancestral North America which was later overlain by passive margin sediments. This ramp has also been detected in seismic reflection surveys from the southern Cordillera and should exist at depth along Line 21 where it is expected to coincide with the gravity gradient as it does in the north (Cook et al., 1991). This ramp may have caused uplift of basement rocks as thrust faults rode over the ramped edge of cratonic North America (Cook et al., 1992). With the exception of one deep seismic study conducted by Tatel and Tuve in 1955 (Tatel and Tuve, 1956), the SNorCLE Corridor 2 region had never been studied through refraction or reflection seismology prior to the SNoRE experiment. Phase V of LITHO-P R O B E involves deep seismic reflection experiments along Lines 21, 22 and 31 of the Chapter 1. INTRODUCTION 18 Figure 1.8: West-facing ramp of cratonic North America's edge (adapted from Cook et al. 1991). Letters A-D and SNoRE Line 11 indicate reflection lines where ramp has been resolved. Location of SNoRE Line 21 is shown. SNorCLE transect (Fig. 1.1) during the winter of 1999-2000 (Clowes, 1997). These re-flection studies, coupled with the refraction data, will provide excellent seismic coverage of the region. Since 1981, the Geological Survey of Canada has operated one seismic station along Line 21 at Muncho Lake, which is located halfway between Fort Nelson, B C , and Watson Lake, Yukon. The station has provided some information about the North American miogeocline (Lowe et al., 1994). As its database of recordings grows, it will inevitably become more useful. Chapter 1. INTRODUCTION 19 Gravity Surveys Prior to 1991, only low resolution gravity measurements had been taken in the northern Cordillera near and around the area of Line 21. These surveys were confined to data collection every 1-2 km along the main highways. However, following a three year program of enhanced gravity mapping which started in 1991, the resolutions have greatly improved with station intervals on the order of 10 km over the region (Lowe et al., 1994). The new maps denote a fairly homogeneous gravitational high in areas east of the Tintina Fault/Northern Rocky Mountain Trench while the trench itself appears as a gravitational low. West of the Tintina, the gravity signature is more varied and contains anomalies with shorter wavelengths. These varying gravitational signatures are said to reflect the transition from the thick craton of ancestral North America in the east to the thin lower crust of the Cordillera in the west. It is important to note that the Tintina's strong gravitational low signal fades in the vicinity of the Yukon/BC border making the extrapolation of the Tintina down into the Northern Rocky Mountain Trench more speculative in terms of its exact location. Magnetic Surveys Between 1957 and 1990, a total of seven high resolution aeromagnetic surveys were con-ducted over the northern Cordillera (Lowe et al., 1994). Despite fragmentary coverage, the aeromagnetic surveys reveal a pattern of large long wavelength anomalies east of the Tintina Fault/Northern Rocky Mountain Trench which are particularly strong in the northeastern corner of British Columbia. These anomalies have been associated with the Rocky Mountains in the Foreland Belt and the Proterozoic craton to the east. West of the Tintina fault/Northern Rocky Mountain Trench, the aeromagnetic signal reveals anom-alies with shorter wavelengths which are believed to reflect shallower depths to magnetic Chapter 1. INTRODUCTION 20 basement (Berry et al., 1971; Lowe et al., 1994). The Tintina Fault itself shows up as a very long quasi-continuous magnetic feature, Magnetotelluric Surveys In the summer of 1996, magnetotelluric studies to determine subsurface electrical resis-tivity were conducted to the northeast of the study area in Corridor 1 of the S N O R C L E transect (Wu et al., 1998). M T data collection along Line 21 in Corridor 2 has been ap-proved for Phase V of the L I T H O P R O B E project (Clowes, 1997). Principal acquisition will take place in 1999-2000. Heat Flow Surveys Over the past few years, heat flow measurements have been collected in crystalline rocks of northern British Columbia and the southern Yukon. These measurements reveal an average measured heat flow for the entire Cordillera of 88 m W / m 2 and heat flow averages of 103.3 ± 2 1 . 0 m W / m 2 for the northern Cordillera (north of 59° latitude). Given that the heat flow measurements for Yellowknife, Ft. Providence and the Nahanni terrane are 54 m W / m 2 , 100 m W / m 2 and 118-136 m W / m 2 respectively (Lewis and Hyndman, 1998; Lewis and Hyndman, 1999), the S N O R C L E transect follows the trend observed in the southern Cordillera, that heat flow increases from ancestral North America into the Cordillera. The differing heat flow values suggest that the thermal transition from Cordillera to ancestral North American craton extends into the Wopmay Orogen at least to Fort Providence, N W T (Lewis and Hyndman, 1998). Such an interpretation of the observations contrasts greatly with those from the southern Cordillera where the thermal transition lies west of the Foreland Belt. However, the high heat flow values also could be due to high values of heat production in the crustal rocks, a suggestion which requires further study. Chapter 1. INTRODUCTION 21 1.4 Experimental Objectives The objectives for this study are three-fold: firstly, to determine the velocity structure of the crust and upper mantle from the Western Canada Sedimentary Basin, across the deformation front of the northern Cordillera and into deformed North American rocks, to the Tintina Fault, the eastern boundary of the accreted terranes; secondly, to resolve the transition at depth between the eastern portion of the Cordillera and the Proterozoic-Paleozoic rifted margin of ancestral North America against which it has been thrust; and thirdly, to integrate the seismic refraction results with other geophysical and geological information to provide greater insight into the tectonic evolution of northeastern British Columbia. Chapter 2 D A T A ACQUISITION A N D P R O C E S S I N G 2.1 The 1997 Slave-Northern Cordillera Refraction Experiment The L I T H O P R O B E Slave-Northern Cordillera Refraction Experiment (SNoRE '97) was conducted in northern British Columbia, the Yukon and the Northwest Territories of Canada during 4 weeks in August-September 1997 (Fig. 1.1). The experiment was de-signed to' determine the velocity structure of the crust and the upper mantle from Great Slave Lake in the Northwest Territories to the coast of the Alaska Panhandle. The main goals of the experiment were 1) to develop a near-continuous velocity structural cross-section from the Archean Slave province to the young accreted terranes of the western Cordillera, 2) to help resolve any-lateral variability of crustal structure along strike in the northern Cordillera, 3) to tie with the NSF-funded A C C R E T E reflection/refraction experiment in the southern Alaska Panhandle, 4) to compare the velocity structure of the southern and northern Canadian Cordillera, and 5) to compare seismic results with those determined for a similar region of Alaska. To achieve these goals, four refraction/wide angle reflection profiles were collected along the three corridors of the SNorCLE tran-sect. A total of 87 scientists from the Geological Survey of Canada (GSC), the United States Geological Survey (USGS), IRIS-PASSCAL, and the universities of British Co-lumbia, Victoria, Saskatchewan, Western Ontario, Manitoba, Alberta, Calgary, Stanford and Wisconsin (Eau Claire) took part in the experiment. 22 Chapter 2. DATA ACQUISITION AND PROCESSING 23 2.1.1 Location The SNoRE '97 refraction/wide-angle reflection survey comprised four refraction/wide-angle reflection profiles numbered 11, 21, 22 and 31 (Fig. 1.1). Line 11 in Corridor 1, which is the main focus of Fernandez Viejo et al. (1999a, 1999b), is a west-southwest line, 600-km-long, running from east of Yellowknife to west of Nahanni Butte in the Northwest Territories. This line coincides with a L I T H O P R O B E near-vertical seismic reflection survey (Cook et al., 1999) which crosses the Archean and Proterozoic components of the SNorCLE transect. Line 22 in Corridor 2, which is the main focus of Hammer et al. (1999), extends 550-km south from Watson Lake, Yukon, to Stewart, B C , across the accreted terranes of the Cordillera where it ties in with the A C C R E T E seismic survey (Morozov et al., 1999). Line 31 in Corridor 3, which is the main focus of Creaser and Spence (1999), is 300-km long and extends northeast to southwest across the Tintina Fault in central Yukon Territory from the miogeocline to the accreted terranes. This thesis focuses on the modeling and interpretation of Line 21 in Corridor 2, which is 460-km-long and extends along the Alaska Highway from Fort Nelson, B C , west-northwest to Watson Lake, Yukon, from the edge of cratonic North America to the western edge of the Cordilleran belts of North American origin. The line coincides with a L I T H O P R O B E near-vertical seismic reflection survey which will be collected in November and December, 1999. 2.2 Data Acquisition The SNoRE data were recorded during four deployments of recording instruments. Line 21 was the first, deployed on August 14 t h and 15 t h , 1997, and recorded on August 15 t h , 1997. A total of 596 portable seismographs were deployed inline along Line 21 (station spacing ~ 1.0 km) and broadside along the northern portion of Line 22 (station spacing Chapter 2. DATA ACQUISITION AND PROCESSING 24 ~ 2.0 km) to record arrivals from 8 shot points (see Table 2.1 for shot details). Shot Latitude Longitude Julian Shot Elevation Charge No. of No. Day Time (m) (kg) holes 2101 58.2000 122.7750 228 6:10:00 563 3000 3 2102 58.7917 123.5333 228 9:00:00 423 2400 2 2104 58.8362 125.0253 228 6:00:00 670 1200 1 2105 59.5761 126.5400 229 6:20:00 548 2000 2 2106 59.7000 127.2833 229 6:30:00 670 2000 2 2107 59.9750 128.2000 228 6:30:00 623 2000 2 2108 60.0750 128.9500 228 9:20:00 630 2400 2 2109 60.2000 129.5167 228 6:40:00 788 3000 3 Table 2.1: Shot point information for Line 21. 2.2.1 Drilling and Loading of Source Locations The drilling and loading of the shot locations took place well before the arrival of the deployment crew to the field. In most cases, conventional drill trucks were used to drill holes to a maximum depth of 60 m although the depths of individual holes were directly dependent on the amount of explosives to be used and on the diameter of the drill bits used (usually 21 cm). The upper parts of the holes were cased where necessary with steel piping to prevent collapse of the drill hole walls. Hydromex T3 explosives were loaded into the lower third of the holes and then the holes were back-filled. In all shot holes, the explosives were inter-layered with up to three booster charges distributed near the base, middle and top of the explosives. These boosters were attached to individual electrical detonating caps with wires that extended to the surface which allowed the shooters to detonate the explosions at a safe distance. For shot locations with greater than 1200 kg of explosives, up to three holes were drilled approximately 30 to 50 meters apart and the charge was distributed between Chapter 2. DATA ACQUISITION AND PROCESSING 25 them. The drilling of multiple closely spaced holes is advantageous for several reasons: • • it is cheaper and easier than deeper drilling, • for large experiments, the closely spaced holes approximate a point source, • multiple hole explosions minimize surficial damage, • multiple hole explosions are more likely to explode uniformly and completely. A total of 36 shots comprising 90000 kg of Hydromex T3 explosive were detonated during the entire experiment. Most shot sites used an average of 1000 kg to 3000 kg of explosives. For the eight inline shots on Line 21, a total of 18000 kg of explosives were used. 2.2.2 Surveying Shots and Receivers Surveying of shot and receiver locations for the SNoRE experiment took place several months before the refraction experiment and was performed using the Global Positioning System and 1:50 000 scale topographic maps. The Global Positioning System, which is funded and controlled by the United States Department of Defense, is a network of 24 satellites that continually orbit the Earth. The satellites are positioned so that at least four satellites can be seen by a receiver at any time and at any point on the Earth (including high latitudes). While the GPS system was developed as a United States military tool capable of precisions on the order of millimeters, worldwide non-US military use of GPS is available through the Standard Positioning Service (SPS) which provides coarse positioning with uncertainties of ~100 m horizontally and ~156 m vertically. Subsequent processing of the SPS field positions using GPS software can reduce these uncertainties to 5 m horizontally and 10 m vertically. Chapter 2. DATA ACQUISITION AND PROCESSING 26 A l l GPS shot and receiver locations determined in the field for the SNoRE experiment were later processed at U B C yielding uncertainties of ~10 m horizontally and ~15 m vertically. The surveying of shot sites for Line 21 of the experiment was mostly performed using 1:50 000 scale topographic maps because of necessary site location changes following earlier GPS site positioning. This procedure yielded uncertainties between 150 and 300 m both horizontally and vertically. Shots 2108 and 2109 were the only shot locations along Line 21 to be re-surveyed using GPS. To reduce the uncertainties associated with the shot points lacking GPS measurements, the uncertain shot locations were updated using local detailed maps, shot site descriptions and improvement of near shot traces. Initial surveying of receiver site locations was performed from May to July, 1997 using GPS. The initial receiver sites were then marked with a stake and annotated with their respective station numbers. A l l information about the site locations was then entered into road logs that were used by deployers to easily locate the sites. For significant site relocations, which were flagged in the road logs by the deployers, site re-surveying using GPS was performed after the experiment. Thirty-six sites for which precise locations were not obtainable due to insufficient satellite information avai-lable at time of measurement were nagged in the location file and the imprecise locations were adjusted based on arrival comparisons with adjacent receiver sites. Another 9 site locations lacked GPS data entirely and their locations were estimated using arrival time comparisons from adjacent receiver sites. 2.2.3 Instrument Deployment The deployment of the seismograph receivers was done primarily by teams of two de-ployers although some single deployers were used following an emergency situation along Chapter 2. DATA ACQUISITION AND PROCESSING 27 Line 21 (bear attack) to ensure that receiver coverage was complete prior to detonation of the shots. For each receiver deployment, deployers were required to follow a series of predetermined steps (taken directly from Gorman and Henstock, 1997): • locate the survey stake corresponding to the desired station number. • select a position within a few meters of the stake where the seismograph instrument box can be positioned free from natural and man-made hazards. • record identification numbers of all instruments left at each station in a field log. • dig a hole or holes to bury the geophone(s). A depth of at least 50 cm or a location on bedrock is optimal. • position and level geophones in the holes and attach them to the instrument box. • carefully bury the geophones to ensure good coupling with the ground. • tidy up the site to ensure that it will not be disturbed by animals or people during the period of data collection. Once deployment was complete, the shots were detonated during the night to reduce the amount of cultural noise. After all shots had been successfully detonated (which at times required a second night of shooting), the receivers were collected by the deployers and returned to the field headquarters for data downloading. 2.2.4 Recording Instruments The 596 portable seismographs used in this experiment (see Table 2.2 for instrument details) were of four types: 178 Scintrex-EDA model PRSls ; 20 Scintrex-EDA model PRS4s; 212 RefTeks and 186 SGRs. Both the PRSls and the SGRs are one-component Chapter 2. DATA ACQUISITION AND PROCESSING 28 seismographs while the RefTeks and PRS4s are three-component seismographs. The SGRs were predominantly used for broadside coverage throughout the experiment while the PRS and RefTek instruments were normally interspersed along the main lines to allow direct comparison between receiver timing and amplitudes from different instruments and to extend the lateral coverage of 3-component instruments. The PRS1 Portable Recording System instruments, provided by the Geological Survey of Canada (GSC), were specifically developed by the GSC to be used in large scale refraction surveys (Gorman and Henstock, 1997). These instruments were programmed at the field headquarters using the LithoSEIS software and then each was deployed with a geophone at the receiver sites. L-4 geophones (see Table 2.3 for details) were used with these instruments. The PRS1 instruments record only one channel (vertical component) and the data is stored in solid state memory which is uploaded once it is returned to the field headquarters. The data from all the PRS instruments is then compiled in the Society of Exploration Geophysicists standard SEG-Y format and archived onto exabyte tapes (Gorman and Henstock, 1997). The PRS4 instruments, also provided by the GSC, allow all the capabilities of the P R S l instruments to be applied to three components of data by recording on three L-4 geophones oriented in the vertical, east-west, and north-south directions. The RefTek three-component seismographs used in the experiment were developed by Refraction Technology Inc. and were provided by the Program for Array Seismic Studies of the Continental Lithosphere (PASSCAL) of the Incorporated Research Institutions for Seismology (IRIS). These instruments are programmed at the field headquarters using a palmtop computer and then each is deployed with a three-component geophone at the receiver sites. L-28 geophones (see Table 2.3 for details) were used with these instru-ments. The data are recorded by the RefTek instrument onto a disk and uploaded onto a Sun Sparcstation for compilation upon return to field headquarters. The data are then Chapter 2. DATA ACQUISITION AND PROCESSING 29 archived onto exabyte tape. The SGR Seismic Group Recorders, which were developed by Amoco and built by Globe Universal Sciences (Gorman and Henstock, 1997), were provided by the United States Geological Survey (USGS). L-4 and L-22 geophones (see Table 2.3 for details) were used with these instruments and the deployment was similar to the PRS1 instrument deployment. Data from the SGR instruments are digitally recorded onto special tape cassettes which are read into an Everex System 1800 microcomputer upon return to the field headquarters (Gorman and Henstock, 1997). These data are then compiled and archived on 9-track tapes. Instrument Type Number used Provided by Band Width Sampling Interval PRS1 178 GSC 0.1 to 100 Hz 8.3333 ms PRS4 20 GSC 0.1 to 100 Hz 8.3333 ms RefTek 212 IRIS-PASSCAL 1 to 125 Hz 8 ms SGR 186 USGS 2 to 200 Hz 2 ms Table 2.2: Seismic Instrument Information. Geophone Type Made by Band Width Sensitivity L-4 L-22 L-28 Mark Products Mark Products Mark Products 1 to 20 Hz 2 to 40 Hz 4.5 to 80 Hz 166.54 V/m/s 88 V/m/s 30.4 V/m/s Table 2.3: Geophone Information. In the sensitivity column, V/m/s signifies volts per meter per second. The inline receiver coverage along Line 21 involved the interspersed deployment of PRS1 instruments at odd-numbered stations (21-121 to 21-447) and RefTek instruments at even-numbered stations (21-148 to 21-586). PRS4 instruments were deployed at even-numbered stations (21-588 to 21-626). The SGR instruments were deployed both inline Chapter 2. DATA ACQ UISITION AND PROCESSING 30 (at odd-numbered stations from 21-449 to 21-625 and at all stations between 21-627 and 21-651) and broadside (at odd-numbered stations from 22-289 to 22-339). The broadside recordings which were designed to provide coverage across the Tintina Fault will not be considered in this thesis. 2.2.5 Timing of Shots and Receivers Since accurate timing records are essential for refraction seismology, shot timing was set by synchronization of the shot clocks with a satellite-linked GOES clock time signal. Drifts of the shot clocks were considered small (20-30 ms over three days) such that no significant corrections were applied. For both PRS and SGR instruments, receiver timing was set by rating the instrument clocks with the standard reference before and after deployment. The resulting timing errors, following corrections for clock drift, averaged between 30 and 40 ms. Receiver timing for RefTek instruments was performed in two ways. For RefTeks equipped with GPS receivers, the clocks were aligned with the GPS time signal and as such, experienced no time drift. The remaining RefTeks were calibrated by rating the instrument clock with a standard before and after deployment. The timing errors for these RefTeks averaged between 10 to 20 ms. Since the accuracy of the receiver clocks degrades with time deployed, on the order of tens of milliseconds per day, the receivers are returned to headquarters as soon as possible following data acquisition so that the clock drift can be determined and corrected for. For all instruments except the GPS-equipped RefTeks, a linear drift in the degradation of the timing was assumed and a linearly interpolated timing correction was applied during processing. It is important to note that the linear drift assumption may not hold and may result in significant errors for long deployments. Chapter 2. DATA ACQUISITION AND PROCESSING 31 The data record lengths recovered from both the PRS and the SGR instruments amounted to 60 seconds of recorded data for each shot. RefTek instruments, which have greater memory capabilities, recorded 160 seconds of data per shot. A l l record lengths were later merged and the PRS and SGR data records were padded to produce a single data set of 160 second record lengths. 2.3 Preparation of Data 2.3.1 Processing Routines The different data formats recorded by the four instrument types used in the experiment are all based on the Society of Exploration Geophysics standard known as SEG-Y. To enable different data formats from different instrument types to be combined into one standardized SEG-Y file format for each shot, the PLOTSEC (Amor, 1996) processing routines were developed. This package of programs allows the user to read, modify and merge different S E G - Y data formats together for each shot. For the initial step in the merging of the different data types, the data were read into the system using the plot sec jrsegy routine which was also required to resample the input PRS data at a sampling interval of 8 ms, using a cubic spline algorithm, to match the other data. 2.3.2 Gain Since different instrument-geophone combinations yield different amplitudes, the gains listed in Table 2.4 were applied to the data types as they were read in. These applied gains allowed an approximation of true-relative amplitudes between the data from different instrument-geophone combinations. The gains listed in Table 2.4 were assigned based on experimental values provided by the geophone and instrument makers and on visual Chapter 2. DATA ACQUISITION AND PROCESSING 32 examination of the amplitudes of traces from neighboring instrument types using the PRS Is as the standard. Instrument Type Geophone Type Gain Applied PRS1 L-4 1 PRS4 L-4 0.5 RefTek L-28 -3.289 x 107 SGR L-4 0.584 x 107 SGR L-22 -1.136 x 107 Table 2.4: Gain applied to each instrument-geophone combination. 2.3.3 Updating of S E G - Y Headers Once the data from all the instrument types were read into the system, the S E G - Y header files were updated using the plotsecupdate command. In most cases, the required updates involved the inclusion of shot information such as time of shot, shot size and shot and receiver survey information so that shot receiver distances could be calculated. (This was done using a SAC subroutine based upon the reference spheroid of 1968 and defined by a major radius of 6378.160 km and a flattening of 0.00335293 (Gorman and Henstock, 1997)). During this stage, all known shot and receiver timing corrections were applied to the data. 2.3.4 Merging Data and Creating Final Data Sets With the update of headers complete, the data from the different instrument types were merged using the plotsec_merge routine into a single standardized file for each shot point. These files were then translated from PLOTSEC format to S E G - Y format using the plotsec_wsegy routine. Chapter 2. DATA ACQUISITION AND PROCESSING 33 2.4 Data Processing In addition to the routines available for merging of different S E G - Y data formats, the PLOTSEC package contains routines to aid in the manipulation of the data once it has been properly merged together. These routines allow the user to merge different data sets according to time and/or offset, update parameters in the S E G - Y header file, stack data, interactively pick phase arrivals, apply different reducing velocities and plot data sections to the screen or on hardcopy according to time and/or offset. Using the plotsecJilt and the plotsec_plot routines, both true relative amplitude and trace normalized filtered data plots (shown in Appendix A) were constructed for all shot points. The trace normalized plots (Figs. A . l a , A.4a,..., A.22a) were normalized to the maximum amplitude along a given trace. 2.4.1 True Relative Amplitude Scaling To account for the decrease in amplitude with increasing offset caused by spherical sprea-ding and attenuation, the traces in the true relative amplitude plots (Figs. A.2,A.5, . . . , A.23) were scaled and individually multiplied by their source-receiver distance raised to a power of 1.85. The choice of 1.85 was determined based on visual improvement of the far-offset arrivals in the true relative amplitude plots. 2.4.2 Filtering The plotted data in Appendix A were minimum phase bandpass filtered from 0.5 to 10 Hz for all shots; particularly noisy traces were subsequently ehminated manually using plotsecJilt. During the picking stage (discussed in the next chapter), the filter was interactively applied and removed to aid in determining the first-breaks of the phase arrivals. Chapter 3 D A T A D E S C R I P T I O N A N D ANALYSIS 3.1 Data Quality The data set from Line 21 is of good quality, although it is the noisiest of the data sets from the SNoRE experiment. This may be due in part to its deployment along the Alaska Highway. Most shots, with the notable exception of shot 2104, have good signal-to-noise (S/N) ratios and primary phases (Pg, PmP and Pn) are observable for all shots except 2104 which does not show any clear Pn. The data quality is particularly good at the ends of the line with shots 2102 and 2108 having the highest S /N ratios. Shear wave arrivals are not distinguishable for most shots, with the exception again of 2104 where the Sg phase can be seen. Due to the sparsity of shear wave data, only the P-wave arrivals will be analyzed in this thesis. Data plots and related results are presented in Appendix A . Figures A.2, A.5, . . . , A.23 i display the data in true relative amplitude format. Figures A.1, A.4,.. . , A.22 show trace normalized plots with the picked phases identified. 3.1.1 Observations of Pg Phase The Pg phase, which represents rays turning in the upper crust, consists of the first arrivals for receivers at near offsets to the shots. The arrival times and amplitudes from this phase provide insight into the velocity and velocity gradient distributions within the upper crust. Visual comparison of Pg arrivals from all shots provides an estimate of 34 Chapter 3. DATA DESCRIPTION AND ANALYSIS 35 lateral variability of velocities and velocity gradients within the upper crust along the refraction line. The Pg phase is visible for all shots and is easily picked due to its relatively high S /N ratio. On all sections, Pg exists as the first arrival for offsets out to ~ 150 km. At near offsets of less than 15 km, the Pg arrivals delineate a surface layer of laterally variable velocities ranging from 2.0 km/s to 5.5 km/s. At offsets greater than 15 km, Pg arrivals follow linear branches indicating that the upper crust in the region does not appear to contain layers of sharply differing velocities or strong velocity gradients. The Pg arrival branches have an average apparent velocity of 6.1 km/s although the apparent velocity for each shot varies significantly along the line, indicating that the study area has strong lateral heterogeneity. Short wavelength variations in arrival times along the Pg branches are believed for the most part to be caused by near-surface features. However, there exists one important exception to this interpretation. Receivers located between shots 2109 and 2108 on the data plot from shot 2108 (Fig. 3.1) demonstrate Pg arrivals with significantly later arrivals over an offset range of 30 km. The location of the delayed arrivals along the refraction line, coincides with the approximate location of the Tintina Fault scar. The feature causing the delay of arrivals at this location appears to also be responsible for a delay in Pg arrivals from shots 2107 and 2109, as well as PmP arrivals from shots 2105, 2106 and 2107 and Pn arrivals from shots 2101 and 2102. The noisiest data set along the line is from shot 2104 (Fig. A.8) at the Cordilleran deformation front where Pg can only be seen for near offsets of up to 100 km. These arrivals show a significant increase in velocity to 6.4 km/s within 5 km of the surface to the northwest. While these velocities are quite high for upper crustal rocks, they are corroborated by fast apparent velocities recorded to the southeast from shot 2105. In terms of amplitudes, the Pg branches do not show prominent amplitude changes, Chapter 3. DATA DESCRIPTION AND ANALYSIS 36 Figure 3.1: Trace normalized plot for shot 2108 showing slowed Pg phase over a limited number of stations. Inset is an enlargement of data within the box outlined by the thin line shown. Thin horizontal bars on each trace indicate the picked traveltimes. indicating that the upper crust does not contain strong velocity gradients. The highest Pg amplitudes occur for shots at the ends of the line with the exception of shot 2109, at the westernmost end on the line, which has very low amplitudes at near offsets up to 50 km to the east (Fig. 3.2a). Pg amplitudes appear muted for all shots between the locations of shots 2102 and 2106 (Fig. 3.2b- 3.2d) which may be linked to structures in the Cordilleran deformation front. Chapter 3. DATA DESCRIPTION AND ANALYSIS 37 Muted Pg amplitudes?' m\fx~m\t^ M l Distance (km) 0 Offset (km) 0 100 150 200 250 300 350 400 450 Distance (km) 0 50 100 150 200 250 300 350 400 450 100 150 200 250 300 350 400 450 offset (km)-200 -150 -100 - 50 0 50 100 1 50 200 250 D f f iS" , f f ° ™ 5 0 ->nn , 0 0 , s n 1 5 ° , n n 2 0 0 , ™ 2 5 0 , o n 3 0 0 4 | 3 5 0 n 4 0 0 4 5 0 0 Distance (km) 0 50 100 150 200 250 300 350 400 450 unset(Km) - J50 -JOU -dx - i l l" -1MJ - T O ou u su Offset (km)-150 100 -50 0 50 100 1 50 200 250 300 Figure 3.2: Example data plots where Pg amplitudes appear muted. Muted regions are indic-ated with black ovals, a) Shot 2109; b) Shot 2105; c) Shot 2102; d) Shot 2106. 3.1.2 Observations of Crustal Reflectors Wide-angle crustal reflections appear as secondary arrivals after the Pg phase. When observed consistently over at least a few shots, they can provide information about the geometry of a subsurface boundary and the velocity contrast across it. When observed sporadically with inconsistent traveltimes, they are less informative but still must be considered. The data set from Line 21 of the SNoRE experiment exhibits few prominent crustal reflections. One prominent exception is reflector R l observed in the easternmost part of the line (Fig. A.4). While R l was initially picked as PmP due to its high amplitudes, Chapter 3. DATA DESCRIPTION AND ANALYSIS 38 modeling results to be presented in the next chapter reveal that it is in fact a lower crustal reflector that masks the later PmP arrivals in its coda. The lack of crustal reflectors indicates that the study region does not contain layers with sharp velocity discontinuities except at the base of the crust in the easternmost part of the line. A few short reflectors (R) were also picked from the data set (Fig. 3.1). These reflectors whose amplitudes were all relatively poor, are seen for a few shots but are not connected and as such are interpreted as being unassociated with velocity discontinuities at layer boundaries. 3.1.3 Observations of PmP Phase PmP arrivals represent rays reflecting from the crust-mantle boundary or Moho. They provide a strong constraint on the geometry of the boundary and a moderate constraint on the velocity contrast across it. PmP reflections are evident for all shots along Line 21 except shots 2101 and 2102 where the PmP arrivals are masked by the R l arrivals. A n example is shown in Figure 3.3 (see appendix A for the complete data set). While shots northwest of 2104 show clear PmP arrivals with strong amplitudes and distinct codas, shots 2104 and those to the southeast exhibit more complicated arrivals. Initially, the PmP phase appeared to be quite clear for all shots. However, upon further investigation using traveltime and amplitude modeling (to be discussed later) it was found that the strong PmP-type reflections for shots 2101, 2102 and the eastern arrivals from 2104 may have resulted from a strong lower crustal reflector (R l ) located a few kilometers above the Moho which served to mask the later PmP reflections in its coda. Given these results, the PmP arrivals to the west of shot point location 2104 provide good constraints on the depth of the Moho and on the velocities in the lower crust while to the southeast of location 2104, lower crustal velocities are constrained by the R l phase Chapter 3. DATA DESCRIPTION AND ANALYSIS 39 Figure 3.3: Trace normalized plot for shot 2107 showing examples of the PmP phase at both ends of the line. Inset is an enlargement of data within the box outlined by the thin line shown. Thin horizontal bars on each trace indicate the picked traveltimes. and Moho depth is constrained by the Pn phase. 3.1.4 Observations of P n Phase The Pn phase represents rays turning in the upper mantle and provides good constraints on upper mantle velocities and Moho depths. It appears as the first arrival for offsets greater than 150 km for all shots except 2104. While the Pn branches are not perfectly linear due to crustal heterogeneity (see Fig. 3.4 for an example), they do show an average Chapter 3. DATA DESCRIPTION AND ANALYSIS 40 apparent velocity of ~ 8.1 km/s for the upper mantle. On a shot by shot basis, however, the Pn branches show highly variable apparent velocities with velocities greater than ~ 8.2 km/s from the end of the line shots to less than ~ 8.0 km/s for more central shots. These apparent velocities indicate that the upper mantle is also strongly laterally heterogeneous. Figure 3.4: Trace normalized plot for shot 2102 showing an example of the Pn phase. Inset is an enlargement of data within the box outlined by the thin line shown. Thin horizontal bars on each trace indicate the picked traveltimes. The amplitudes of the Pn phase vary along the line from shot to shot. The end shots (shots 2101, 2102 and 2108) generally exhibit strong amplitudes at offsets greater than 300 km while the amplitudes decrease slightly at nearer offsets. Shot 2108 shows a highly Chapter 3. DATA DESCRIPTION AND ANALYSIS 41 variable Pn amplitude trend with amplitudes appearing muted between offsets of 200 km and 300 km (Fig. 3.5). The more central shots (shots 2105, 2106 and 2107) demonstrate lower Pn amplitudes overall, which may be due to the smaller shot sizes detonated at the central shot locations. Figure 3.5: True relative amplitude plot for shot 2108. Note the decrease in Pn amplitudes between offsets of 200 km and 300 km. 3.1.5 Summary of phase observations The good quality data from Line 21 had a high enough S /N ratio such that the main phases (Pg, PmP and Pn) were easily identifiable for most shots. The offsets to which all Chapter 3. DATA DESCRIPTION AND ANALYSIS 42 of the observed phases are identified for each shot are listed in Table 3.1. The cor-responding data plots with the phases identified are presented in Appendix A from Fig. A . l a , A.4a,..., A.22a. 3.2 Data Interpretation and Picking To enable the data set collected along Line 21 to be modeled into a velocity model that can reproduce the observed arrivals, the first arrivals of each of the phases described in the last few sections were located and picked from the data. The picking process is tremendously important and must be undertaken with care to reduce picking uncertainties since the traveltimes calculated for the constructed model will be directly compared to these picked traveltimes. Problems can arise during picking from a number of sources. Firstly, the individual phases must be located on the data plot and their first arrivals or first breaks must be picked. While this process can be quite straightforward for the Pg and Pn arrivals which are the first arrivals for near and far offsets respectively, the arrival of later reflections from the Moho and from crustal reflectors can be difficult to pinpoint since their first arrivals are superposed with the codas of the Pg arrivals. During picking, it is also important that the picks for a given phase be chosen at the same point along each trace such that if the picks are aligned, the peaks and troughs of the first few cycles will fine up. Once the phases have been properly correlated, uncertainties can still arise due to skipping of cycles with negligible amplitudes. Skipping of these cycles causes the arrivals to falsely appear delayed. With all the uncertainties involved in the picking process and the resolved model's overall dependence on the observed picks, pick determination is perhaps the most crucial step in modeling refraction data and as such, the picks are reevaluated iteratively as Chapter 3. DATA DESCRIPTION AND ANALYSIS 43 the modeling proceeds. To enable the data to be interpreted and manipulated (filtered, scaled, traces muted, traces killed, plotted, etc.), the P L O T S E C routines developed at U B C (Amor, 1996) were employed. The specific routines used at different stages of the analysis are specified when appropriate. Once the data sets for each instrument type had been merged together using proce-dures described in Chapter 2, the first arrivals for the primary phases were picked using large paper plots of varying scale and the interactive plotsec_pick routine. For each pick, the pick uncertainty was calculated using the plotsecamppk routine which assigns an uncertainty value based on the signal to noise ratio over a time window of 0.1 s before and after the pick. The SNoRE 97 experiment along Line 21 resulted in a total of 3908 traces. From these traces, 1360 traces provided traveltime picks for rays turning in the crust (Pg phase), 660 provided picks for rays reflecting off the Moho (PmP phase) and 1023 provided picks for rays turning in the upper mantle (Pn phase). Crustal reflectors are mostly absent from the data set with two notable exceptions. Firstly, a strong lower crustal reflector (Rl ) masked the PmP arrivals in the eastern portion of the line and provided 307 picks. Secondly, a few disconnected small scale reflectors (R) provided 361 picks. The number of picks for each shot and each phase and their corresponding average uncertainties are displayed in Table 3.2. Chapter 3. DATA DESCRIPTION AND ANALYSIS 44 Phase Shot Point Pg 2101 2102 2104 2105 2106 2107 2108 2109 PmP 2104 2105 2106 2107 2108 2109 Pn 2101 2102 2105 2106 2107 2108 2109 R 2102 2104 2105 2106 2107 2108 2109 R l 2101 2102 2104 N W 160 - 40 135 - 15 70 - 0 120 - 0 150 - 0 100-0 55 - 5 20 - 10 140 - 80 200 - 115 150 - 105 100 - 75 470 - 160 390 - 155 SE 5 - 60 0 - 130 0-100 0 - 135 0 - 130 0 - 145 5 - 170 85 - 155 85 - 165 70 - 170 95 - 175 80 - 170 150 - 75 50 - 30 75 - 15 75 - 30 225 - 80 175 - 90 195 - 255 190 - 300 150 - 360 190 - 405 195 - 435 45 - 130 35 - 55 60 - 140 90 - 150 Table 3.1: The source-receiver offsets in kilometers at which each phase was visible for a given shot point! NW indicates offsets to the northwest of the specified shot while SE indicates offsets to the southeast of the specified shot. Chapter 3. DATA DESCRIPTION AND ANALYSIS 45 Shot No. P g R R l PmP Pn Total 2101 91 0 180 0 250 521 (119) (0) (116) (0) (142) (129) 2102 150 97 109 0 204 560 (76) (105) (120) (0) (132) (82) 2104 127 31 18 67 0 243 (104) (106) (151) (135) (0) (117) 2105 195 47 0 151 18 411 (84) (141) (0) (148) (142) (116) 2106 263 51 0 123 68 505 (75) (110) (0) (133) (136) (98) 2107 210 91 0 128 156 585 (65) (133) (0) (111) (148) (181) 2108 182 16 0 95 166 459 (43) (30) (0) (104) (125) (61) 2109 142 28 0 96 161 427 (52) (132) (0) (116) (135) (69) Total 1360 361 307 660 1023 3711 (81) (103) (120) (130) (137) (112) Table 3.2: Number of observations (from traveltime picks) for each shot point with correspond-ing average traveltime uncertainties (in milliseconds) in brackets. The bottom row gives the total number of observations (with corresponding traveltime uncertainties) for each phase. Chapter 3. DATA DESCRIPTION AND ANALYSIS 46 3.3 Modeling Algorithms and Techniques Since the 1970s, analysis and modeling of refraction data sets has involved the use of forward two-dimensional ray tracing algorithms such as those developed by Cerveny et al. (1977), Whittall and Clowes (1979), McMechan and Mooney (1980), Spence et al. (1984), and Luetgert (1992). These early 2-D algorithms which performed ray tracing between sources and receivers with a trial and error approach, were limiting due to their unquantifiable non-uniqueness and their high demand on human interaction, especially for multi-shot seismic experiments with many receivers. Since the early 1980s, forward modeling ray tracing algorithms which incorporate inversion techniques have been deve-loped by Firbas (1981), Spence et al. (1985), Huang et al. (1986), White (1989), White and Clowes (1990), Lutter et al. (1990), Lutter and Nowack (1990), and Zelt and Smith (1992). These forward modeling/in version algorithms perform more efficient ray tracing, allow more effective updating of velocity models and reduce the non-uniqueness of the resolved velocity models. The algorithms also reduce the amount of time required to develop a suitable velocity model by combining user controlled forward modeling and computer controlled (inversion) procedures to establish the most appropriate model. The advantages of using inversion algorithms are that the data can be fit according to a specific norm, the model parameter resolution can be estimated and the seismologist obtains an appreciation of the uncertainty and the non-uniqueness of the problem (Zelt and Smith, 1992). Prior to the formulation of R A Y I N V R by Zelt and Smith (1992), inversion algorithms for refraction data used travel time calculations from forward mode-ling algorithms that were not designed with the inversion approach in mind. R A Y I N V R was chosen as the modeling routine for this study because it is an inversion algorithm that uses a forward modeling procedure designed according to the specific needs of the inversion. These needs are a flexible model parameterization, a minimum number of Chapter 3. DATA DESCRIPTION AND ANALYSIS 47 independent parameters and an efficient robust method of ray tracing (Zelt and Smith, 1992). Once the phase arrivals had been picked for all shot points, R A Y I N V R forward and 2-D inversion modeling were done to construct a 2-D model that would reproduce the observed traveltime arrivals. To account for the crooked shot-receiver arrays of Line 21, the data were projected onto a straight line connecting shot points 2101 and 2109 whilst maintaining their correct offset distances. Both the shot point locations and the elevation along Line 21 were projected perpendicularly onto the straight line. Projection of the data onto a straight line allowed the arrivals to be modeled in a 2-D plane. While Zelt and Zelt (1998) have determined that in most cases, out of plane effects contribute only relatively small errors to 2-D velocity models, significant 3-D geometrical effects can sometimes lead to conflicting traveltime arrivals from different shots that cannot be fully reconciled with a 2-D model. Nonetheless, due to Line 21's overall linearity, the available analysis procedures, and the lack of data to consider any variations outside of the vertical plane, interpretation as a 2-D model was considered the only viable approach. 3.3.1 Tau-p Analysis To obtain an initial 2-D model from the data to be used as an input model for R A Y I N V R , wave-field r — p analysis was used (Gorman and Clowes, 1999). Wave-field r — p analysis involves transforming the data from t — x (traveltime-distance) space to T — p space where T is the time intercept for a specific slowness and p is the slowness. This transformation was carried out according to the equation /oo P(T+ px,x)dx. (3.1) •CO P{r+px,x) is the observed seismogram wavefield and S(r,p) is the transformed wavefield. Chapter 3. DATA DESCRIPTION AND ANALYSIS 48 This transformation treats all the refraction arrivals in t—x space as plane wave elements with slope p and time-intercept r . Once the initial transformation has been executed, the T — p data are downward continued to z — p space using a user specified 1-D velocity function v(z) that is applied on a trace by trace basis and that can be iteratively amended by the user until a suitable fit is achieved. To ascertain the appropriateness of the v(z) function chosen for a given shot, 1-D ray-theoretical forward modeled time-distance curves were calculated using the W K B J technique (Chapman, 1978). These time-distance curves were then overlain on the ori-ginal data set to check their fit. While a perfect fit was not expected since a 1-D model cannot reproduce laterally varying 2-D data, an approximate fit provides a 1-D velocity model that can approximately reproduce the data from a specific shot. Using the r — p technique, 1-D velocity profiles were constructed for each shot. From these velocity profiles, which represent changes of velocity with depth (v(z)), 2-D velocity functions (v(x,z)) were constructed by incorporating the offset of the refracted ray as it turns within the 1-D velocity field. The 2-D profiles were then superimposed on a grid and contoured using G M T contouring routines (Wessel and Smith, 1995) to construct a quasi-2d velocity model to be used as an initial model in R A Y I N V R (Fig. 3.6). The great advantages to this technique are that the results are not based on any kind of picking routine and that a rough initial quasi-2d model can be constructed in a short amount of time. The main difficulty associated with the use of this procedure is that the individual phases lose their clarity when represented in r — p space making the choice of the 1-D velocity function extremely non-unique. Chapter 3. DATA DESCRIPTION AND ANALYSIS 49 WNW ESE 50 100 150 200 250 300 350 400 450 Distance (km) (km/s) 1.9 3.0 4.0 5.8 6.0 6.2 6.4 6.8 7.0 7.5 8.1 8.4 8.9 Figure 3.6: Initial r — p Velocity Model for Line 21 plotted with a vertical exaggeration of 5 to 2. The shot numbers and locations are plotted along the top of the 2-D velocity model. Dotted black lines in the model represent the regions where contour control exists from the 2-D velocity profiles. 3.3.2 R A Y I N V R Modeling The Forward Problem The forward problem in refraction seismology inversion is, simply stated, to determine the amount of time it takes to travel along a ray in a given velocity field between a source and a receiver. The time it takes to travel a given distance in a given constant velocity field is simply the distance traveled divided by that constant velocity. For a bending ray, the travel time along a ray can be determined by integrating the slowness (1/velocity) field along a ray according to the equation: Chapter 3. DATA DESCRIPTION AND ANALYSIS 50 t = JL v(x, z) (3.2) where L is the ray path. This equation can be discretized for practical applications such To solve the forward problem, it is essential to construct an initial velocity model and to calculate ray paths within the velocity model. Once the rays have been calculated, their travel times can be determined by integrating over the ray path using the trapezoidal rule. Model Parameterization The model parameterization for R A Y I N V R discretizes the velocity model in terms of both velocity and boundary nodes which divide up the model into trapezoids (Fig. 3.7). By assigning different velocity values to each of the velocity nodes defining the trapezoid, velocity gradients can be incorporated into each model trapezoid. The use of gradients allows the individual trapezoids to be quite large and simplifies the parameterization since fewer trapezoids are required to model a laterally complex region. Earlier ray tracing algorithms used constant velocity blocks which simplified the calculations but increased the number of blocks required to model a complicated velocity model. The velocity model parameterization consists of a two-dimensional layered model discretized into blocks of trapezoidal shape where the medium is assumed to be isotropic and laterally homogeneous in a direction normal to the plane of the model (Zelt and Ellis, 1988). The layers of the model are required to extend from one side of the model to the other without crossing any other layers. As such, the layers can be parameterized that: (3.3) Chapter 3. DATA DESCRIPTION AND ANALYSIS 51 Distance -V V Figure 3.7: Schematic diagram of RAYINVR parameterization. Solid dots represent bound-ary nodes and inverted triangles represent velocity nodes. Solid lines define horizontal layer boundaries while dashed lines represent vertical boundaries assigned by RAYINVR. to model pinch-outs and isolated bodies by simply assigning the layer's top and bottom boundaries to exist at the same depth. The trapezoidal blocks of the parameterization can have sloped boundaries at the top and bottom of the layer which can be used to model the topography of any layer. Boundary nodes are plotted in the model at points where the layer topography is specified. For layers with flat topography, only one boundary node is required to specify the depth of the layer. In addition to boundary nodes, all boundaries of the blocks in a layer can be assigned velocity gradients or can be of a fixed velocity depending on the needs of the user. These are plotted in the form of velocity points which appear in the model wherever the velocity is specified. The velocities along boundaries between velocity points are calculated using linear interpolation. Velocity gradients across boundaries can be discontinuous if the velocity defined for the bottom of a layer is different from the velocity at the top of the underlying layer. The velocity can also be made continuous across a boundary by simply Chapter 3. DATA DESCRIPTION AND ANALYSIS 52 assigning the same velocity value both to the bottom of the layer and to the top of the underlying one. The boundaries along the sides of the blocks are required to be vertical and are auto-matically assigned by the R A Y I N V R program. These vertical boundaries are assigned whenever a velocity or a boundary node is assigned along the top of a given model layer. The velocities assigned are P-wave velocities from which the S-wave velocities can be calculated directly using assigned Poisson ratios for each trapezoidal block. Within a trapezoidal block with four boundaries, the P-wave velocity is interpolated at all points as explained by Zelt and Ellis (1988). The complex interpolation provides a velocity field that varies linearly with position along all four boundaries of the block and ensures that all points in the velocity model are defined. Ray Tracing Algorithm Rays are traced through the model using zero-order asymptotic ray theory by numerically solving two sets of first-order differential equations (Cerveny et al., 1977). This system of equations is solved using the Runge-Kutta method (Sheriff and Geldart, 1983) and is subsequently verified to ensure that Snell's law is satisfied at all points where the rays encounter a model boundary. Blocky model parameterizations with velocity gradients and velocity discontinuities have the disadvantage of causing scattering and focusing of ray paths. To reduce this instability, R A Y I N V R performs a boundary smoothing simulation during ray tracing. The smoothing simulation has little effect on the ray tracing except to alter the take-off angles at the smoothed boundary. This operation is advantageous as it ensures that more rays will reach the surface to provide more traveltime data. The ray tracing algorithm within R A Y I N V R employs a variable step length, A , during Chapter 3. DATA DESCRIPTION AND ANALYSIS 53 ray tracing to ensure that the ray is appropriately sampled. The step length is defined as (3.4) where a is specified by the user, vx is the partial derivative of v in the x direction and vz is the partial derivative of v in the z direction. If a is set to a value between 0.025 and 0.1, total traveltime errors for ray paths can be as low as ±0.002 — 0.01.S (Zelt and Ellis, 1988). A varying step length ensures that bending rays in velocity fields with gradients are sampled more often than straight rays traveling through trapezoids of constant velocity. The Inverse Problem The inverse problem in seismic refraction is to determine a velocity model that can reproduce a given data set of traveltimes. Since the ray paths used to calculate the traveltimes are directly dependent on the velocity that they pass through, the problem is said to be non-linear. This inherent non-linearity in the seismic problem requires that the problem be linearized by considering perturbations in the velocity model rather than velocities themselves. This approach necessitates the use of an initial model and an iterative approach to deal with successive perturbations. The mathematics involved in this inverse approach are elaborated on in Appendix B. The Iterative Approach Nonlinear relationships require the implementation of an initial model and an iterative approach. Firstly, R A Y I N V R traces the specified rays through the initial model and then calculates the traveltimes for each ray. These traveltimes are subsequently compared with the observed traveltimes to calculate traveltime perturbations. The perturbations Chapter 3. DATA DESCRIPTION AND ANALYSIS 54 in the velocity model caused by the travel time perturbations are then computed at eachnteration step using a damped-least squares approach. These model perturbations are subsequently added to the velocity model which can then be compared against the stopping criterion. The Stopping Criterion The stopping criterion for an iterative inversion approach is satisfied when the best model is believed to have been produced. Typically, this best model is that which provides the desired trade-off between the travel time residuals (c/>A and the model resolution (c/>m) and which allows rays to be traced to all observation points (Zelt and Smith, 1992). However, the stopping criterion for a given nonlinear inversion problem is specific to the particular problem and to the needs of the individual seismologist. The damped least-squares technique, used in the inversion step of R A Y I N V R , amends the initial velocity model by adding in the perturbations it has calculated. R A Y I N V R is then run in the amended velocity model to trace rays and to calculate the R M S (root-mean-square) residual traveltimes and the %2 factor. If these have not reached the desired fit to the observed traveltimes (i.e. have not attained the stopping criterion), the damped least-squares approach can be reapplied to the traveltimes of the amended velocity model. This process is rerun iteratively until the desired fit is attained. Commonly, the stopping criterion are three-fold, 1) to minimize the RMS residual traveltime, 2) to minimize the \ 2 factor, and 3) to trace rays to all observation points. The R M S residual traveltime gives an average uncertainty for the model by comparing the observed and calculated traveltimes for each arrival. So as not to over-fit the obser-vations, an acceptable RMS residual traveltime should be on the order of the uncertainty of the picks. Meanwhile, the statistical %2 approach evaluates the differences between Chapter 3. DATA DESCRIPTION AND ANALYSIS 5 5 the observed and calculated traveltimes (the misfit) with respect to their uncertainties, according to X2 = E[i *' }i ( 3 . 5 ) i=l "i where o~i is the standard deviation in the traveltimes. This x2 factor which is subsequently normalized by the number of observations, should equal 1 for a model that reproduces the observed traveltimes within the uncertainty bounds of the picks, x2 values less than 1 indicate that the data have been over-fit and that the model contains structure that is not required by the data. Conversely, x2 values greater than 1 normally mean that the data have sampled heterogeneities in the region that have not been resolved by the model. Other factors such as 3 -D effects, the non-linearity of the source-receiver arrays, misidentification of arrivals and the use of an inappropriate modeling technique (e.g. ray theory) may also contribute to x2 values greater than 1 (Zelt and Forsyth, 1 9 9 4 ) . In such cases, a x 2 factor as close to 1 as possible is sought (Zelt, 1 9 9 9 ) . 3.3.3 Amplitude Modeling When modeling refraction/wide-angle reflection data, it is important that the final velo-city model satisfy both the traveltime and amplitude information provided by the data. However, due to the difficulty in recording well-calibrated amplitudes in land-based expe-riments and the sensitivity of amplitudes to small-scale heterogeneity (Zelt and Forsyth, 1 9 9 4 ) , the model should be constructed so that it fits the traveltime information first and foremost. For this study, amplitude modeling was used solely as a secondary check to constrain velocity gradients in regions with poor ray coverage. The R A Y I N V R algorithm deals with the traveltime information and provides the user with a model that satisfies those constraints. To obtain as accurate a velocity model as Chapter 3. DATA DESCRIPTION AND ANALYSIS 56 possible, the amplitude information should be incorporated into the modeling procedure to fine-tune the model once the traveltimes have been suitably modeled. This step can be done using the T R A M P algorithm, also developed by Zelt and Ellis (1988). The T R A M P algorithm is an efficient 2-D amplitude forward modeling program that allows the user to construct synthetic seismograms using the same model parameteriza-tion as R A Y I N V R . The T R A M P algorithm calculates the amplitudes for the synthetic seismograms using a quality factor (Q) which accounts for attenuation of the seismic pulse through the transmitting medium and which is equal to 27r/fraction of energy lost per cycle. For this study, Q was assigned a value of 600 based on a scattering model for earthquakes (Singh and Herrmann, 1983) which determined Q values of 600-800 for the foreland region of the Cordillera. Using the routine plot sec _amppk, amplitude values were calculated for the observed data by summing the amplitudes over a 100 ms window following each pick. The time window for the summation was determined by examining clear wavelets from shot 2108. These wavelets were windowed and the first three cycles of the wavelets were stacked to provide a source wavelet for the synthetic seismograms. P L T S Y N then convolved the source wavelet with an impulse response and scaled the resulting wavelet by the amplitude values calculated by T R A M P to produce synthetic seismograms. Chapter 4 D E V E L O P M E N T O F T H E V E L O C I T Y M O D E L 4.1 Modeling the Refraction Data Modeling of refraction data using R A Y I N V R requires the implementation of a number of key steps which are illustrated in the flow chart presented in Fig. 4.1. Firstly, a preliminary 2-D velocity model is constructed based on results from preliminary modeling techniques (r — p analysis) and based on a priori information from surface geology. The preliminary model is then discretized according to the required parameterization for R A Y I N V R . This discretization involves dividing the 2-D model into layers defined by boundary nodes and then assigning velocities at points along the top and bottom of each of the layers. Once the preliminary model is parameterized, the R A Y I N V R procedure is run using a layer-stripping approach. This approach involves ray tracing in the top layer of the model (closest to the surface) to attempt to fit the near offset Pg observed arrivals first. In order to fit the observed arrivals, the inversion routine is run and the model perturbations required to fit the arrivals are added to the model. To ensure that the resulting perturbations are geologically reasonable, the new model is examined by the seismologist and adjusted, if necessary, by adding or removing boundary and/or velocity nodes or by changing their resolved values by hand. The inversion is then rerun for the layer and the resulting model is readjusted by the seismologist until an acceptable fit is achieved. The iterative use of the inversion step and the subsequent forward modeling 57 Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 58 Picking of first arrivals Tau-p analysis Parameterization of Initial RAYINVR Model Proceed to next layer X Forward Ray Tracing Damped Least Squares Inversion Pick Reassessment f Produce AVO Plots and Synthetic Seismographs FINISHED! YES Figure 4.1: Flow chart for modeling of refraction data. Dashed lines indicate optional steps. adjustments prevent the addition of geologically unreasonable model perturbations and reduce the amount of time required to achieve a reasonable velocity distribution. Once the rays traced through the top layer reproduce the observed traveltimes to a satisfactory degree, the first layer is kept fixed and the next deepest layer is adjusted to fit the later observed arrivals. This technique is repeated for each layer all the way down to the base of the model. The use of a layer-stripping approach is very important as the PmP and Pn arrivals can be significantly affected by the upper crustal velocity distribution which must be well, resolved before the deeper velocities can be adequately modeled. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 59 Since the results from R A Y I N V R are directly dependent on the observed traveltimes picked by the seismologist, the accuracy of the picks is of supreme importance. As the modeling procedure advances, pick reassessment is often required as information from neighboring shots can reveal inaccuracies in the choice of the original picks. Although the adjustment of picks can be done at any point during the modeling procedure, the practice should be kept to a minimum to avoid forcing the data to fit the model. Once a suitable velocity model is constructed using the traveltime information, T R A M P is used to calculate the amplitudes of rays reaching the surface of the model. Using the calculated amplitudes, amplitude values are interpolated for receivers spaced at 1 km in-tervals along the top of the model and the amplitude information is then used to generate synthetic seismograms using P L T S Y N . The resulting synthetic seismograms are shown in Appendix A from Fig. A.3a, A.6a,..., A.24a. Since the velocity model required to fit the traveltimes for Line 21 was so laterally complex, few changes in velocity gradients could be made without adversely affecting the traveltime fits. For this reason, velocity gradients were only altered to satisfy the amplitude data in the lower crust where ray coverage was poor. To obtain a somewhat quantitative view of the effectiveness with which the amplitudes are modeled, a graphical comparison of the amplitudes for the observed picks and for the synthetic seismograms can be employed. In this study, amplitude vs. offset (AVO) plots (presented in Appendix A from Fig. A.3b, A.6b,..., A.24b) were created for each phase and for both the observed picks and the calculated seismograms. This technique provided a quantitative means of determining how well the constructed velocity model reproduced the real amplitude data. Given that only trends in amplitude with offset are required for this type of comparison, the AVO data were smoothed using a gentle smoothing subroutine provided by Ulrych (Ulrych, 1999; personal communication). To avoid high frequency variations in the A V O plot trends, the smoothing routine was iteratively applied Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 60 to the data by a number of times that was arbitrarily made equal to the number of picked traces. Despite the lack of fine adjustments to the velocity model during this stage and the limited ability of ray theory to reproduce amplitude data, the trends from the A V O plots for all shots (except 2104 where the trends appear to be opposite) are generally quite similar. The best matches for all phases occur at shot 2109 at the northwestern end of the line with trend similarity diminishing southeastward. 4.1.1 Detailed modeling results To provide a detailed examination of how the velocity structural model for Line 21 was constructed, the following sections describe the modeling procedure for each portion of the model. The model parameterization of the final model is shown in Fig 4.2 for reference. It should be noted that whenever the inversion routine was employed, a damping factor of 1.0 was used to ensure that the calculated traveltimes were fit as closely as possible within the uncertainty bounds of the observed traveltimes. The maximum al-lowable calculated uncertainties for each iteration were assigned as 0.1 km/s for the velocity nodes and 1.0 km for boundary nodes. These values were kept constant for all inversions throughout the modeling procedure. A . Upper Crustal Layers The uppermost layer of the model is a thin layer (< 1 km) used to represent the surficial Phanerozoic sediments in the region. The initial model consisted of 8 velocity nodes at the top and the bottom of the layer situated immediately beneath the shots points where the ray coverage was concentrated. The top of the layer was assigned 26 boundary Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 61 ; ! : ! ! • ! ! ! ' ! • ! !! Vertical Exaggeration 4.25:1 70 J 1 —t ! —r* ' — i — ' 1—! f '[ ' 1 l I 1 i 0 50 100 150 200 250 300 350 400 450 Distance (km) Figure 4.2: Nodal parameterization of final velocity model for Line 21. Velocity nodes are designated by white inverted triangles and boundary nodes appear as darkened circles. Black inverted triangles indicates points where both velocity and boundary nodes are specified. Layer numbers are indicated in circles for each layer. Solid lines denote the horizontal layer boundaries which connect the boundary nodes while dashed lines are the vertical boundaries assigned by RAYINVR. nodes at a spacing of ~ 20 km to model the topography, which varied up to 2.5 km over the length of the profile in the region; these boundary nodes were kept constant throughout the modeling procedure. The top layer was modeled using only inversion; however, subsequent modeling of deeper layers required the addition of more velocity nodes to the top layer to account for rapid variations in the Pg branches which were caused by near-surface features. Final modeling results for Layer 1, which were based solely on traveltime fits, indicated that the layer was strongly laterally heterogeneous with a range of surface velocities from 1.9 km/s to 5.45 km/s and a range of vertical velocity gradients from 0.2 s _ 1 to 3.5 s - 1 . In addition, a total of 10 boundary nodes were required to keep the layer relatively thin Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 62 with respect to the varying topography. Layer 2 was designed to model the sediments lying between the Phanerozoic cover and the Precambrian crystalline basement. For the easternmost part of the line, the model generally followed the cross-section of sedimentary layers shown in Fig. 1.7. The base of Layer 2 was fixed beneath the locations of shots 2101 and 2102 at 2.8 km below sea level. Immediately to the west of the location of shot 2102, the Bovie Lake Fault is interpreted to offset vertically the crystalline basement from 2.5 km to 5 km with respect to the basement at shot 2102's location. For this reason, the base of Layer 2 was fixed at 4.5 km below sea level between the locations of shots 2102 and 2104. Since no other basement structure depth measurements were readily available for the rest of the line, Pg arrivals at near offsets for the remaining shots were used to determine a reasonable depth for the western portion of layer 2. From these results, 7 boundary nodes were assigned such that the layer thickness remained ~ 1 km. The drastic shallowing of the base of Layer 2 beneath shot location 2104 is required by the high apparent velocities from near-offset Pg arrivals (Fig. A.7). These arrivals suggest rocks with velocities of 6.4 km/s within 5 km of the surface indicating uplifted or upthrust high velocity basement material. Pg arrivals at greater offsets to the east arrive significantly later, having traveled through the increased thickness of slower sediments. In terms of velocity nodes for Layer 2, the initial model consisted of 8 velocity nodes along both the top and bottom of the layer. These nodes were located immediately below the shot points where ray coverage is densest. While this layer was also fit using only the inversion algorithm, modeling of deeper layers required several additional adjustments to the nodal density for Layer 2. In the final velocity model, Layer 2 demonstrates a high degree of lateral heterogeneity. The parameterization of the layer involves 18 velocity and 10 boundary nodes along the top of the layer and 21 velocity and 12 boundary nodes along the bottom of the layer. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 63 Velocities in the layer vary between 2.42 km/s and 6.40 km/s while vertical velocity gradients vary between 0 s _ 1 to 1.8 s _ 1 . Extending ~ 9 km below sea level, Layer 3 represents the deepest layer within the upper crust of the model parameterization. Rays traveling through this layer account for 1074 calculated Pg arrivals; full Pg ray coverage is shown in Fig. 4.3. As with the shallower layers, the parameterization for Layer 3 in the initial model involved 8 velocity and depth nodes across both the top and bottom of the layer. When the inversion was run, it was found that the lateral heterogeneity in the Pg traveltimes could not be fit by the coarse initial model produced by r — p analysis. i 1 r 0 50 100 150 200 250 300 Distance (km) 350 400 450 Figure 4.3: Ray coverage for Pg phase. Shot locations are numbered and shown as stars across the top of the model. Model layers are shown and numbered. In particular, the slowing of Pg arrivals to the west for shot 2108 over a narrow offset range was found to be too drastic to be accounted for by near-surface features. Consequently, the modeling of the slowed Pg arrivals required a complex parameterization of a narrow low velocity region that extended up to the surface. Given the overlap of ray coverage from neighboring shots, a low velocity notch between the locations of shots Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 64 2108 and 2109 was constrained by Pg arrivals from shots 2106, 2107, 2108 and 2109. Due to R A Y I N V R ' s difficulty in dealing with near-vertical model boundaries, the inversion process became quite unstable when attempting to model the notch. For this reason, the western portion of Layer 3 was modeled using primarily forward modeling with some inversion. To account for Pg arrival branches from shots 2101, 2102 and 2104, in the eastern portion of Layer 3, a sharp vertical juxtaposition of fast (6.4 km/s) and slow (5.8 km/s) material was required just east of shot location 2104 at the Cordilleran deformation front. This juxtaposition was resolved by the inversion and fine-tuned with forward modeling to ensure that rays were traced to all observations. In the final model, Layer 3 is parameterized by 20 velocity and 12 boundary nodes along its top and 16 velocity and 8 boundary nodes along its base. Velocity values vary between 4.74 km/s in the low velocity notch to 6.49 km/s at the Cordilleran deformation front. Vertical velocity gradients within the layer are very low. Amplitude constraints were considered once a suitable fit of traveltimes had been achieved for rays turning in the top three layers of the model parameterization. However, due to the complex nodal density required to reproduce the traveltimes in these layers, few adjustments could be made to satisfy the amplitude constraints without adversely affecting the traveltime fits or without preventing rays from being traced to all observation locations. B. Middle Crustal Layers The middle crust in the velocity model is modeled by Layer 4 which is approximately 8 km thick. Traveltimes in this layer are constrained by 104 far-offset Pg arrivals from the shots at the ends of the line. Few Pg constraints are available for Layer 4 at the Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 65 middle of the profile so amplitude constraints were used to control velocity gradients in that region. The boundary between layers 3 and 4 is modeled as a continuous velocity boundary with a slight change in velocity gradient across the boundary to account for a slight bend in Pg arrivals at far offsets. 'While the initial model for Layer 4 involved 8 velocity and boundary nodes across the top and bottom of the layer, the final model resulted in 16 velocity nodes across the top of the layer. This increase in nodal density was required to maintain continuity in velocities with the 16 velocity nodes along the base of Layer 3. At the base of the layer, the number of velocity nodes was reduced to 15 as the extra node provided no statistical improvement in the traveltime fits. In terms of boundary nodes, these were reduced from 8 at the top to 5 at the base since no middle crustal reflections were available to justify a more complicated basal topography. Velocities in the final model range between 5.62 km/s to 6.49 km/s for Layer 4; velocity gradients appear insignificant. Given the overall lack of Pg constraints for this layer, velocities were constrained by the available Pg arrivals, amplitude modeling results and the rays reflected from the Moho. To avoid adding unnecessary structure to the model, velocity variations in Layer 4 were kept to a minimum. C. Crustal Reflectors in the Model The disconnected short reflections, R, recorded from most of the shot points, are shown in Fig. 4.4 with their corresponding ray coverage. These reflections did not appear to be laterally continuous across any of the data sections and as such, could not have resulted from a laterally extensive crustal reflector. To incorporate the structural information pro-vided by these anomalous reflections into the velocity model, the reflectors were assumed to not be associated with any specific velocity discontinuities at the model boundaries. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 66 Thus, they were modeled as floating reflectors that could occur at any point in the model, regardless of layer boundaries. These floaters provide some insight into structures in the crust but are not associated with velocity discontinuities in the model. 7 0 J 1 1 — i 1 1 1 1 1 1 <— 0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 4 0 0 4 5 0 Distance (km) Figure 4.4: Ray coverage for R phases. Shot locations are numbered and shown as stars across the top of the model. Model layers are shown and numbered. D. Lower Crustal Layers and the Moho Layer 5 is approximately 15 km thick and extends to the base of the Moho from the loca-tion of shot 2104 westward and down to a pinched out transition zone at the easternmost part of the line. Velocity constraints from Pg arrivals for Layer 5 are restricted to only 14 Pg arrivals, most of which are generated by shot 2101. Lacking other traveltime con-straints such as lower crustal reflections, the boundary between layers 4 and 5 was again modeled as a continuous velocity boundary and PmP amplitudes were used to constrain the velocity gradients. To ensure the continuity in velocities across the layer boundary, the top of Layer 5 was assigned the same number of, and values for, velocity nodes as those assigned at the base of Layer 4. To determine the depth to the base of Layer 5 and the velocities along that boundary, PmP arrivals (ray coverage shown in Fig. 4.5) were used to invert both velocity and Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 67 60 4 ® 70 300 350 400 450 0 50 100 150 200 250 Distance (km) Figure 4.5: Ray coverage for PmP phase. Shot locations are numbered and shown as stars across the top of the model. Model layers are shown and numbered. boundary nodes. Since the reflector R l at the base of Layer 5 in the eastern portion of the line is believed to have masked the later PmP arrivals (to be discussed in next subsection), the velocity and boundary nodes for the base of Layer 5 in the easternmost portion of the line were resolved by inverting the R l arrivals (ray coverage shown in Fig. 4.6). With the lack of clear PmP arrivals east of shot location 2104, the Moho could only be firmly constrained for the western part of the line and the eastern portion had to be constrained by Pn arrivals. While it is common practice to invert for Moho depths using both PmP and Pn arrivals, only PmP arrivals were used for the inversion for Moho depths along Line 21. This practice was done to avoid biasing the Moho depths with the laterally varying Pn traveltimes. The base of Layer 5 in the final velocity model was constrained by 12 velocity nodes and 14 boundary nodes, from which 10 boundary nodes outlined the Moho in the west. Velocity values at the base of the crust varied between 6.50 km/s and 7.08 km/s. The vertical velocity gradient across Layer 5 remained constant at 0.1 s _ 1 . Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 68 70 J 1 1 - 1 1 1 ' 1 1 1 0 50 100 150 200 250 300 350 400 450 Distance (km) Figure 4.6: Ray coverage for R l phase. Shot locations are numbered and shown as stars across the top of the model. Model layers are shown and numbered. E . Upper Mantle Layer From Line 21, 1023 Pn arrivals (full ray coverage shown in Fig. 4.7) were picked from all shots except 2104 where they were not distinguishable. These arrivals were clear at offsets up to 435 km at the ends of the line, providing excellent ray coverage to considerable depth in the upper mantle. Given these many arrivals and the confidence in the crustal model, the inversion of the Pn arrivals was expected to be very straightforward. During the initial stages of modeling the lower crust and the upper mantle along Line 21, the R l reflector at the easternmost end of the line had been interpreted and picked as a Moho reflection (PmP). The upper mantle was parameterized as one layer with 11 velocity nodes along the top of the layer and 9 velocity nodes along the base. The top of the mantle was assigned a constant velocity of 8.0 km/s while the base was assigned a constant value of 8.4 km/s. From here, the inversion was set into motion. Complications arose when the initial inversion results revealed that the top of the mantle beneath shot location 2104 required velocities of 8.5 km/s and higher while the top of the mantle beneath the locations of shots 2101 and 2102 required upper mantle Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 69 1 1 i i i i 1 i 0 50 100 150 200 250 300 350 400 450 Distance (km) Figure 4.7: Ray coverage for Pn phase. Shot locations are numbered and shown as stars across the top of the model. Model layers are shown and numbered. velocities as slow as 7.3 km/s. In addition to these highly variable velocity results, the model was incapable of tracing the clear Pn arrivals seen from shot 2105 (Fig. A . 10). Since rays must be traced to all observations for a model to be acceptable, an extra high velocity lower crustal layer (Layer 6) was added to the model beneath the eastern-most portion of the line to attempt to obtain the Pn rays from shot 2105. With the addition of a 3-4 km thick pinched out layer with an average velocity of 7.3 km/s, the Pn rays were finally traced although their arrivals were still faster than the observed Pn phase (Fig. 4.8) since they still had to travel through the higher upper mantle veloci-ties beneath shot location 2104. These higher velocities beneath the Foreland Belt were required to fit the Pn arrivals from all the other shots. The addition of a pinched-out high velocity lower crustal layer at the eastern end of the line allowed the velocities at the top of Layer 7 to be increased to more reasonable values of 7.7 km/s in the east. The geometry of the pinched-out zone was determined by inverting the lower boundary nodes of Layer 6 using the Pn arrivals from Layer 7. By definition, velocities less than 7.6 km/s cannot exist below the Moho. Thus, Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 70 14 12 10J co8 Q 2H PmP mmiugajia**" Pn ^ Pg Distance (km) 0 50 100 150 200 250 300 350 400 450 Offset (km).200 -150 -100 -50 0 50 100 150 200 250 Figure 4.8: Comparison of observed (vertical bars with length equal to twice the pick un-certainty) and calculated (thin lines) traveltimes from shot 2105. A l l other identified phases are labeled. Note how the Pn traveltimes calculated from the final velocity model still appear before the observed phase even with the low velocity zone in the upper mantle between shot locations 2102 and 2104. the PmP-type reflections from shots 2101, 2102 and 2104 had to be reassessed as being caused by a lower crustal reflector (Rl ) that existed a few kilometers above the Moho and that masked the later PmP arrivals in its coda (example from Shot 2101 shown in Fig. 4.9). This interpretation is consistent with the complex PmP-type reflections seen on the easternmost shots. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 71 Figure 4.9: Trace normalized plot for shot 2101 showing PmP phase masked by earlier R l arrivals. Inset is an enlargement of data within the box outlined by the thin line shown. Thin horizontal bars on each trace indicate the picked traveltimes. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 72 4.1.2 Spatial Resolution and Absolute Parameter Uncertainty The velocity model that was developed in this thesis was arrived at using a minimum-parameter/prior-structure model approach (Zelt, 1999) which signified that a priori in-formation (e.g. well information on crustal layer depths) was used to construct a coarsely parameterized model with a non-uniform nodal spacing. While the results show that the final velocity structural model fit the observed data to a satisfactory degree, it is im-portant to note that the final model does not represent a unique solution to the inverse problem. This non-uniqueness arises from the model's dependence on the parameteri-zation chosen and also from the inevitable uncertainty of the traveltime picks for the data being modeled. The best direct way to evaluate the non-uniqueness of a given model is to construct several alternate models that fit the data equally well but that have different parameterizations (Zelt, 1999). While this direct technique provides the modeler with a strong impression as to what features are required in the model to fit the data, the excessive amount of time required to construct several models precludes this approach. While direct methods are the most widely used techniques to address the issue of non-uniqueness, indirect methods can be employed to assess the reliability of model parameters in the final velocity model. Zelt and Smith (1992) suggest two tests that can be applied to the model parameters to quantify spatial resolution and absolute parameter uncertainty. The spatial resolution test examines how well an individual node is resolved in the model by determining whether perturbations of that node can be corrected by the in-version step without affecting neighboring nodes. To execute the test, a specific node is chosen and perturbed by an amount on the order of its uncertainty. Once the node is perturbed, rays are traced through the perturbed model and new traveltimes are calcu-lated. These calculated traveltimes are then interpolated to the receiver locations from Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 73 the experiment and assigned the uncertainties of the corresponding observed picks. The perturbed node is then reset to its value from the final model and the calculated tra-veltimes are inverted with respect to the selected node and the other model parameters that were determined at the same time during the inversion for the final velocity model. The spatial resolution of the selected node is determined based on the extent to which neighboring nodes are affected by the inversion step. If the node was perfectly resolved in the final model then the neighboring nodes will not have been affected by the inversion. If that was not the case, then the perturbation will have been smeared into neighboring model parameters. The lateral extent of the smearing corresponds to the lateral spatial resolution of the model. Zelt and Smith (1992) suggest using perturbations large enough to cause significant traveltime anomalies but small enough so that the ray tracing is not strongly affected. For this study, the velocity nodes were perturbed by 0.1-0.2 km/s and the boundary nodes were perturbed by 1 km. The absolute parameter uncertainty test examines how much a particular node can be perturbed before the perturbation affects the inversion so much that the resulting model is significantly different from the final model. For this test, a node is selected and perturbed by an amount on the order of its uncertainty. The perturbed node is then kept fixed while the other model parameters that were determined at the same time as the selected node during the inversion for the final model are inverted with respect to the original observed traveltimes. The resulting model is then compared with the final model to determine whether the models are significantly different based on a) whether rays can still be traced to all observations and b) a comparison of the x2 values and the number of data points using the statistical F-test. If the models are not found to be significantly different, the perturbation is increased and the process is repeated. The minimum perturbation required for the resulting model to be significantly different from the final model represents an estimate of the absolute parameter uncertainty. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 74 The main drawback to the tests suggested is that they assume that the model was achieved using solely inversion. For laterally heterogeneous regions with fairly complica-ted parameterizations such as the region modeled in this thesis, the inversion step can become quite unstable and the inversion must be monitored and adjusted to maintain a model that is geologically reasonable. While the forward modeling approach used to keep the model "geologically honest" often improves the traveltime fits considerably, the resulting addition of nodes to certain sections of the model can cause rays to not be traced to all observations following subsequent inversions. This was particularly the case for the western portion of the model where the low velocity notch was modeled, for the eastern portion of the model where the drill hole information was used to position the boundary nodes at the base of layer 2 and for the upper mantle where the Pn arrivals indi-cated complex lateral heterogeneity. The inability to ray trace in these regions following inversion caused the parameter tests to be inconclusive. Despite these complications, both the spatial resolution and the absolute parameter uncertainty tests were applied successfully throughout the rest of the final model. Given the significant amount of time that would be required to test all model parameters, 2 to 3 representative model parameters were tested for each layer as suggested by Zelt and Smith (1992) and Zelt (1999). The results from both tests are presented in Table 4.1. While the results in Table 4.1 provide an estimate for the reliability of the final velocity model based on traveltimes, the tests applied are unable to take into account the extra constraints added to the model from amplitude forward modeling results. Similarly, due to the overall lateral heterogeneity of the final model, lateral resolution and absolute parameter uncertainty estimates were highly variable for each layer. Thus, the values presented in Table 4.1 should be regarded as worst case scenario estimates. Subsequent forward modeling tests were run on the final velocity model to reassess the lateral resolution and absolute parameter uncertainty while taking into account the Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 75 added constraints from traveltime and amplitude forward modeling. These tests involved assessing the effect of small adjustments of individual velocity and boundary nodes on ray tracing, statistical traveltime fits and amplitude vs. offset plots. For regions of the model where forward modeling was the main modeling approach, the forward modeled lateral resolution and absolute parameter uncertainties are listed in the table with the Zelt and Smith (1992) test results presented in brackets. Node Location Layer Approximate lateral resolution (km) Approximate absolute uncertainty V Surface 1 (T) <15(40)-110 0.1-1.0 (0.1-1.6) km/s V Top of Upper Crust 2 ( T ) <15(40)-110 0.1-0.5 km/s V Within Upper Crust 3 ( T ) 75 0.1-0.2 km/s V Top of Middle Crust 4 ( T ) 60 0.1-0.2 km/s V Base of Lower Crust 5 ( B ) 80 (120) 0.3 km/s B Moho west of 2104 5 ( B ) 60 1.5-2.5 km B R l Reflector 5 ( B ) 100 1.5-2.5 km B Moho east of 2104 6 ( B ) 100 3.5 km V Top of Upper Mantle 7 (T) 90 0.2-0.25 (0.25-0.4) km/s V Base of model 7 (B) 140 0.3 km/s Table 4.1: Results of spatial resolution and absolute parameter uncertainty tests for the final velocity structural model. V and B in the node column refer to velocity and boundary nodes respectively. T and B in the layer column refer to the top and base of the layer. Bracketed values represent results from Zelt and Smith (1992) tests that have been reassessed using forward traveltime and amplitude modeling constraints. Some general conclusions can be drawn from the model assessment test results. Firstly, the lateral resolution of the model in terms of velocity nodes improved from the surface down to the top of the middle crust due to the increase in Pg ray coverage. This improvement in lateral resolution was accompanied by a decrease in absolute pa-rameter uncertainties. Conversely, the velocity nodes at the base of the crust provided lateral resolutions between 80 and 120 km and an increased parameter uncertainty of 0.3 Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 76 km/s based on the PmP traveltimes. Secondly, while the PmP arrivals provided a lateral resolution of 60 km on the boun-dary nodes of the Moho west of shot 2104, fewer reflected arrivals were available for the R l phase resulting in poorer resolution of the R l reflector. Despite this lateral variation in resolution for the base of Layer 5, absolute uncertainty values remained between 1.5-2.5 km across the layer. East of shot location 2104, the lack of PmP arrivals required the depth of the Moho to be constrained solely by Pn arrivals which resulted in a poorer resolution of 100 km and a high absolute uncertainty of 3.5 km for the Moho boundary nodes. Pn arrivals from 7 shots provided sufficient ray coverage in the upper mantle to allow for a lateral resolution at the top of the upper mantle of 90 km which degraded down to 140 km at the base of the final model due to a decrease in ray coverage with depth. Absolute parameter uncertainties from forward modeling tests varied between 0.2 km/s and 0.25 km/s throughout the upper mantle and did not seem to exhibit a depth dependence. The higher uncertainties of 0.25 to 0.40 km/s from the Zelt and Smith (1992) tests for the upper mantle resulted because the inversion step insisted on inserting geologically unreasonable upper mantle velocities of 8.5 km/s and higher immediately beneath the Foreland Belt. Using forward modeling, these velocities were reduced to 8.3 km/s. 4.2 The Velocity Structural Model for Line 21 The velocity structural model for Line 21 (Fig. 4.10) was achieved through an interac-tive combination of traveltime forward modeling and 2-D inversion as well as forward amplitude modeling. The results reveal that the data from Line 21 sampled a laterally heterogeneous region of considerable complexity. It is important to note that the resolved Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 77 model (Fig. 4 . 1 0 ) bears little resemblance to the model parameterization (Fig. 4 . 2 ) , sig-nifying that the modeling procedure was not biased by the parameterization chosen. The detailed velocity model for Line 2 1 was accepted as the best model due to its overall RMS residual traveltime of 1 2 0 ms, its normalized x 2 value of 2 . 9 9 1 and its ability to trace rays to all observation points. The statistical fits achieved for all phases and for all shots are displayed in Table 4 . 2 and Table 4 . 3 , respectively. The ray coverage for the model, which was quite dense, is displayed in Fig. 4 . 1 1 . Despite numerous attempts, a %2 value of 2 . 9 9 1 was the lowest achieved for Line 2 1 that allowed rays to be traced to all observations. This may have been due to mispicking of a phase due to the background noise level and/or the coda of earlier arrivals or due to 3 -D variations in the region that could not be accounted for with a 2-D model, although Zelt and Zelt ( 1 9 9 8 ) suggest that, in general, this should not be a significant problem. While a more complex parameterization or an increase in pick uncertainty may have lowered the \ 2 factor closer to 1 , I considered that the 1 2 0 ms RMS residual traveltime was so close to the overall pick uncertainty of 1 1 2 ms that the introduction of more complexity would have incorporated unnecessary structure into the model. # o f . # o f Pick RMS Travel Normalized Phase Picked Calculated Uncertainties Time x2 Traveltimes Traveltimes (ms) Residual (ms) Misfit Pg 1 3 6 0 1 3 3 6 8 1 1 0 6 4 . 8 0 1 PmP 6 6 0 6 5 8 1 3 0 1 0 5 1 . 2 1 7 Pn 1 0 2 3 1 0 1 8 1 3 7 1 5 2 2 . 4 9 4 R 3 6 1 2 0 7 1 0 3 4 9 1 . 0 7 5 R l 3 0 7 3 0 2 1 2 0 1 1 9 2 . 2 5 3 Total 3 7 1 1 3 5 2 1 1 1 2 1 2 0 2 . 9 9 1 Table 4 . 2 : Model fitting statistics for each phase. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 78 b) WNW 0 50 100 150 200 250 300 350 Distance (km) 400 450 E 0 -20 Omineca Belt 9 8 | 7 Foreland Belt 6 5 Ancestral N.A. 4 I 2 g-40 -60 H 150 200 250 300 Distance (km) ESE (km/s) 1.9 3.0 4.0 5.8 6.0 6.2 6.4 6.8 7.0 7.5 8.1 8.4 8.9 Figure 4.10: Final velocity structural model for Line 21, a) plotted with no vertical exaggera-tion; b) plotted with a vertical exaggeration of 5:2. Shot locations are plotted with stars along the top of the model. Thick lines signify regions constrained by wide-angle reflections. The M indicates the location of the Moho while the dotted M line represents the Moho depths not constrained by wide-angle reflections. Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 79 # o f # o f Pick RMS Travel Normalized Shot Picked Calculated Uncertainties Time x2 No. Traveltimes Traveltimes (ms) Residual (ms) Misfit 2101 521 511 129 112 1.354 2102 560 482 82 165 4.505 2104 243 228 117 104 2.271 2105 411 402 116 137 3.546 2106 505 484 98 108 3.439 2107 585 543 181 98 2.794 2108 459 449 61 114 3.709 2109 427 422 69 103 2.130 Table 4.3: Model fitting statistics for each shot point. 4.2.1 Main Features of the Velocity Model The final velocity structural model for Line 21, which is shown in Fig. 4.10, is strongly laterally heterogeneous and can be roughly subdivided into three main blocks: the edge of cratonic North America, the Foreland Belt and the Omineca Belt. These three main blocks demonstrate differing velocity signatures that provide clues to their individual evolution. Since the resolved model appears independent of the model parameterization used in its development (as it should), the crustal divisions (e.g. upper crust) used in this section to describe the features of the model are independent of the model parameterization layers. As such, regions of similar velocities within the final velocity model are used to divide the model into upper, middle and lower crustal blocks. The upper crust of the model, which extends down to ~ 9 km below sea level, is characterized by a range of velocities from 1.9 km/s to 6.5 km/s. While the range in velocities for the upper crust is quite extensive, the extreme velocities in the range are Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 80 WNW ESE Omineca Belt 9 8 7 0 -20 Foreland Belt 6 5 Ancestral N.A. J--40 -60 A 0 50 100 150 200 250 300 350 400 450 Distance (km) lntensity=function(ray density) Full Colour=64 rays White=0 rays to Q64 £ l!9 3!o 4!o 5!8 6.0 6l2 6l4 8.8 7!0 7l5 8.\ 8A 8.9 (km/s) Figure 4.11: Ray coverage of velocity model for Line 21. Intensity of color ranges from areas of full color which indicate regions traversed by 64 rays or more to white regions which indicate regions with no ray coverage. concentrated in specific blocks of the model. The edge of cratonic North America, which lies at the southeasternmost part of the line, has an upper crust with an average velocity of 5.8 km/s. This region of relatively low velocities thickens westward (from 5 km to 15 km) toward the Foreland Belt where it is abruptly cut off. Meanwhile, the juxtaposed eastern front of the Foreland Belt is characterized by anomalously high velocities of ~ 6.4 km/s within 5 km of the surface. Further to the west in the Foreland Belt, the upper crustal velocities return to more average velocities of 6.0 km/s. To account for the significant slowing of Pg arrivals between shot locations 2108 and 2109, the upper crust of the Omineca Belt required the addition of a trench of low velocity Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 81 material. It is at the surface of this trench that the slowest velocities in the model (1.9 km/s) are found. This trench which extends down to the base of the upper crust, has steep sides and is ~ 15 km wide. The middle crustal region of the model contains velocities ranging between 6.0 km/s and 6.5 km/s. While its thickness at the edge of cratonic North America is approximately 10 km, the middle crust thickens beneath the eastern front of the Foreland Belt to 15 km, thins to 10 km in the western portion of the Foreland Belt and thickens again significantly in the Omineca Belt to almost 20 km. The thickening that occurs in the Omineca Belt is accompanied by a 10 km offset of the 6.4 km/s contour in the middle crust beneath shot 2108 at the approximate location of the Tintina Fault. The velocities in the lower crust of the model range between 6.5 km/s and 7.3 km/s. With the exception of the high velocity (7.3 km/s) pinched out layer beneath the edge of cratonic North America, high velocities at the base of the crust do not appear to be associated specifically with the three main blocks of the model. These high velocity (> 6.8 km/s) regions are located beneath the eastern portion of the edge of cratonic North America, at the Cordilleran deformation front and beneath shot location 2107. : PmP arrivals west of shot 2104 constrain a Moho that exists at ~ 34 km below sea level in the Omineca Belt, deepens slightly to 35 km below shot 2107, shallows again to 34 km below the western portion of the Foreland Belt and deepens to almost 38 km below the eastern front of the Foreland Belt before shallowing again to 34 km below shot 2104. While PmP arrivals are not available east of shot 2104, upper mantle refractions through the high velocity lower crustal pinched out layer beneath the edge of cratonic North America suggest a Moho that deepens to 40 km below sea level at the southeastern end of the line. The upper mantle of the model is highly laterally heterogeneous with a range of velocities between 7.7 km/s and 8.9 km/s and is modeled down to depths of 70 km below Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 82 sea level. The lowest upper mantle velocities are located beneath the edge of cratonic North America where they underlie the high velocity lower crustal layer that pinches out from east to west. Immediately to the west, the upper mantle beneath the Foreland Belt is characterized by anomalously high velocities up to 8.3 km/s. Meanwhile, the velocities in the upper mantle beneath the Omineca/Foreland Belt transition exhibit values ranging between 7.8 km/s and 8.1 km/s. 4.2.2 Velocity Anomaly Calculation Once a suitable velocity structural model is achieved, insight can be derived by calculating a velocity anomaly model. The velocity anomaly model outlines the detailed structure of differing velocity blocks which can help resolve block boundaries (Winardhi and Mereu, 1997). The anomaly model is obtained by averaging the final velocity model at a given depth and subtracting that average value from all points at that same depth in the final velocity model. Complications can arise near the Moho when it has topography. These problems can be avoided by amending the averaging algorithm such that it calculates two averages at a given depth, one lower crustal average and one upper mantle average. This technique works as long as the maximum lower crustal velocity is less than the minimum upper mantle velocity. While the lateral heterogeneity in the velocity model for Line 21 is readily apparent from the model itself, a velocity anomaly model was generated to stress this lateral complexity (Fig. 4.12). Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL 8 3 WNW -5.4 -0.5 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.5 1.8 Figure 4.12: Velocity anomaly model for Line 21. Shot locations are plotted with stars along the top of the model. Thick lines signify regions constrained by wide-angle reflections. The M indicates the location of the Moho while the dotted M line represents the Moho depths not constrained by wide-angle reflections. Chapter 5 P O T E N T I A L F I E L D D A T A M O D E L I N G The final velocity structural model developed in this thesis provides a cross-sectional image of the distribution of P-wave velocities with respect to depth and lateral position along Line 21. While the P-wave velocities allow some speculation on the rock types making up the field area, the use of other independent geophysical information can im-prove the interpretation by providing extra constraints on other physical parameters of the crust. While few geophysical surveys have collected data in the region of Line 21 (near-vertical reflection and magnetotelluric data will be collected in the winter of 1999), detailed gravity and aeromagnetic field measurements are available. Analyses of these potential field data can provide density and magnetic susceptibility constraints to aid in the interpretation of Line 21. 5.1 Gravity Data and Modeling The Bouguer-corrected gravity data with a crustal rock density of 2.67 g/cm3 were ac-quired in digital format from the Geophysical Data Centre (Geological Survey of Canada, Ottawa) and were gridded at ~ 4 km spacings using a continuous curvature surface grid-ding algorithm (Fig. 5.1). The main feature of the Bouguer gravity map is the prominent change from relative high values above the W C S B and underlying Precambrian domains to the relative low values associated with the Cordillera. Two additional features of 84 Chapter 5. POTENTIAL FIELD DATA MODELING 8 5 relevance to this study are: 1) the Tintina Fault which is defined by a gravity expres-sion which shows a narrow increase in gravity (relative to values to the N E and SW) associated with the trace of the fault; and 2) the region of the western end of Line 21 where higher values are observed. These higher values extend E-W along latitude 60°N to the Cordilleran deformation front and appear distinct from the Cordillera and from the region to the east. Figure 5.1: Bouguer gravity map of region surrounding Line 21. The refraction line and all shot points are overlain on the plot for reference. The dark grey box within the map delineates the portion of the data used in the GRAV3D inversion. The traditional approach to modeling the Bouguer gravity signal along a refraction line has been to sample the gravity values along the line and model them using a 2.5-d modeling algorithm such as S A K I (Webring, 1985). To obtain meaningful results from such a procedure, the extracted line of gravity values should extend perpendicular to Chapter 5. POTENTIAL FIELD DATA MODELING 86 the strike of the main anomalies and the density blocks making up the model should be extended either side of the line to a finite or infinite extent based on a priori information from other geological and geophysical sources. Simple examination of the Bouguer gravity map (Fig. 5.1) reveals that Line 21 was not an optimal candidate for 2.5-d modeling due to its orientation oblique to gravitational strike and the lack of a priori regional information about underlying structures. For this reason, preliminary 3-d gravity inversion modeling was carried out to provide an indication of density variations and determine the offline extent of anomalous densities at depth. 5.1.1 3-d Inversion Modeling using G R A V 3 D The gravity field measured at a specific location is the linear combination of the vertical gravity effects generated by bodies of anomalous density at depth. The traditional ap-proach to interpreting deep density distributions has involved the visual examination of gravity maps to infer depth and lateral extent of the density blocks that caused specific gravitational anomalies. Due to the inherent lack of depth resolution in gravity data, this approach results in density models that are highly non-unique. Without a priori knowledge of deep structure from other techniques (e.g. seismic refraction, seismic reflection, geological information, drilling information), visual inter-pretation of gravity maps can only resolve large scale density trends and can provide little information about the depth and lateral extent of density blocks. As such, little meaningful interpretation of deep density structure can be inferred by such an approach. The GRAV3D modeling algorithm (Li and Oldenburg, 1998), used for the preliminary modeling of deep density blocks in the region, inverts the observed gravity data to obtain a 3-d density distribution below the observation locations. Since no depth resolution can Chapter 5. POTENTIAL FIELD DATA MODELING 87 be inferred from gravity data, a weighting function is applied to the resulting density distribution to account for the natural geometric decay of the resolution kernels with depth and as such, to prevent the inversion from concentrating the density anomalies at the surface of the model. This weighting function approximates the inverse distance-squared decay of the resolution kernel. The inversion is formulated as an optimization problem which allows a specific type of model (in this case, the smoothest model with the least amount of structure) to be constructed while fitting the data within its error bounds. While the resulting density model does not qualify as a unique solution, it does produce a physical model consistent with the gravity observations which can be used to provide constraints on deep density structure. The observed gravity data used in the inversion were a subset of the data shown in Fig. 5.1. These windowed data, which lay between latitudes 57°N and 61°N and between longitudes 131°W and 122°W, were projected onto an eastings and northings grid system (600 km by 600 km) whose origin lay at 57°N and 131°W. Since the regional field was unknown, a planar trend was fit to the data and a constant regional field of -103.6 mGal was further removed from the Bouguer corrected data. To ensure that the model was not underparameterized, all gravity values within a given cell were averaged and the averaged gravity values were then reassigned to the observation locations in the cell (Fig. 5.2). The mesh onto which the density distribution is modeled is composed of rectangular prisms with a constant density anomaly assigned to each prism. The mesh used to construct the prehminary 3-d density distribution consisted of 42 cells in the easting direction by 42 cells in the northing direction by 20 cells deep. The dimensions of the model cells were 15 km by 15 km by 5 km deep with several notable exceptions. Firstly, the cells at the edges of the mesh were assigned widths of 100 km to act as padding cells; secondly, to incorporate topographical information for the region, the Chapter 5. POTENTIAL FIELD DATA MODELING 88 500 400 ^ 300 200 100 100 200 300 400 Eastings (km) 500 600 3pa| (mgals) -160 -90 -40 -30 -20 -10 0 10 20 30 40 120 Figure 5.2: Averaged observed Bouguer gravity map with regional field removed. The refrac-tion fine is overlain on the plot for reference. Individual cells are indicated by the small squares. top 5 cells of the mesh were each only 500 m deep; and lastly, the deepest layer of cells in the mesh had dimensions of 15 km by 15 km by 220 km deep. These deeper cells acted as padding cells to prevent any influences from densities at greater depths. Given the sparsity of observation locations in the northernmost row of cells, the resolved anomalies in these cells had less constraints than those to the south. To avoid interpreting these unconstrained anomalies, all plots are presented on a 500 km by 600 km mesh corresponding to the region with the most data constraints. To incorporate the topography of the region into the mesh, GRAV3D requires the construction of a special file that lists how many blocks at the top of the mesh must be removed to mimic the topography. In other words, if we assume that the top of the mesh corresponds to a height of 2500 meters above sea level and the top 5 cells of the mesh are each 500 meters deep, a region of the model where the topography is 2000 meters Chapter 5. POTENTIAL FIELD DATA MODELING 89 requires the removal of one ceU off the top of the mesh. Fig. 5.3 shows the topography for the region and the corresponding number of cells removed from the top of the mesh. b), 100 200 300 400 500 600 I (blocks removed) 200 300 400 500 600 0 I | | q (meters) ITT" -2200 0 400 700 900 1100 1300 1500 1800 2600 0 * j Figure 5.3: Topographical information for the Line 21 region; a) the topographic map of the region; b) the number of cells removed from the GRAV3D mesh to model the topography. The refraction line is overlain on the plots for reference. Preliminary inversion results produced a density distribution that reproduced the observed gravity values well (compare Fig. 5.2 with Fig. 5.4), as the algorithm requires. However, examination of the predicted density model (Fig. 5.5) revealed that even with a depth weighting applied, the algorithm concentrated higher density bodies near the surface. These results contradict the physical requirement of deep crustal scale studies that density generally should increase as a function of depth. This discrepancy, which results from the fact that the L i and Oldenburg (1998) code was designed for small scale (less than 10 km) shallow surveys in the mineral exploration industry, can be overcome by using a reference model. To ensure that the model generated by the GRAV3D code would respect the physical requirement that density increase as a function of depth, a reference model (Fig. 5.6) Chapter 5. POTENTIAL FIELD DATA MODELING 90 200 300 400 Eastings (km) I (mgals) -160 -90 -40 -30 -20 -10 0 10 20 30 40 120 Figure 5.4: Initial predicted Bouguer gravity map. The refraction line is overlain on the plot for reference. was created against which the resolved model would be tested. The incorporation of a reference model forces the algorithm to construct a model that fits the observed data while minimizing differences between the reference and the constructed models. The reference model for the Line 21 region was constructed such that blocks at and below 35 km depth (i.e., approximately below the Moho) were assigned density anomalies of 0.53 g/cmz with respect to a reducing density of 2.67 g/cm3 (i.e. the upper mantle was assigned a density of 3.20 g/cm3). Within the crust, the assigned density anomalies increased linearly with depth from -0.37 g/cm3 at the surface up to 0.33 g/cm3 at the base of the crust corresponding to a density range of 2.30-3.00 g/cm3 for the crust. This range of densities was selected using the range in P-wave velocities resolved along Line 21 which were translated into densities using the P-wave velocity/density relationships presented in Barton (1986). Chapter 5. POTENTIAL FIELD DATA MODELING 91 50 100 150 200 250 300 350 400 450 500 550 600 Distance (km) (g/cm**3) -0.085 -0.015 -0.010 -0.006 -0.003 0.000 0.003 0.006 0.010 0.015 0.038 Figure 5.5: Slice through initial density anomaly model at 300 km northing. 50 100 150 200 250 300 350 400 450 500 550 600 Distance (km) (g/cm**3) -0.7 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Figure 5.6: Reference density anomaly model used in the inversion. Chapter 5. POTENTIAL FIELD DATA MODELING 92 Final Inversion Results The initial run of the inversion using the reference model reached 34 iterations before the program determined that the stopping criteria (a \ 2 factor of 1.0) of the inversion couldn't be reached. To avoid altering the uncertainties in the observed data which were provided by the Geological Survey of Canada, the desired x2 factor of the misfit was simply increased to a value of 8.0 which is on the order of the misfit achieved during the 34th iteration of the inversion. Once the %2 factor was readjusted, the inversion was rerun. With the increased %2 factor, the algorithm converged after 32 iterations. The pre-dicted data (Fig. 5.7b) generated by the inversion successfully reproduced the observed gravity data, as required by the algorithm, while also creating a geologically reasonable density anomaly model. Examination of the map of differences between the observed and predicted data (Fig. 5.7c) shows a fairly random distribution around most of the region surrounding Line 21. This randomness indicates that, for most of the region of interest, the inversion successfully reproduced the observations without introducing large spurious anomalies and that the inversion was not fitting the noise in the data. Figure 5.7: Comparison of observed and predicted Bouguer gravity maps from the GRAV3D inversion; a) observed Bouguer gravity map; b) predicted Bouguer gravity map; c) map of the differences between the observed and predicted Bouguer gravity maps. The refraction line is overlain on all plots for reference. Chapter 5. POTENTIAL FIELD DATA MODELING 93 While the calculated density anomaly model does not provide strict density con-straints on the region, it does provide insight into the relative variations of the density distribution at depth as shown by the depth slices of Fig. 5.8. Within 5 km of the surface, the density model exhibits a density anomaly contrast between the eastern portion of the region which lies on the edge of ancestral North America and the less dense Cordillera to the west. This density contrast persists at all depths to 65 km, suggesting that the crust and upper mantle of ancestral North America and those of the northern Cordillera are quite distinct. Curiously, the eastern portion of Line 21 extends onto the edge of cratonic North America into a zone of anomalously low densities compared to the surrounding density anomalies, a characteristic which may reflect the uneven geometry of the rifted margin of ancestral North America. A vertical slice through the model at the exact location of the modeled seismic re-fraction line (Fig. 5.9) was taken to examine whether the resolved velocity blocks cor-responded with distinct regions of anomalous densities. Comparison of the final velocity structural model (Fig. 4.10) with the GRAV3D derived density model (Fig. 5.9) reveal encouraging results. The west facing ramp of higher velocity material along the edge of ancestral North America in the velocity model coincides with a west facing wedge of high density anomalies in the density model. Similarly, a column of low density anomalies in the upper crust between shots 2108 and 2109 mimics the low velocity wedge resolved in the velocity model. At greater depths, the upper mantle exhibits lower density anomalies beneath shot 2102 which correspond with the location of low upper mantle velocities in the velocity structural model. The GRAV3D inversion was intended as a preliminary step in resolving the lateral extent of density anomalies. Even without incorporating the seismically resolved Moho depths, the inversion results show general consistency with the final velocity structural model. Chapter 5. POTENTIAL FIELD DATA MODELING 94 I (g/cm"3) 3.7 -0 .5 -0.4 -0.3 -0 .2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Figure 5 8- Depth Slices through the density anomaly model generated by G R A V 3 D . The shallowest slice appears at the upper left corner while the deepest slice appears at the lower right corner. The refraction line is overlain on all the slices for reference. Chapter 5. POTENTIAL FIELD DATA MODELING 95 0 100 200 300 400 500 600 Eastings (km) P P W H P P P - T - ' J ^ (mgals) -160 -90 -40 -30 -20 -10 0 10 20 30 40 120 W N W Figure 5.9: Slice through GRAV3D density anomaly model along modeled refraction line. The top map shows the location of the slice (black line) with respect to the actual refraction line (grey). The Moho depth determined from the refraction data modeling is overlain on the density anomaly model. Areas with no ray coverage in the corresponding velocity structural model have been masked out to allow visual comparison with Fig. 4.10. Chapter 5. POTENTIAL FIELD DATA MODELING 96 5.1.2 2.5-d Forward Modeling using SAKI To determine whether the observed gravity field values along the refraction line could be reproduced with a density model based on lithological, density anomaly and struc-tural constraints derived from the seismic refraction model and the 3-d gravity inver-sion, 2.5-d forward gravity modeling was carried out using the S A K I modeling program (Webring, 1985). The final velocity structural model (Fig. 4.10) was translated into a density model using reasonable experimental P-wave velocity/density relationships (Bar-ton, 1986), while attempting to maintain the relative density anomalies resolved from G R A V 3 D . Both the surface topography and the Moho topography determined from PmP wide-angle reflections were included in the model as constraints. The density model parameterization used in SAKI involves a polygonal parameteri-zation similar to that of Zelt and Smith (1992). This parameterization divides the model into polygons of finite or infinite lateral extent (perpendicular to the plane of the section), each of which are assigned a density. To avoid edge effects, the model was extended an extra 2000 km at either end and 200 km at its base. The calculated gravity response of the density model was compared with the observed gravity data along the line to assess the reliability of the density model. The observed gravity data to be modeled were extracted from Fig. 5.1 along the same line as that used to model the refraction data. To assess whether the chosen line represented the local field values well, eight additional parallel lines, four to the northeast and four to the southwest, were derived from the gravity map. Data along these lines were averaged point by point to determine the average gravity values along the swath. The trends shown in Fig. 5.10 demonstrate how the main line approximates the average quite well, and to some extent represents the optimal choice given that the averaged line does not resolve the narrow gravity highs and lows in the western portion of the line Chapter 5. POTENTIAL FIELD DATA MODELING 97 which coincide with the narrow wedge of low velocity material interpreted in the velocity model between shot locations 2108 and 2109. Forward Modeling Procedure and Results Forward modeling of gravity data using S A K I is a slow and non-unique procedure given that there exists a wide range of experimentally acceptable densities for any given velo-city. As such, a variety of physically reasonable densities parameterized into any number of polygons can be constructed into a model that will reproduce the observed Bouguer gravity profile. Recognizing this inevitable non-uniqueness, the S A K I 2.5-d model was constructed by initially breaking the model up into polygons that roughly corresponded with areas of specific velocities in the velocity structural model (Fig. 4.10). Once these polygons were created, they were assigned density values close to the mean of experimen-tally derived densities for their given velocity range (Barton, 1986). In order to reproduce the observed Bouguer gravity profile, the polygonal structures and density values were iteratively varied using small changes until an overall fit was achieved. Once the main trends were fit, small near-surface anomalies were inserted in the model to reproduce any short wavelength anomalies in the Bouguer gravity profile. The resulting 2.5-d density model is shown in Fig. 5.11; the comparison of observed and calculated values is shown above the model. Chapter 5. POTENTIAL FIELD DATA MODELING 98 -160.0 ' 1 1 ' • 1 ' ! ' '—1 -40.0 60.0 160.0 260.0 360.0 460.0 Model Distance (km) Figure 5.10: Comparison of gravity trends along parallel lines; a) the contoured Bouguer gravity map overlain by the actual refraction line (thick grey line), the modeled refraction line (thickest black line) and the 8 extra parallel lines (thin black lines) that were sampled; b) a graphical comparison of the Bouguer gravity values extracted along all 9 lines. In the legend, northl signifies the nearest line sampled to the north of the modeled refraction line while north4 is the northernmost sampled line and so forth. Shot point locations are indicated along refraction line for reference. Chapter 5. POTENTIAL FIELD DATA MODELING 99 -80 -90 |-100 E -110 -120 ^ 20 E o. 40 5 • 60 ii ii n ii H II in * Observed - Calculated 1 0 0 Model Distance (km) 200 300 400 B l 2.73 275 2.80 2.82 2.85 j 2.90 2.95 3.00 3.10 3.15 3.35 g/cc Figure 5.11: 2.5-d Density model derived for Line 21; top plot shows match between the observed (+ symbols) and calculated (black line) Bouguer gravity profile; bottom plot shows the density model. Labeled shot locations from refraction line are indicated as stars across the top of the model. The density structural model demonstrates that the gravitational low beneath the Foreland Belt (model distances 220-300 km) is not caused solely by the deepening of the Moho in that region but requires lower densities in the overlying lower and middle crust than those at equivalent depths on either side. This is evident specifically beneath the upthrust higher density block corresponding with the higher velocities in the upper crust beneath shot locations 2104 and 2105. The structure of higher density regions in the lower Chapter 5. POTENTIAL FIELD DATA MODELING 100 crust corresponds well with the regions of higher velocities in the velocity structural model for Line 21. In particular, the west-facing ramp of higher density material, also resolved from seismic refraction results and by the GRAV3D inversion, successfully reproduces the increase in Bouguer gravity values at the eastern end of the line. At the western end of the line, a narrow low density wedge of material reproduces the strong negative anomaly in the same location as the wedge of low velocities in the velocity model between shot locations 2108 and 2109. The results shown in Fig. 5.11 signify that for reasonable P-wave velocity to density conversions, a density model that successfully reproduces the observed Bouguer gravity profile can be constructed for Line 21. 5.2 Isostatic Balance Considerations A first-order comparison between the observed gravity (Fig. 5.2) and the topography (Fig. 5.3a) in the study area reveals that areas with higher topography generally have strongly negative Bouguer gravity values while lower areas have higher Bouguer values. This trend arises mostly from isostatic compensation of the higher topography with a deepening of the base of the crust. The crustal roots replace denser upper mantle material, which results in a mass deficit for the region and thus a negative Bouguer gravity value. The 2.5-d density model determined for Line 21 successfully reproduces the observed Bouguer gravity values while respecting the constraints provided by the seismic refraction and 3-d gravity inversion results. To provide further insight for the interpretation of both the velocity and density models, the density model was examined to determine whether it was in isostatic equilibrium, a state that exists when the pressure exerted at depth by the Earth's crust and upper mantle across a laterally varying density distribution are all compensated lithostatically to the same degree. This compensation takes place at Chapter 5. POTENTIAL FIELD DATA MODELING 101 a depth below which there exists no lateral variations in density. For the 2.5-d model derived using S A K I , that depth exists at 70 km. To determine the degree of isostatic compensation, the state of stress at the depth of compensation can be calculated. For a column of rock with constant density which experiences no lateral stresses, the lithostatic stress or pressure at a given depth can be calculated (Turcotte and Schubert, 1982) according to crzz = pgz (5.1) where o~zz is the normal stress acting at the base of the column due to the overlying rock, p is the column's density, g is the gravitational acceleration (9.81 m/s2) and z is the depth below the surface. The weight of the column at the depth of compensation can be calculated by simply multiplying both sides of Eq. 5.1 by the cross sectional area of the column, SA. The mass of the column can then be obtained by simply dividing the weight by g. For a column of rock where density varies with depth, isostatic compensation calculations become slightly more complicated as the thickness of particular density layers must be determined and the weight of all layers must be summed. To determine the state of isostatic compensation for the density model for Line 21, the technique outlined in the Appendix of Stacey (1973) was employed. Firstly, the model is subdivided into vertical columns which are then divided into three main blocks (Fig. 5.12): an upper mantle block, a crustal block and a block for parts of the model with densities less than the Bouguer reduction density of 2.67 g/cm3. Next the depths of the five main boundaries, z\ to z 5 are determined; see, Fig. 5.12 for definitions. Finally, the equation {zi + z2)pb + (z4 - z2)pc•+ (z5 - z4)pm = k (5.2) Chapter 5. POTENTIAL FIELD DATA MODELING 102 is solved, where k is the mass per unit area in kilograms per square meter at the base of each column. -1 Z4| Pc •m Figure 5.12: Model for isostatic balance calculations. pb, pc and pm represent the main blocks making up each column. zx is the height of topography, z 2 is the depth below which all densities are greater than or equal to the Bouguer reducing density of 2.67 g/cm3, zz is the depth to the base of the crust if zx = z0, z 4 is the depth to the base of the crust including the crustal root, and z 5 is the maximum depth for lateral changes in density. A n algorithm was constructed to divide the density model derived from S A K I into columns and to calculate the mass per unit area of each column at the depth of compen-sation. Two important assumptions are made in this algorithm; firstly, that the lateral extent of the modeled density bodies beyond the width of the column does not affect the isostatic balance calculation and secondly that thermal isostasy effects in the region are negligible. Results from application of the algorithm are shown in Fig. 5.13. Fig. 5.13b shows what densities would be required in the crust (thin dashed line), compared to the average crustal density for each column (thick solid line), to achieve isostatic balance if Chapter 5. POTENTIAL FIELD DATA MODELING 103 the thickness of the blocks remained constant and if all adjustments had to be taken up by the crust. Fig. 5.13c shows what densities would be required in the upper mantle (thin dashed line), compared to the average upper mantle density for each column (thick solid line), to achieve isostatic balance if again the thickness of the blocks remained constant and if all adjustments had to be taken up by the upper mantle. The results of the isostatic balance calculations are quite revealing. In the portion of the model where the topography increases and the Moho deepens (between model distances of 175 km and 300 km), a mass deficit, which caused the low Bouguer gravity values in the region, exists. Directly to the east of this mass deficit, a mass excess exists where the Moho depth decreased and the higher density upper mantle material took its place. At the easternmost end of the line, the mass deficit increases as the line crosses onto the edge of the ancestral North America where Moho depths start increasing once again towards the ancient craton. Based on Fig. 5.13b and Fig. 5.13c, restoring isostatic balance along Line 21 requires a general decrease in densities along most of the line. These results indicate one of several possibilities (partially taken from Woollard (1969)): 1) Line 21 is not in isostatic equi-librium; 2) thermal isostasy may play an important role in achieving equilibrium; 3) the lateral extent of the density bodies does affect the isostatic calculations; or 4) the base of the compensated crust does not correspond with the depth of the seismically determined Moho. Based on available information, all four possibilities are equally feasible. Chapter 5. POTENTIAL FIELD DATA MODELING 104 w 2.06x10' cd S 3 3 b) 2.9 S 2 8 CP - 2.7 « 2.6 3 o QfJ 3.4 c) 1 3.3 CO C CD 3.2 Q jU 3.1 C (0 S 200 Distance along the line (km) 400 TfT T F - ' i . ' ' ii 11 1 l| L r - ' • ' l I ' L ._ J f t , I | 1 II r I - , f _ l r 1 I ,1 If j K 200 Distance along the line (km) 400 - 1 " 1 . 1 1 i r 1 — 1 ' r. F ~ II h i i i ; : •I1 11 II -1 r' 1 r i i i - r f .'. ^ 1 i i " i ' -- 1 . i ; |i J \ - J ~ 1 " • 1 1 1 i 'I'I ' i 1 Z !! » ' i n r - - J l1 r 1 1 | 1 l l 1 ^ - > 1 • J - I 1 i i 1 r 200 Distance along the line (km) 400 Figure 5.13: Results from calculations of isostatic compensation; a) the mass per unit area in kilograms per square meter of each 1 km wide column of the density model calculated at the base of each column of the density model; b) average density of the model's crust (thick solid line) against the model densities that would be required in the crust to achieve isostatic balance (thin dashed line); c) the average density of the model's mantle (thick solid line) against the model densities that would be required in the mantle to achieve isostatic balance (thin dashed line). Chapter 5. POTENTIAL FIELD DATA MODELING 105 5.3 Aeromagnetic Data Interpretation The aeromagnetic data (Fig. 5.14) for the Line 21 region were acquired in digital format from the Geophysical Data Centre (Geological Survey of Canada, Ottawa) and were generated from a 400 m grid of the profile data. The data set represents a combination of data from both digital and analog surveys (digitized from contour maps) flown in the area and, as such, contains mixed data qualities and flight line spacings ranging between 800 m and 6400 m. White areas of the map correspond to regions where no data have been collected or where the data are not yet available to the public (this is the case for 4 surveys collected north of 60°N in the Line 21 region; they will only be available in May 2000). Examination of Fig. 5.14 reveals a pattern of large long wavelength anomalies east of the Tintina Fault, which are particularly strong in the northeastern corner of British Columbia and correspond with the Nahanni and Fort Simpson terranes which underlie the thick sediments of the W C S B (Cook et al., 1992). West of the Tintina Fault, the map shows more chaotic anomalies with shorter wavelengths. The Tintina Fault is expressed as a very long, quasi-continuous magnetic feature. The magnetic anomalies measured at the surface of the Earth are directly related to the proportions of magnetic minerals (e.g. magnetite) in the rocks of the subsurface. These magnetic minerals are most often found in magmatic rocks which have cooled from a melt. During cooling, the magnetic minerals align themselves with the Earth's mag-netic field and retain that magnetization as long as they are not heated above their Curie temperature (about 575°Celsius for magnetite). With increasing depth in the Earth's crust, both temperature and pressure increase. These trends result in the loss of mag-netization of deeper crustal rocks. For this reason, magnetic anomalies measured at the Earth's surface rarely contain any contributions from the middle crust and below. Given Chapter 5. POTENTIAL FIELD DATA MODELING 106 -760 -200 -100 -50 0 60 120 180 250 350 1700 Figure 5.14: Aeromagnetic map of region surrounding Line 21. The refraction line as well as all shot points are overlain on the map for reference. White regions of the map represent areas where either data has not been collected or where the data is still proprietary. The dark grey box indicates the region over which the reduction to pole transformation was applied. the lithospheric focus of this thesis, I considered that results from aeromagnetic modeling would not contribute significantly to the conclusions and so only a visual interpretation was attempted. 5.3.1 Reduction to pole For gravity surveys, an anomaly tends to be located directly over the mass that is causing the anomaly. However, this is not always the case for magnetic surveys. For regions where the magnetization and the magnetic field are not vertical, magnetic anomalies can often appear to be shifted laterally or distorted due to phase contributions from the Chapter 5. POTENTIAL FIELD DATA MODELING 107 non-vertical dipole (Blakely, 1996). To aid in the interpretation of aeromagnetic data, a transformation can be applied to the data to correct for this distortion and realign the anomalies with their respective sources. This transformation is known as the reduction to pole as it readjusts the data to what they would look like if measured at the Earth's pole where the magnetization and ambient field are both vertical. Subroutine B.26 in the appendix of Blakely (1996) provides an efficient technique for reducing aeromagnetic data to pole. The only required inputs for the routine are the original inclinations and declinations of the magnetization and of the ambient field. The ambient field in the region of Line 21 has an inclination of 78° and a declination of 30°. Since the inclination and declination of the magnetization in the region were unknown, these were assigned the same values as the ambient field. The results of the reduction to pole are shown in Fig. 5.15. Once the transformation was applied, both the original and the reduced to pole magnetic data were gridded using a continuous curvature surface gridding algorithm. As such, regions with no aeromagnetic data coverage (correspon-ding to white regions in the map in Fig. 5.14) appear filled in Fig. 5.15. Anomalies in these filled areas are artifacts from the gridding algorithm and, as such, should not be interpreted as real features. Results from the reduction to pole transformation, shown in Fig. 5.15, indicate that although the anomalies were not displaced laterally due to their inclination and declina-tion, their signatures appeared smeared. The transformation thus succeeded in refining the aeromagnetic anomalies enabling a more accurate qualitative interpretation. 5.3.2 Interpretation of Magnetic Anomalies The long and fairly continuous positive magnetic anomaly associated with the Tintina Fault appears as a striking feature in both the original and the reduced-to-pole maps. Chapter 5. POTENTIAL FIELD DATA MODELING 108 Original Data Reduced to pole Figure 5.15: Comparison of aeromagnetic data before and after reduction to pole. The refrac-tion line has been overlain on both plots for reference. This positive anomaly may be the result of realignment of magnetic minerals from sedi-ments filling in the fault scar due to metamorphic remelting. At the Yukon/British Co-lumbia border where the fault has no surface expression, an enlargement of the reduced-to-pole aeromagnetic map provides one possible means of indicating its exact location (Fig. 5.16). The map reveals that the inferred location of the Tintina Fault lies ~ 10-15 km east of the magnetic high associated with that same fault. This aeromagnetic high also coincides approximately with the location of significantly low velocities at the surface of the velocity structural model between shots 2108 and 2109. These lines of evidence suggest that the inferred location for the Tintina Fault is inaccurate and it should be relocated 10-20 km to the west. Chapter 5. POTENTIAL FIELD DATA MODELING 109 -760 -200 -100 -50 0 60 120 180 250 350 1700 Figure 5.16: Close-up of reduced to pole aeromagnetic map at the Yukon/BC border. White portion of the map corresponds to region with no aeromagnetic coverage. Both the refraction line and the relevant shots points are overlain for reference. The inferred location of the Tintina Fault is indicated with a thin black line. This study indicates that the fault should be located about 20 km west of its inferred location. To the east of the Tintina Fault, large long wavelength magnetic anomalies dominate the aeromagnetic map while anomalies with shorter wavelengths are muted out due to the significant thickness of sediments blanketing the eastern regions (Fig. 5.15). The two main positive anomalies in this area correspond to the Nahanni and the Fort Simpson terranes which have been interpreted as island arcs due to their high content of magnetic minerals. Lacking significant sedimentary cover, the shallower magnetic structures of the accreted Cordilleran terranes result in a more chaotic distribution of short wavelength anomalies west of the Tintina Fault (Cook et al., 1995). Chapter 6 DISCUSSION 6.1 Comparison with other Geophysical Studies in the Region This section compares the results from SNoRE Line 21 with those from other studies conducted in the surrounding region and in the southern Cordillera to assess similarities and thus improve the overall interpretation. While survey results from other geophysi-cal techniques (e.g. seismic reflection and magnetotellurics) along Line 21 are not yet available, nearby studies help to determine the lateral continuity of structures imaged in this thesis. More remote studies provide clues to the large-scale significance of the thesis results. 6.1.1 Comparison with Refraction Studies in the Northern Cordillera Seismic refraction surveying in the northern Cordillera has increased dramatically over the past 5 years due primarily to the efforts of the L I T H O P R O B E project. Although the other results from the most recent refraction experiment (e.g. SNoRE '97) have not yet been published, preliminary results are available for comparison. Early experiments Prior to the SNoRE '97 experiment, only one other refraction experiment had been undertaken near the British Columbia/Yukon border. This study, conducted by Tatel and Tuve (1956), was recorded during the summer of 1955 and consisted of 9 shots 110 Chapter 6. DISCUSSION 111 detonated in the College Fjord of Alaska and another 13 shots detonated near Skagway, Alaska. The Skagway shots were recorded by instruments deployed along the B C / Y u k o n border eastward to the western end of Line 21 of the S N o R E experiment. Despite the coarse receiver spacing and the lack of shot reversals, Tatel and Tuve were able to draw some conclusions about large scale crustal structure for the northern BC/sou the rn Yukon area. Their results revealed that the crust in the region is relatively uniform with an average upper crustal P-wave velocity of 6.1 km/s . A lower crustal reflector between 31 and 34 k m which is not apparent i n the Line 21 data set, is also mentioned i n their findings. Whi le Moho reflections were detected from their survey, no depth estimates were suggested and they did not observe a significant increase in velocity with depth in the lower crust. The P R A S E 85 Experiment Southeast of Line 21 in northeastern Br i t i sh Columbia and northwestern Albe r t a be-tween 54°N and 59°N and 122°W and 115°W, the Peace River A r c h Seismic Experiment ( P R A S E ) took place during the summer of 1985. This experiment, which consisted of four perpendicular refraction/wide-angle reflection lines, was designed to determine the velocity structure of the crust and upper mantle of the region surrounding the Peace River A r c h , a regional geological structure with anomalous vertical displacements compared to the rest of the W C S B (Zelt and El l is , 1989). Whi le the P R A S E experiment's proximity to Line 21 provides some hope that di-rect comparison between its results and those from Line 21 could be insightful, the experiment's location southeast of the Great Slave Lake Shear Zone sheds doubt on the significance of any crustal similarities given that the G S L S Z juxtaposes different base-ment terranes (see F i g . 1.3). Nonetheless, the profiles provide good comparisons between Chapter 6. DISCUSSION 112 large-scale surficial features and Moho structure which are less affected by the basement terranes. The western end of line B of the P R A S E refraction experiment lies along the same geological strike as the eastern end of Line 21 and trends roughly E-W, from just east of the Cordilleran deformation front, oyer the axis of the Peace River Arch. Despite the differing basement structures between the two fines, I note that the western portion of fine B resolves a thickening from 1.95 km to 4.90 km of low velocity sediments toward the west, similar to the westward thickening low velocity package imaged at the eastern end of Line 21. In terms of Moho depth and character, fine B of the P R A S E experiment displays a diffuse transitional zone for the crust-mantle boundary at a depth of 40 km at the western end of the fine which deepens to 43.5 km and becomes a sharper discontinuity to the east. The western transitional zone coincides with both the depth and character of the lower crustal transition zone modeled at the eastern end of Line 21, indicating that the zone modeled along Line 21 could extend along strike of the entire ancient cratonic margin. The SNoRE '97 Experiment As mentioned in Chapter 2, the seismic refraction data from Line 21 represents only a quarter of the data collected during the SNoRE '97 experiment. Although the work on the other three fines is ongoing, comparisons with preliminary results can be made. Line 11 (Fig. 1.1) which lies to the northeast of Line 21, extends 600 km from east of Yellowknife to west of Nahanni Butte in the Northwest Territories and crosses, from east to west, the Archean Slave craton, the Great Bear magmatic arc, the Hottah terrane, the Fort Simpson terrane and the Nahanni terrane. As such, Line l l ' s westernmost extent Chapter 6. DISCUSSION 113 crosses similar basement domains as the easternmost end of Line 21, allowing for an along-strike offset of ~ 200 km. Although research is ongoing on Line 11, recent publications (Fernandez Viejo et al., 1999a, 1999b) shed some light on the velocity and Moho structures in the region. To compare the portions of the velocity models that he across the same basement domains, Fig. 6.1 shows the westernmost section of the Line 11 model and the easternmost edge of the model derived in this thesis. In the upper crust, both models resolve a west facing ramp of higher velocity material overlain by a westward thickening package of low velocities. Fernandez Viejo et al. (1999a) correlate this westward thickening package with a deep Proterozoic sedimentary basin also imaged with multichannel reflection techniques (Cook et.al., 1999). On the whole, resolved Moho depths at the western end of Line 11 are much shallower (~ 30 km) than those resolved for Line 21 (~ 37 km). Normal-moveout-corrected and stacked PmP reflections along this portion of Line 11 indicate that the Moho is relatively diffuse beneath the Nahanni and Fort Simpson terranes (Fernandez Viejo et al., 1999b). This diffusive character may explain why the high velocity lower crustal transition zone that masked Moho reflections along Line 21 is not resolved in the velocity model for Line 11. However, it is also possible that the transition zone, also resolved to the south in the P R A S E experiment, does not extend to the north. In terms of upper mantle velocities, Line 11 results show a low velocity (~ 7.7 km/s) upper mantle beneath the Fort Simpson/Hottah transition zone which is attributed to collision and subsequent delamination of the Fort Simpson terrane's lower crust beneath the Hottah (Fernandez Viejo et al., 1999a). Since the eastern extent of Line 21 does not extend to this same transition zone, the low velocity upper mantle resolved beneath the Nahanni terrane portion of Line 21 probably has a different origin. Figure 6.1: Comparison of velocity models from Line 11 and Line 21; a) location map with relevant portion of Line 11 shown in green and relevant portion of Line 21 shown in red; b) velocity model for relevant portion of Line 11; c) velocity model for relevant portion of Line 21. Shot locations are numbered and indicated with stars along the surface of both models. Model regions with no ray coverage have been blocked out. Relevant terranes are labeled along the top of both model plots. Chapter 6. DISCUSSION 115 Line 22, which lies to the west of Line 21 and intersects Line 21 at shot 2108, extends 460 km south from Watson Lake, Yukon, to Stewart, British Columbia. Hammer et al. (1999) have modeled the results from this profile, which runs south from the Omineca Belt through the Intermontane Belt to the Coast Belt. The intersection of Line 22 with Line 21 at shot 2108 (Fig. 1.1) provides a unique opportunity for direct comparison of models constructed from independent data sets. The preliminary velocity structural model derived for Line 22 (Hammer et al., 1999) was meshed directly with the final velocity model from this thesis (Fig. 6.2). Since shots 2108 and 2208 were detonated at the same location, the models were meshed together at the location of that cross-over point. Figure 6.2 reveals that the models for Line 21 and Line 22 are well matched. Both models resolve extremely slow near-surface velocities (~ 2 km/s) in the location of and either side of their cross over point at shot 2108 (2208). In addition, the thickening of the lower crust at the western end of Line 21 near the location of the Tintina Fault coincides with the thickened lower crust resolved at the northern end of Line 22. Moho depths across Fig. 6.2 reveal a good match with a slight shallowing of the Moho beneath the Tintina Fault. This Moho expression beneath the fault lends credence to the idea that the Tintina extends down through the entire crust and affects the crust-mantle boundary. In terms of upper mantle velocities, the models agree at their cross-over point and an overall decrease in upper mantle velocities can be observed from the Foreland Belt westward. Given that the resolved models for lines 21 and 22 are not perpendicular to geological strike, Fig. 6.2 gives a false impression of the perpendicular to strike width of the imaged features. To obtain a better idea of the actual perpendicular to strike widths for the different tectonic blocks, both models were projected onto a line that ran perpendicular to geological strike between shot location 2101 (just south of Fort Nelson) and Stewart, Chapter 6. DISCUSSION 116 Br i t i sh Columbia. The projected model (Fig . 6.3) indicates that the west-facing ramp of higher velocities at the eastern end of Line 21 and the corresponding westward thickening package of low velocities had widths that were confined to less than 100 km. Similarly, the upthrust higher velocity middle crustal material at the eastern front of the Foreland Belt affected less than 100 k m of the width of the Cordilleran Orogen. Line 31 (Fig. 1.1) which lies to the northwest of Line 21 and extends 300 k m N E to S W in the central Yukon, crosses, from east to west, the deformed miogeoclinal rocks of the Omineca Belt , the Tint ina Fault, the displaced continental margin and the ac-creted terranes of the Intermontane Belt . Despite the Tin t ina Fault dividing different morphogeological belts than it does along Line 21, comparisons can be made between the resolved velocities of the Omineca Belt and the significance of the Tin t ina Fault along both lines. Preliminary results from Creaser and Spence (1999) indicate a general increase in up-per and middle crustal velocities westward across the Tint ina Fault and a corresponding increase i n Moho depths from 32 k m to 37 km. The low upper crustal velocities seen near the surface of the Omineca Belt along Line 31 agree with the overall low velocities at the surface at the western edge of Line 21. However, the modeled increase in Moho depths across the Tint ina Fault along Line 31 cannot be resolved along Line 21 due to the fault's location at the westernmost end of the line where P m P coverage diminishes. During the S N o R E '97 experiment, broadside deployment of instruments for both lines 21 and 22 was undertaken to address the question of crustal thickness on either side of the T in t ina Fault. The analysis and modeling of these broadside data wi l l begin once the inline profile analyses for lines 21 and 22 are complete. aptei6. DISCUSSION (iu>l) qjdsa Chapter 6. DISCUSSION 118 b) - 1 3 0 - 1 2 5 ' -120" a) - 1 3 0 125" INTERMONTANE OMINECA FORELAND ANCESTRAL N.A. I ~ Stikinia |CC CA Cassia SM j FS 100 200 300 Distance (km) 400 500 P-wave Velocities (km/s) 1.9 3.0 4.0 5.8 6.0 6.2 6.4 6.8 7.0 7.5 8.1 8.4 8.9 Figure 6.3: Velocity models for lines 21 and 22 projected perpendicular to geological strike; a) location map with the red line indicating the line onto which both velocity models were projected; b) projected velocity models for lines 21 and 22. The Moho is indicated by the thick black line underlain by a thicker white line. Across the top of the plot, both the corresponding terranes and morphogeological belts are indicated. The abbreviations SM, CA, CC and FS refer to the Slide Mountain, Cassiar, Cache Creek and Fort Simpson terranes respectively. Shot point locations are indicated with yellow stars for Line 21 and with white stars for Line 22. Chapter 6. DISCUSSION 119 6.1.2 Comparison with Refraction Studies in the Southern Cordillera Several extensive seismic refraction surveys have been conducted in the southern Cor-dillera of British Columbia as part of the L I T H O P R O B E project and as part of the earlier Project Edzoe. As some of these studies transect the same morphogeological belts as those in the northern Cordillera, comparison of survey results along geological strike with the Line 21 results can provide information related to continuity of crustal and mantle structures along the entire Canadian Cordillera. Early experiments A n early refraction study (Chandra and Cumming, 1972) was collected during the sum-mers of 1962-1969 from Greenbush Lake, British Columbia, to Swift Current, Saskatche-wan. The 700 km long refraction profile extended from the Omineca Belt of the southern Cordillera across onto the edge of ancestral North America and provides perhaps the best comparison in terms of location along geological strike for the results from Line 21. Unfortunately, since the shot spacing was fairly coarse and not all shots were reversed, only large scale comparisons can be made reliably. The coarse refraction model constructed by Chandra and Cumming (1972) was achie-ved using the method of delay time solutions and was subsequently adjusted based on coincident reflection results and gravity modeling. The portion of their model relevant to this study displays low upper crustal velocities in the Foreland Belt which increase slightly onto the edge of cratonic North America. The lower crust also exhibits a slight increase in velocities eastward. A deepening of the Moho to ~ 45 km beneath the Fore-land Belt reveals the location of a 5-km-thick crustal root beneath the Rocky Mountains which is accompanied by upper mantle velocities of 8.20 km/s. These upper mantle ve-locities increase to 8.33 km/s beneath ancestral North America. While the Moho depth Chapter 6. DISCUSSION 120 and upper mantle velocities resolved by Chandra and Cumming (1972) generally agree with the more recent results in the southern Cordillera (Zelt and White, 1995) and with coarse features from the velocity model for Line 21, the upper crustal model features from Chandra and Cumming (1972) are too coarse to provide significant insight for this thesis study. A subsequent refraction profile for Project Edzoe was collected during the summer of 1969 (Mereu et al., 1977) from Greenbush Lake, British Columbia, to Fort McMurray, Alberta. This 800 km long profile extended northeast from the Rocky Mountain Trench, across the Foreland Belt to the Alberta Plains of ancestral North America. As with the Chandra and Cumming (1972) study, shot spacing was fairly coarse and not all shots were reversed so only large scale comparisons can be made reliably. Pn arrivals from the data set indicate apparent upper mantle velocities ranging between 8.5 km/s and 8.9 km/s, which are much higher than those modeled beneath the Foreland Belt for Line 21. These anomalously high upper mantle velocities are interpreted by Mereu et al. (1977) as arising from a high velocity upper mantle layer existing between 55 km and 70 km depth. Such results support the existence of high velocity material in the upper mantle at the base of the velocity model for Line 21 (Fig. 4.10). The SCoRE '90 Experiment The most recent experiment to image the Omineca and Foreland belts of the southern Cordillera took place in 1990 as part of the L I T H O P R O B E Southern Cordillera Refrac-tion Experiment (SCoRE). Data from Line 9 of the experiment were modeled by Zelt and White (1995) to provide a velocity structural model for the Omineca and Foreland belts of the southern Cordillera. While Line 9 did not extend east of the Cordilleran deformation front, direct comparisons with the results from the western portion of Line Chapter 6. DISCUSSION 121 21 can be made since the receiver spacings and modeling techniques used by Zelt and White (1995) are similar to those used in this thesis. From the Omineca Belt eastward, the velocity model for Line 9 resolves a deepening of the Moho from 34 km to 41 km, this increase in crustal thickness being accompanied by a general decrease in lower crustal velocities from 6.6-6.7 km/s to 6.3-6.4 km/s in the Foreland Belt. Upper mantle velocities increase from 7.9 km/s to 8.0 km/s from the Omineca Belt eastward. While the velocity model for Line 21 also resolves an increase in Moho depths beneath the Foreland Belt and a corresponding decrease in lower crustal velocities, the absolute changes are not as large. Conversely, upper mantle velocities for Line 21 increase from 7.9 km/s to 8.3 km/s beneath the Foreland Belt, a more dramatic increase than in the southern Cordillera. The variations along geological strike of Moho depths beneath the Foreland Belt, as resolved from lines 9 and 21, indicate an overall eastward thickening of the crust toward cratonic North America; more significant thickening is evidenced in the southern Cordillera where the topography is more elevated. The higher upper mantle velocities beneath the Foreland Belt in the northern Cordillera suggest a colder upper mantle than that which exists in the south. However, such an interpretation disagrees with recent high heat flow values collected in the northern Cordilleran Foreland Belt (see subsection 6.1.4). Clowes et al. (1995), following an earlier suggestion by Gough (1986), propose a thin lithosphere (~ 55-60 km) for the southern Cordillera below the western Omineca and Intermontane belts, thickening considerably to the east. The thin lithosphere is attributed to upwelling of warm mantle below the Omineca Belt. There is, however, no clear evidence of a similar mechanism acting in the northern Cordillera. Chapter 6. DISCUSSION 122 6.1.3 Comparison with Reflection Studies in the Northern Cordillera Seismic reflection profiles are being collected during the winter of 1999 along lines 21, 22 and 31 of the SNorCLE transect. These will eventually provide further constraints on the interpretations put forth in this thesis. Unfortunately, the only other reflection survey presently available for comparison in the northern Cordillera exists northward along strike from the eastern portion of Line 21. Cook et al. (1999) present the results from a seismic reflection line that was collected in 1996 along Line 11 of the SNorCLE transect. Once again, since the western edge of Line 11 is offset along strike from the eastern end of Line 21, some general comparisons can be attempted. Figure 6.4 shows a direct comparison between the interpreted reflection results and the velocity structural model resolved in this thesis. The relevant portions of lines 11 and 21 are shown on the location map in Fig. 6.1a. The westernmost portion of the Line 11 reflection profile is dominated by west-dipping reflections overlying two more prominent west-dipping reflections ( R l and R2) that flatten near the base of the crust to the west (Fig. 6.4). The latter are interpreted as the western edge of the Fort Simpson terrane while the overlying reflections define the Fort Simpson Basin (FSB) of Proterozoic metasediments that developed as a result of lithospheric extension following the collision of the Fort Simpson and the Hottah terranes at 1.84 Ga (Cook et al.,1999). The thinned crust of the Fort Simpson terrane (FST) lies below R l . The west-facing ramp overlain by a 10-15 km package of Proterozoic sediments imaged from the Line 11 reflection results coincides with the west-facing ramp of high velocities overlain by westward thickening low velocity material resolved along Line 21. Similarly, the shallowing of the Moho along Line 11 which is interpreted as being due to crustal extension during the Proterozoic, agrees with the shallowing of the Moho towards the Foreland Belt resolved on Line 21. Chapter 6. DISCUSSION 123 from Cook et al., 1999 Figure 6.4: Comparison between interpreted reflection results from Line 11 and velocity model for Line 21. Shot locations for Line 21 are numbered and indicated with stars along the surface of the top plot. Regions with no ray coverage along Line 21 are blocked out from both plots. The Moho is indicated by the thick black line underlain by a thicker white line. Relevant terranes are labeled along the top of both plots. On the plot of the reflection data, Rl and R2 denote two prominent west-dipping reflections. Abbreviations FSB and FST refer to the Fort Simpson Basin and the Fort Simpson terrane respectively. South of Line 21, a seismic survey along the Alaska Highway recorded by industry in 1991 samples the Fort Simpson terrane. The results from that study (Cook and Van der Velden, 1993) resolve the same west-facing ramped structure overlain by a considerable thickness of sediments. Also, as mentioned in Chapter 1, two other deep reflection lines recorded further north in the Northwest Territories (Cook et al., 1991; Clark and Cook, 1992) across the western margin of ancestral North America have revealed a regionally extensive west-facing ramp at depth which has been confirmed by borehole data (Fig. 1.8). Chapter 6. DISCUSSION 124 The ramp imaged in all these surveys (illustrated in Fig. 1.8) has been interpreted as the rifted Proterozoic margin of ancestral North America which was later overlain by passive margin sediments. The wide distribution of this ramp along the cratonic edge suggests that it is a feature that exists along the entire western margin of ancestral North America. 6.1.4 Comparison with Heat Flow Studies Extensive heat flow studies conducted in the southern Cordillera (Hyndman and Lewis, 1995; Hyndman and Lewis, 1999) and recent studies in the northern Cordillera (Lewis and Hyndman, 1998; Lewis and Hyndman, 1999) have resolved a westward increase in heat flow values from cratonic North America to the Cordillera. For most of the craton, heat flow values range between 40-60 m W / m 2 while in the Cordillera they range between 80-100 m W / m 2 (Hyndman and Lewis, 1999). Although the northern Cordilleran heat flow database is far from complete, differences in heat flow character along the strike of the Cordillera-craton transition are readily apparent. In particular, the transition from cold craton to hot Cordillera in the southern Cordillera exists 100 km west of the Cordilleran deformation front at the location of the Rocky Mountain Trench, but the transition in the northern Cordillera appears to exist considerably east of the deformation front. On the craton just east of Fort Providence, N W T , below which lies the eastern part of the Hottah terrane, the heat flow is ~ 100 m W / m 2 ; west of Yellowknife, N W T , situated on the Slave Craton, the heat flow is 53 m W / m 2 (Lewis and Hyndman, 1998). The cause of the unusually high heat flow near Fort Providence is a subject of current research. Recent heat flow measurements collected from hydrocarbon exploration wells, west of Fort Providence, N W T and within the Hottah, Fort Simpson and Nahanni terranes, Chapter 6. DISCUSSION 125 have recorded heat flow values of 89 ± 17 m W / m 2 across the terranes and up to 118-136 m W / m 2 in the Nahanni terrane itself (Lewis and Hyndman, 1999). These values are significantly higher than cratonic heat flow measurements collected in the southern Cordillera, suggesting that the thermal regime varies significantly along the length of the Cordillera-craton transition. The average heat flow for British Columbia and the Yukon north of 59°N is 103.3 ± 21.0 m W / m 2 (Lewis and Hyndman, 1999), an extremely high value. A recent study in Australia, where similar high heat flow values have been measured within a craton, concludes that the source of the extra heat is anomalously high heat production within the crust (Sandiford and Hand, 1998). Recent work in the Foreland Belt of the SNorCLE Transect suggests a similar origin for the high heat flow values in the region of Line 21 (Hyndman, 1999; personal communication). Using radioactive heat generation measurements from petroleum well basement samples, heat generation within the Foreland Belt has been estimated at 4.5 m W / m 3 . This value suggests that 45 m W / m 2 of the heat flow in the region originates within the upper 10 km of the crust. These results offer some explanation as to how the upper mantle beneath the Foreland Belt could have remained cold (inferred from the high Pn velocities) in an area of overall high heat flow. The location of the high heat flow values in the Nahanni terrane corresponds with a region of low velocities and low densities resolved within the upper mantle along Line 21. In this location, a significant portion of the crustal thickness (~ 20 km) is composed of metasediments. While it is possible that the low velocities and densities in the up-per mantle were emplaced beneath the edge of the craton by some uncertain tectonic mechanism, the overlying significant thickness of sediments suggests a crustal source. If the amount of radiogenic potassium in the overlying sediments was particularly high, the resulting high heat production in the crust could have heated the upper mantle from Chapter 6. DISCUSSION 126 above, yielding the low velocities and densities resolved in this thesis. Based on the prominent heat flow transition in the southern Cordillera, which occurs 100 km west of the Cordilleran deformation front, the Foreland Belt east of the Rocky Mountain Trench and the edge of cratonic North America are interpreted as comprising a strong lithosphere to depths of 100 km (Hyndman and Lewis, 1999). One significant result of such a scenario is that deformation from Late Mesozoic and Early Tertiary tectonics was limited to the overlying sedimentary thrust sheets and did not affect the stronger lower crust and upper mantle. However, in the northern Cordillera, where the thermal transition is more gradual, the higher heat flow values extending onto the craton suggest a weaker lithosphere overall. To the north of Line 21, this weakened lithosphere may have allowed Cordilleran deformation to extend eastward onto the craton, causing the bulge in the Cordilleran deformation front (Hyndman and Lewis, 1999; see Fig. 1.2). 6.2 Final Velocity Model within the Context of the Regional Tectonics 6.2.1 Tectonic Implications of Upper and Middle Crustal Results V The edge of cratonic North America, which lies at the southeasternmost part of Line 21, is characterized by a region of relatively low velocities which thicken westward (from 5 km to 15 km) toward the Foreland Belt over a west-facing ramp of higher velocities in the middle crust. These crustal features are interpreted as the ramped edge of cratonic North America overlain by the passive margin sediments that developed following the Proterozoic rifting events which opened the proto-Pacific Ocean. These features and the accompanying interpretation agree with results from refraction and reflection studies (mentioned above) that have modeled the same west-facing ramp as a continuous feature along the edge of cratonic North America. The eastern front of the Foreland Belt, which is characterized by anomalously high velocities of ~ 6.4 km/s within 5 km of the surface, Chapter 6. DISCUSSION 127 is interpreted to have developed following upthrust of middle crustal material over the edge of cratonic North America during the cofhsional tectonics that caused the formation of the northern Rocky Mountains. At the western end of Line 21, the steeply-sided trench of low velocity material, that appears to have a velocity expression to the base of the upper crust, is interpreted as sediments and brecciated surface rocks filling in the Tintina Fault scar. The location of the low velocity trench, 15 km to the west of shot location 2108, is coincident with the location of a strong aeromagnetic feature that extends north across the BC/Yukon border from the northern Rocky Mountain Trench. These results suggest that the inferred location of the Tintina Fault to the east of shot location 2108 is in error; it should be relocated ~ 20 km to the west. The middle crustal region of the velocity model for Line 21 contains velocities ranging between 6.0 km/s and 6.5 km/s with variable thicknesses. While its thickness at the edge of cratonic North America is approximately 10 km, the middle crust thickens beneath the eastern front of the Foreland Belt to 15 km, thins to 10 km in the western portion of the Foreland Belt and thickens again significantly in the Omineca Belt to almost 20 km. The thickening that occurs in the Omineca Belt is accompanied by a 10 km offset of the 6.4 km/s contour in the middle crust beneath shot location 2108. This dramatic offset is interpreted as resulting from the juxtaposition of different crustal blocks due to northwestward displacement of the cratonic margin along the Tintina Fault. The Tintina is thus interpreted as a crustal-scale transform fault that extends to the base of the crust. Since the offset middle crust occurs beneath shot location 2108 and the surface exposure of the Tintina is resolved 15 km to the west of shot location 2108, it is possible that either a) both features are associated with the Tintina Fault which must have an eastward dip of ~ 68°, or that b) the Tintina is composed of a system of faults similar to those making up the San Andreas Fault System in California. If the Tintina is indeed a system of Chapter 6. DISCUSSION 128 faults, the middle crustal offset in velocities may have been caused by a different fault than that evidenced at the surface to the west of shot location 2108. 6.2.2 Tectonic Implications of Lower Crustal and Upper Mantle Results The velocities in the lower crust of the model range between 6.5 km/s and 7.3 km/s with high velocity (> 6.8 km/s) regions located beneath the eastern portion of the edge of cratonic North America, at the Cordilleran deformation front and beneath shot location 2107: The high velocity (7.3 km/s) pinched out layer beneath the edge of cratonic North America, whose base (the Moho) produced reflections which were hidden in the coda of the reflections resulting from the overlying R l reflector, is interpreted as being composed of mafic material which may have resulted from underplating of mantle material during the Proterozoic rifting episodes. As North American basement is inferred to extend to the western edge of the Omineca Belt, the high velocity regions in the lower crust at the Cordilleran deformation front and beneath shot location 2107 may have resulted from upwelling of mafic mantle material through the extended and thinned North American basement during the Proterozoic rifting events. It should be noted that due to poor lower crustal resolvability, these two high velocity regions may however be insignificant. The upper mantle along Line 21 is characterized by highly variable velocities which are constrained by excellent Pn ray coverage. At the western end of the line, typical Cordilleran velocities of 7.9-8.0 km/s are resolved beneath the Omineca Belt while ano-malously high velocities of 8.3 km/s are seen beneath the front edge of the Foreland Belt. These high upper mantle velocities, which may represent deeper upwelled mantle mate-rial, are juxtaposed against low velocity (7.7 km/s) upper mantle material beneath the Nahanni terrane of the cratonic edge. The lower upper mantle velocities at the eastern end of the line, while anomalous, correspond with a localized region of low densities in the Chapter 6. DISCUSSION 129 upper mantle resolved from the 3-D Bouguer gravity data inversion, and with the locus of high heat flow values measured along the cratonic edge. From Fig. 1.8, it is interesting to note that the deformation front and the west-facing crustal ramp interweave with each other the entire length of the Canadian Cordillera. In the region of Line 21 and north into the southern Yukon, the deformation front is outboard of the ramp. This geometry is responsible for reducing the amount of deformation imposed on the miogeoclinal sedi-ments in the northern Cordillera around Line 21 and may also have had an impact on lower crustal geometry as the crust would have remained relatively thin in the Nahanni terrane without the added load of the Foreland Belt. These lines of evidence suggest that beneath the edge of cratonic North America in northeastern British Columbia, hot, rela-tively low velocity upper mantle material may have been preferentially emplaced beneath the cratonic edge in a pool confined between the Cordilleran deformation front and the ramped edge of cratonic North America. 6.2.3 Interpretation of Variations along-strike of the Cordillera-Craton Tran-sition Zone As mentioned in Chapter 1, the rocks of the North American miogeocline that make up the Foreland Belt have been thrust over western ancestral North America as far as 50 km in the northern Cordillera and 200 km in the southern Cordillera (Gabrielse et al., 1991). These differing amounts of overthrust and differing lithospheric characters of the northern and southern Cordillera provide a framework for the interpretation of along-strike variations in the Cordillera-Craton transition zone. In the southern Cordillera, the heat flow values at the Cordillera-craton transition are much lower than in the northern Cordillera. Despite the dichotomy that arises due to the lower upper mantle velocities at the transition in the southern Cordillera, the lower Chapter 6. DISCUSSION 130 crust and upper mantle are inferred to be colder and thus, stronger. Since the west-facing ramp of the cratonic edge of North America in the southern Cordillera presently lies west of the Cordilleran deformation front, during the formation of the Foreland Belt of the Cordillera, the miogeoclinal sediments were detached from the strong undeformed crustal basement and upthrust over the ramped cratonic edge. As both the lower crust and upper mantle were strong, the overloading of the edge of the craton had little effect on the upper mantle and only served to thicken the crust in the region (Fig. 6.5). Southern Cordillera / \ WCSB V js North American Craton Blfes*-*-Northern Cordillera - — J W C S B ^ ^ - ^ y ^ North American Craton Figure 6.5: Interpreted Foreland-craton transitions for the southern and northern Cordilleras. The top panel illustrates the interpreted Foreland-craton transition in the southern Cordillera with a thickened crust overlying a strong, undeformed upper mantle. The lower panel illustrates the interpreted transition in the northern Cordillera where the overloading of sediments onto the cratonic margin resulted in depression of the Moho and deformation within the upper mantle. Chapter 6. DISCUSSION 131 Proterozoic Low Velocity Upper Mantle I High Velocity Upper Mantle Middle Jurassic Pass ive— Cratonic North America b - - i - * . Low Ve lo^ i i y - ^ - , ~-- . Upper Mantle High Velocity Upper Mantle Figure 6.6: Interpreted evolution of Foreland-craton transition in the northern Cordillera. Top panel illustrates how the Foreland-craton transition may have appeared at the end of the Proterozoic with passive margin sediments overlying the rifted margin of cratonic North America. The lower panel illustrates the impact of the accretion of the western morphogeological belts on the Foreland-craton transition as the passive margin sediments were upthrust onto and overloaded the cratonic edge to form the Foreland Belt. In the northern Cordillera, the higher heat flow values at the Cordillera-craton tran-sition suggest a hotter and weaker lower crust and upper mantle. As the miogeoclinal se-diments were thrust eastward onto the extended crust of cratonic North America (Fig. 6.5 (lower panel) and Fig. 6.6), overloading of the edge of the craton took place outboard of the west-facing ramped structure. Since the lower crust and upper mantle were weaker in the northern Cordillera, this overloading may have caused the depression in the Moho beneath the Foreland Belt imaged in the velocity model for Line 21. Since for high heat flows, the lower crust and upper mantle are decoupled, this depressed Moho may have displaced the low velocity upper mantle material to either side of the deepened Moho. As this lower velocity upper mantle material was moved away, deeper high velocity man-tle material may have moved up to take its place, resulting in the higher upper mantle Chapter 6. DISCUSSION 132 velocities imaged beneath the Foreland Belt along Line 21. Since the ramp of the cratonic edge lies east of the deformation front, a region of the miogeocline remained unaffected during the overloading of the western edge of cratonic North America. As such, the Moho below that region was not displaced downward. As the low velocity upper mantle material was displaced away from beneath the Foreland Belt, it may have been unable to descend beneath the eastward thickening cratonic edge due to its lower density and as such, became pooled beneath the Nahanni terrane. This pooling could have resulted in the low upper mantle velocities and densities in the region resolved in this thesis and offers a possible source for the high heat flow values recorded by Lewis and Hyndman (1999). 6.2.4 Interpreted Lithospheric Structure for Northeastern British Columbia Using the interpretation of the velocity structural model for Line 21 and the interpreted reflection results from Line 11 (Cook et al., 1999), an interpreted cross-section of the crust and upper mantle (Fig. 6.7) was constructed from the Archean Slave Craton in the east to the western edge of the Cordilleran morphogeological belts of cratonic North American origin. From this cross-section, the evolution of the transition at depth between the eastern portion of the Cordillera and the Proterozoic rifted margin of ancestral North America can be addressed. Following the Wopmay Orogen (described in Chapter 1), the western edge of cratonic North America developed into a west-facing ramped rifted margin which was subsequently overlain by the Proterozoic-Paleozoic sediments of the Fort Simpson Basin. During these rifting events, the cratonic edge of North America was thickened due to underplating of mafic upper mantle material. Chapter 6. DISCUSSION 133 SNorCLE Reflection Line 11 SNoRE Line 21 Tintina Fault O.B. I Foreland Belt Ancestral North America Fort Nahanni : Simpson Terrane i Terrane Proterozoic Mantle North American Lower Lithosphere F i g u r e 6.7: Interpreted cross-section for the Cord i l le ran-Nor th Amer i can craton t ransi t ion. Locat ions of S N o R E Line 21 and S N o r C L E reflection line 11 are shown across the top of the diagram. The lateral extents of cratonic N o r t h A m e r i c a , the Foreland Bel t and the Omineca Bel t ( O . B . ) are shown across the top of the diagram. The inferred dip of the T i n t i n a Faul t is indicated w i t h a thick dashed l ine at the western end of L ine 21. D u r i n g the compress iona l events tha t resul ted f rom the accre t ion of the western C o r -d i l l e r an belts d u r i n g the M i d d l e Jurass ic (ca. 170 M a ) , these passive m a r g i n metasedi -ments a n d the u n d e r l y i n g N o r t h A m e r i c a n m i d d l e crust were up th rus t on to the c ra ton ic m a r g i n to f o r m the Fo re l and B e l t . However , i n contrast to the southern C o r d i l l e r a where the r a m p e d edge of c ra ton ic N o r t h A m e r i c a p layed a direct role i n the up th rus t of passive m a r g i n sediments , the r a m p e d edge of c ra ton ic N o r t h A m e r i c a l ay east of the n o r t h e r n R o c k y M o u n t a i n s of the Fo re l and B e l t and as such, several hundreds of k i lomete rs of the F o r t S i m p s o n B a s i n were left undeformed . T h e u p t h r u s t i n g of passive m a r g i n sediments to f o r m the th i cke r crust of the Fo re l and B e l t resu l ted i n a depression i n the M o h o due to a weaker crust and upper man t l e i n the n o r t h e r n C o r d i l l e r a . T h i s M o h o depression m a y have caused the p a r t i t i o n i n g of the lower ve loc i ty uppe r man t l e m a t e r i a l to ei ther side of the Fo re l and B e l t and the u p w e l l i n g of Chapter 6. DISCUSSION 134 higher velocity upper mantle material directly beneath the Foreland Belt. By the early Eocene, the formation of large scale right-lateral transcurrent faults such as the Tintina Fault, had resulted in the north-westward displacement of the cratonic margin along the Tintina Fault (system(?)) and thus, the juxtaposition of crustal blocks of varying thicknesses that are resolved at the western end of Line 21. 6.3 Conclusions Through extensive 2-D forward and inverse traveltime modeling and 2-D forward amp-litude modeling, a velocity structural model for Line 21 of the SNoRE '97 refraction/wide-angle reflection experiment was constructed from the Western Canada Sedimentary Basin, across the deformation front of the northern Cordillera and into deformed North Ameri-can rocks, to the Tintina Fault, the eastern boundary of the accreted terranes. The upper crust of the resolved velocity model revealed a westward thickening (from 5 km to 15 km) of low velocity material overlying a west-facing ramp of higher velocity lower crustal material at the southeastern end of the line. From comparison with other refraction and reflection surveys along geological strike, this feature is interpreted as the ramped rifted margin of cratonic North America overlain by passive margin sediments. At the western end of the line, a trench of low velocity material between shot locations 2108 and 2109 is inferred to represent sediments and brecciated surface rocks filling in the Tintina Fault. These results suggest a location for the Tintina Fault which differs from the present inferred location which is 20 km to the east. At depth beneath shot 2108, the juxtaposition of differing thicknesses of lower crustal blocks suggests that the fault may dip to the east or be made up of a system of faults which extend down to, at least, the base of the crust. Beneath the Foreland Belt, high upper crustal velocities resolved within 5 km of the Chapter 6. DISCUSSION 135 surface, which are juxtaposed to the east against the westward thickening low upper crustal velocities of the passive margin metasediments, are interpreted as middle crustal material that was upthrust over the edge of the craton during the accretion of the western Cordilleran belts. A deepening of the Moho beneath the Foreland Belt, accompanied by anomalously high upper mantle velocities flanked by low upper mantle velocities, is inferred to have resulted from deformation of a weak lower crust and upper mantle due to overloading of the cratonic edge in a region of high heat flow. Using results from a combination of 3-D Bouguer gravity inversion, 2.5-D Bouguer gravity forward modeling, aeromagnetic comparisons and heat flow analysis, the inter-pretation of the velocity modeling results were constrained and a framework for the tectonic evolution of northeastern British Columbia was suggested within the context of the velocity, density and heat flow constraints. The tectonic evolution of northeastern British Columbia appears to have been di-rectly affected by the geometry of the rifted western margin of cratonic North America. Following the Proterozoic rifting events, a passive margin sequence was overlain on the west-facing ramp of the cratonic edge. These passive margin sediments and underlying middle crustal material were subsequently upthrust over the edge of the craton during the formation of the Foreland Belt. Since the upthrusting of material occurred outboard of the ramped edge of cratonic North America, over 100 km of passive margin sediments were unaffected by the Cordilleran Orogen. As the Moho was depressed in isostatic re-sponse to the thickened crust outboard of the ramped cratonic edge, the crust remained relatively thin compared to the southern Cordillera and a section of thinner crust re-mained unaffected between the Foreland Belt and the thicker craton to the east. Due to either downward heating from high heat production within the overlying crust or from preferential emplacement of hot upper mantle material beneath the Nahanni terrane, Chapter 6. DISCUSSION 136 a localized pocket of low velocities and densities developed beneath the western mar-gin of ancestral North America. The formation of large scale right-lateral transcurrent faults such as the Tintina Fault during the Eocene later resulted in the northwestward displacement of the cratonic margin at the western end of Line 21. References Aitken, J .D. and M . E . McMechan, 1991. Middle Proterozoic Assemblages, in Ge-ology of the Cordilleran Orogen in Canada, edited by H . Gabrielse and C.J . Yorath, Geological Survey of Canada, Geology of Canada No.4, Vol.G-2, p.99-124. Amor, J . , 1996. P L O T S E C : Generalized Software for Seismic Refraction Data, L S P F Newsletter, 23, p.1039-1043. Barton, P.J., 1986. The relationship between seismic velocity and density in the continental crust - a useful constraint?, Geophysical Journal of the Royal As-tronomical Society, Vol.87, p. 195-208. Berry, M . J . , W.R. Jacoby, E.R. Niblett, and R . A . Stacey, 1971. A review of Geophy-sical Studies in the Canadian Cordillera, Canadian Journal of Earth Sciences, Vol.8, p.788-801. Blakely, R.J . , 1996. Potential Theory in Gravity and Magnetic Applications, Cam-bridge University Press, Cambridge. Cecile, M.P. , D.W. Morrow, and G.K. Williams, 1997. Early Paleozoic (Cambrian to Early Devonian) tectonic framework, Canadian Cordillera, Bulletin of Ca-nadian Petroleum Geology, Vol.45, p.54-74. 137 References 138 Cerveny, V . , L A . Molotkov, and I. Psenvcik, 1977. Ray Method in Seismology, University of Karlova, Prague, Czechoslovakia. Chandra, N . N . and G.L. Cumming, 1972. Seismic Refraction Studies in Western Canada, Canadian Journal of Earth Sciences, Vol.9, p.1099-1109. Chapman, C .H . , 1978. A new method for computing synthetic seismograms, Geo-physical Journal of the Royal Astronomical Society, Vol.54, p.481-518. Clark, E . A . and F . A . Cook, 1992. Crustal-scale ramp in a Middle Proterozoic orogen, Northwest Territories, Canada, Canadian Journal of Earth Sciences, Vol.29, p.142-157. Clowes, R . M . (editor), 1997. L I T H O P R O B E Phase V Proposal - The Evolution of a Continent Revealed, L I T H O P R O B E Secretariat, The University of British Columbia, Vancouver, p.44-48. Clowes, R . M . , C A . Zelt, J.R. Amor, and R . M . Ellis, 1995. Lithospheric structure in the southern Canadian Cordillera from a network of seismic refraction lines, Canadian Journal of Earth Sciences, Vol.32, p.1485-1513. Coney, P.J. , D .L . Jones, and J .W.H. Monger, 1980. Cordillera suspect terranes, Nature, Vol.288, p.329-333. Cook, F . A . and A . J . van der Velden, 1993. Proterozoic crustal transition beneath the Western Canada sedimentary basin, Geology, Vol.21, p.785-788. Cook, F . A . , J .L. Varsek and E . A . Clark, 1991. Proterozoic craton to basin crustal References 139 transition in western Canada and its influence on the evolution of the Cordillera, Canadian Journal of Earth Sciences, Vol.28, p.1148-1158. Cook, F . A . , M . Dredge and E .A . Clark, 1992. The Proterozoic Fort Simpson struc-tural trend in northwestern Canada, Geological Society of America Bulletin, Vol.104, p.1121-1137. Cook, F .A . , J .L. Varsek and J .B. Thurston, 1995. Tectonic significance of gravity and magnetic variations along the L I T H O P R O B E southern Canadian Cor-dillera Transect, Canadian Journal of Earth Sciences, Vol.32, p.1584-1610. Cook, R .A . , A . J . van der Velden and K . W . Hall, 1999. Frozen subduction in Canada's Northwest Territories: L I T H O P R O B E deep lithospheric reflection profiling of the western Canadian Shield, Tectonics, Vol.18, p. 1-24. Creaser, B . and G. Spence, 1999. Crustal Seismic Velocity Structure of the North-ern Cordillera, Southern Yukon Territory - L I T H O P R O B E SNoRE Line 31, L I T H O P R O B E Report No.69, p.163-167. Fernandez Viejo, G. , R . M . Clowes, J.R., Amor, and the SNoRE Working Group, 1999a. Lithospheric Velocity Structure beneath the Slave Province and Wop-may Orogen from the SNorCLE Refraction Experiment, L I T H O P R O B E Report No.69, p.73-80. Fernandez Viejo, G., R . M . Clowes, and J.R. Amor, 1999b. Imaging the lithospheric mantle in northwestern Canada with seismic wide-angle reflections, Geophysical Research Letters, Vol.26, p.2809-2812. References 140 Firbas, P., 1981. Inversion of travel-time data for laterally heterogeneous velocity structure - linearization approach, Geophysical Journal of the Royal Astrono-mical Society, Vol.67, p. 189-198. Gabrielse, H . , 1967. Tectonic evolution of the northern Canadian Cordillera, Ca-nadian Journal of Earth Sciences, Vol.4, p.271-298. Gabrielse, H . , 1985. Major dextral transcurrent displacements along the North-ern Rocky Mountain Trench and related lineaments in north-central British Columbia, Geological Society of America Bulletin, Vol.96, p. 1-14. Gabrielse, H . , 1991. Structural Styles. Part F: Transcurrent Faults, Tintina and Northern Rocky Mountain Trench System, in Geology of the Cordilleran Oro-gen in Canada, edited by H . Gabrielse and C.J . Yorath, Geological Survey of Canada, Geology of Canada No.4, Vol.G-2, p.651. Gabrielse, H . and R .B . Campbell, 1991. Upper Proterozoic Assemblages, in Geology of the Cordilleran Orogen in Canada, edited by H . Gabrielse and C .J . Yorath, Geological Survey of Canada, Geology of Canada No.4, Vol.G-2, p. 127-150. Gabrielse, H . , J .W.H. Monger, J.O. Wheeler and C.J . Yorath, 1991. Tectonic frame-work. Part A : Morphogeological belts, tectonic assemblages, and terranes, in Geology of the Cordilleran Orogen in Canada, edited by H . Gabrielse and C.J . Yorath, Geological Survey of Canada, Geology of Canada No.4, Vol.G-2, p.16-28. Gorman, A .R . and T.J . Henstock, 1997. Southern Alberta Refraction Experiment References 141 ( S A R E X ) and Deep Probe Refraction Experiment 1995: Field Acquisition and Preliminary Data Processing Report, L I T H O P R O B E Report No.52, p.6-22. Gorman, A .R . and R . M . Clowes, 1999. Wave-field tau-p analysis for 2-d velocity models: Application to western North American lithosphere, Geophysical Re-search Letters, Vol.26, p.2323-2326. Gough, D . L , 1986. Mantle upflow tectonics in the Canadian Cordillera, Journal of Geophysical Research, Vol.91, p.1909-1919. Hammer, P.T.C. , R . M . , Clowes, and R . M . Ellis, 1999. Crustal structure of N W British Columbia and SE Alaska from seismic wide-angle studies: Coast Plu-tonic Complex to Stikinia, Journal of Geophysical Research, in press. Hoffman, P.F., 1989. Precambrian geology and tectonic history of North America, The Geology of North America - A n Overview, p.491-493. Huang, H . , C. Spencer, and A. Green, 1986. A method for the inversion of refraction and reflection travel times for laterally varying velocity structures, Bulletin of the Seismological Society of America, Vol.76, p.837-846. Hyndman, R.S. and T.J . Lewis, 1995. Review: The thermal regime along the south-ern Canadian Cordillera L I T H O P R O B E corridor, Canadian Journal of Earth Sciences, Vol.32, p.1611-1617. Hyndman, R.S. and T.J . Lewis, 1999. Geophysical consequences of the Cordillera-craton thermal transition in southwestern Canada, Tectonophysics, Vol.306, .References 142 p.397-422. Lewis, T .J . and R.S. Hyndman, 1998. Heat Flow Transitions along the SNorCLE transect: Asthenospheric Flow, L I T H O P R O B E Report No.64, p.92. Lewis, T . J . and R.S. Hyndman, 1999. High heat flow along the SNorCLE transect, L I T H O P R O B E Report No.69, p.83-84. L i , Y . and D.W. Oldenburg, 1998. 3D inversion of gravity data, Geophysics, Vol.63, p.109-119. Lowe, C , R . B . Horner, J .K . Mortensen, S.T. Johnson and C F . Roots, 1994. New geophysical data from the northern Cordillera: preliminary interpretations and implications for the tectonics and deep geology, Canadian Journal of Earth Sciences, Vol.31, p.891-904. Luetgert, J . , 1992. Macray-interactive two-dimensional seismic raytracing for the Macintosh, U.S. Geological Survey, Open File Report 92-356. Lutter, W . J . and R.L . Nowack, 1990. Inversion for crustal structure using reflections from the P A S S C A L Ouachita experiment, Journal of Geophysical Research, Vol.95, p.4633-4646. Lutter, W.J . , R .L . Nowack, and L .W. Braile, 1990. Seismic imaging of upper crustal structure using travel times from the PASSCAL Ouachita experiment, Journal of Geophysical Research, Vol.95, p.4621-4631. McMechan, G.A. and W.D. , Mooney, 1980. Asymptotic ray theory and synthetic References 143 seismograms for laterally varying structures: theory and application to the Imperial Valley, California, Bulletin of the Seismological Society of America, Vol.70, p.2021-2035. Mereu, R.F . , S.C. Majumdar, and R .E . White, 1977. The structure of the crust and upper mantle under the highest ranges of the Canadian Rockies from a seismic refraction survey, Canadian Journal of Earth Sciences, Vol.14, p. 196-208. Monger, J .W.H. , 1993. Canadian Cordilleran Tectonics: from geosynclines to crustal collage, Canadian Journal of Earth Sciences, Vol.30, p.209-230. Morozov, L B . , S.B. Smithson, and L. Hollister, 1999. Wide-angle seismic imaging across accreted terranes, southeastern Alaska and western British Columbia, Tectonophysics, Vol.299, p.281-296. Monger, J .W.H. and R . A . Price, 1979. Geodynamic evolution of the Canadian Cordillera - progress and problems, Canadian Journal of Earth Sciences, Vol.16, p.770-791. Mossop, G.D. and Shetsen, I. (compilers), 1994. Geological Atlas of the Western Canada Sedimentary Basin, Calgary, Canadian Society of Petroleum Geologists and Alberta Research Council, 510 p. Ross, G . M . , 1991. Precambrian basement in the Canadian Cordillera: an introduc-tion, Canadian journal of Earth Sciences, Vol.28, p.1133-1139. References 144 Ross, G . M . , R.R. Parrish, and D. Winston, 1992. Provenance and U-Pb geochrono-logy of the Mesoproterozoic Belt Supergroup (northwestern United States): im-plications for the age of deposition and pre-Panthalassa plate reconstructions, Earth and Planetary Science Letters, Vol.113, p.57-76. Sandiford, M . and M . Hand, 1998. Controls on the locus of intraplate deformation in central Australia, Earth and Planetary Science Letters, Vol.162, p.97-110. Sheriff, R .E . and L.P. Geldart, 1983. Exploration Seismology, Vol.2: Data Pro-cessing and Interpretation, Cambridge University Press, Cambridge, U K . Singh, S. and R .B . Herrmann, 1983. Regionalization of Crustal Coda Q in the Continental United States, Journal of Geophysical Research, Vol.88, p.527-538. Spence, G.D., K .P . Whittall, and R . M . Clowes, 1984. Practical synthetic seismo-grams for laterally varying media calculated by asymptotic ray theory, Bulletin of the Seismological Society of America, Vol.74, p.1209-1223. Spence, G.D., R . M . Clowes, and R . M . Ellis, 1985. Seismic structure across the active subduction zone of western Canada, Journal of Geophysical Research, Vol.90, p.6754-6772. Stacey, R .A . , 1973. Gravity anomalies, crustal structure, and plate tectonics in the Canadian Cordillera, Canadian Journal of Earth Sciences, Vol.10, p.615-628. Sweeney, J.F. , R . A . Stephenson, R .G . Currie, and J . M . DeLaurier, 1991. Tectonic References 145 Framework. Part C: Crustal Geophysics, in Geology of the Cordilleran Oro-gen in Canada, edited by H . Gabrielse and C.J . Yorath, Geological Survey of Canada, Geology of Canada No.4, Vol.G-2, p.39-58. Tatel, H.E. and M . A . Tuve, 1956. Seismic crustal measurements in Alaska, EOS, Transactions, Vol.37, p.360. Tempelman-Kluit, D.J . , 1977. Thrust relations in the northern Omineca Crystalline Belt and evidence for strike slip on the Tintina Fault (abstract), Geological Associations of Canada, Program with Abstracts, Vol.2, p.51. Thompson, R.I., E . Mercier and C F . Roots, 1987. Extension and its influence on Canadian Cordillera passive margin evolution, Continental Extensional Tec-tonics, Coward, Dewey and Hancock (editors), Geological Society of London Special Publication, Vol.28, p.409-417. Turcotte, D.L. and G. Schubert, 1982. Geodynamics: Applications of Continuum Physics to Geological Problems, John Wiley and Sons, Toronto, p.75. Villeneuve, M . E . , R . J . Theriault and G . M . Ross, 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, Vol.28, p.1003-1008. Webring, M . , 1985. SAKI: A Fortran program for generalized linear inversion of gravity and magnetic profiles, United States Department of the Interior, Geo-logical Survey, Open File Report 85-122. References 146 Wessel, P. and W.H.F . Smith, 1995. New version of the Generic Mapping Tools released, EOS Transactions, American Geophysical Union, Vol.76, p.329. White, D.J . , 1989. Two-dimensional seismic refraction tomography, Geophysical Journal International, Vol.97, p.223-245. White, D.J . and R . M . Clowes, 1990. Shallow crustal structure beneath the Juan de Fuca Ridge from 2-D seismic refraction tomography, Geophysical Journal International, Vol.100, p.349-367. Whittall, K .P . and R . M . Clowes, 1979. A simple, efficient method for the calculation of traveltimes and ray paths in laterally inhomogeneous media, Journal of the Canadian Society of Exploration Geophysicists, Vol.15, p.21-29. Winardhi, S. and R.F. Mereu, 1997. Crustal velocity structure of the Superior and Grenville provinces of the southeastern Canadian Shield, Canadian Journal of Earth Sciences, Vol.34, p.1167-1184. Woollard, G.P., 1969. Regional variations in gravity, American Geophysical Union Monogr. 13, p.320-341. Wu, X . , I. Ferguson and A . Jones, 1998. Electrical resistivity structure between the Nahanni terrane and Slave Province, L I T H O P R O B E Report No!64, p.93. Zelt, C .A . , 1999. Modeling strategies and model assessment for wide-angle seismic traveltime data, Geophysical Journal International, in press. Zelt, C A . and R . M . Ellis, 1988. Practical and efficient ray tracing in 2-D media References 147 for rapid traveltime and amplitude forward modeling, Canadian Journal of Exploration Geophysics, Vol.24, p. 16-31. Zelt, C A . and R . M . Ellis, 1989. Seismic structure of the crust and upper mantle in the Peace River Arch region, Canada, Journal of Geophysical Research, Vol.94, p.5729-5744. Zelt, C A . and R .B . Smith, 1992. Seismic traveltime inversion for 2-D crustal velo-city structure, Geophysical Journal International, Vol.108, p. 16-34. Zelt, C A . and D.A. Forsyth, 1994. Modeling wide-angle seismic data for crustal structure: Southeastern Grenville Province, Journal of Geophysical Research, Vol.99, p.11687-11704. Zelt, C A . and D.J . White, 1995. Crustal structure and tectonics of the southeast-ern Canadian Cordillera, Journal of Geophysical Research, Vol.100, p.24,255-24,273. Zelt, C A . and B .C . Zelt, 1998. Study of out-of-plane effects in the inversion of refraction/wide-angle reflection traveltimes, Tectonophysics, Vol.286, p.209-221. Appendix A Data and Related Results Appendix A displays all data and related results of the analysis for SNoRE Line 21. Traveltime-distance (T-X) plots have been bandpass filtered (0.5-10 Hz) and are shown with velocity reductions of 8.0 km/s [i.e., plotted as (T-X/8.0) versus X where X is shown as both offset distance and model distance (km)]. Synthetic seismogram displays use a wavelet derived from the waveform of a prominent well-defined Pg arrival from shot 2108. N W and SE on data plots signify northwest and southeast respectively. The appendix includes: • true relative amplitude record sections with a gain factor of X 1 8 5 [Figs. A.2, A.5, A.8,. . . , A.23]; • trace normalized record sections [Figs. A . l a , A.4a, A.7a,..., A.22a]; • comparison of traveltimes picked from the data and traveltimes calculated from the final velocity model [Figs. A . l b , A.4b, A.7b,..., A.22b]; • synthetic seismograms with a gain factor of X 1 ' 8 5 calculated for the final velocity model [Figs. A.3a, A.6a, A.9a,..., A.24a]; • amplitude versus offset plots for all phases for the observed data and synthetic seismograms [Figs. A.3b, A.6b, A.9b,..., A.24b]; 148 Appendix A. Data and Related Results 149 2 i Distance (km) 0 ~~50 100 150 200 250 300 350 400 450 Offset (km) -450 -400 -350 -300 -250 -200 -150 -100 -50 0 Figure A . l : a) Trace normalized plot for shot 2101; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data and Related Results 150 Figure A.2: True relative amplitude plot for shot 2101. Appendix A. Data and Related Results 151 Pg ( c a l c u l a t e d ) R l ( ca l cu la ted ) P n ( ca l cu l a ted ) distance (km) distance (km) distance (km) Figure A.3: a) Synthetic seismogram record section for shot 2101; b) Amplitude versus offset plot for observed and synthetic plots. Appendix A. Data and Related Results 152 NW SNoRE 2102 SE o J 1 1 1 1 1 Distance (km) rj 50 100 150 200 250 300 350 400 450 Offset (km) -350 -300 -250 -200 -150 -100 -50 0 50 Figure A.4: a) Trace normalized plot for shot 2102; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data and Related Results 153 Figure A.5: True relative amplitude plot for shot 2102. Appendix A. Data and Related Results 154 Figure A.6: a) Synthetic seismogram record section for shot 2102; b) Amplitude versus offset plot for observed and synthetic plots. Appendix A. Data and Related Results 155 NW SNoRE 2104 SE Distance (km) o 50 100 150 200 250 300 350 400 450 Offset (km) .300 -250 -200 -150 -100 -50 0 50 100 150 Figure A.7: a) Trace normalized plot for shot 2104; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data and Related Results Figure A.8: True relative amplitude plot for shot 2104. Appendix A. Data and Related Results 157 Figure A.9: a) Synthetic seismogram record section for shot 2104; b) Amplitude versus offset plot for observed and synthetic plots. Appendix A. Data and Related Results 158 SNoRE 2105 Distance (km) 0 50 100 150 200 250 300 350 400 450 Offset (km)-200 -150 -100 -50 0 50 100 150 200 250 Distance (km) 0 50 100 150 200 250 300 350 400 450 Offset (km)-200 -150 -100 -50 0 50 100 150 200 250 Figure A. 10: a) Trace normalized plot for shot 2105; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data, and Related Results 159 Figure A.11: True relative amplitude plot for shot 2105. Appendix A. Data and Related Results 160 Figure A. 12: a) Synthetic seismogram record section for shot 2105; b) Amplitude versus offset plot for observed and synthetic plots. 10. Figure A.13: a) Trace normalized plot for shot 2106; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data and Related Results Figure A.14: True relative amplitude plot for shot 2106. Appendix A. Data and Related Results 163 Figure A . 15: a) Synthetic seismogram record section for shot 2106; b) Amplitude versus offset plot for observed and synthetic plots. Appendix A. Data and Related Results 164 12 PmP PmP R ~ * l ! | W i * * * 1 ^ " " Pn JiJ l 50 100 150 Distance (km) 0 0 0 200 Offset (km)-100 -50 0 50 100 150 200 250 300 350 Figure A. 16: a) Trace normalized plot for shot 2107; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. 250 150 300 200 350 400 450 Appendix A. Data and Related Results 165 Figure A. 17: True relative amplitude plot for shot 2107. Appendix A. Data and Related Results 166 Figure A . 18: a) Synthetic seismogram record section for shot 2107; b) Amplitude versus offset plot for observed and synthetic plots. Appendix A. Data and Related Results 167 SNoRE 2108 Distance (km) 0 Offset (km) -50 Distance (km) 0 50 100 150 200 250 300 350 400 450 Offset (km) -50 0 50 100 150 200 250 300 350 400 Figure A.19: a) Trace normalized plot for shot 2108; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data and Related Results 168 Figure A.20: True relative amplitude plot for shot 2108. Appendix A. Data and Related Results 169 0 50 100 150 200 250 300 350 400 450 DISTANCE (km) b) Pg (observed) P m P (observed) Pn (observed) P g ( ca l cu l a ted ) P m P (ca l cu la ted ) P n ( c a l c u l a t e d ) distance (km) distance (km) distance (km) Figure A.21: a) Synthetic seismogram record section for shot 2108; b) Amplitude versus offset plot for observed and synthetic plots. Appendix A. Data and Related Results 170 12-10 co Q 2H 0 \ Distance (km) 0 50 100 150 200 250 300 350 400 450 Offset (km) 0 50 100 150 200 250 300 350 400 450 Figure A.22: a) Trace normalized plot for shot 2109; on the enlarged segment, picks are overlain as horizontal bars; b) Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes. Appendix A. Data, and Related Results Distance (km) 0 50 100 150 200 250 300 350 400 450 Offset (km) 0 50 100 150 200 250 300 350 400 450 Figure A.23: True relative amplitude plot for shot 2109. Appendix A. Data and Related Results 172 Figure A.24: a) Synthetic seismogram record section for shot 2109; b) Amplitude versus offset plot for observed and synthetic plots. Appendix B Inverse Problem Mathematics The traveltime along a ray path is denned as which is a nonlinear relationship given that the ray path is dependent on the velocity. This nonlinearity complicates the inverse problem and requires that the problem be linearized in order to be efficiently solved. To linearize the problem, the traveltime (t0) of a ray traveling through a continuous reference velocity field (v0) is considered such that t0(v0(x,z))= f j~~ rdl (B.2) JLo V0(X,Z) where Lo is the newly considered ray path. If we suppose that the velocity field in our problem, v(x, z), is simply a perturbation, 6v(x, z), of the reference field, v0(x, z), such that v(x, z) = v0(x, z) + 6v(x, z), the travel time t(v(x, z)) can be written as a Taylor series expansion about v0(x, z): t(v(x,z)) = t(v0(x,z)) + — 6v(x,z)+— — 6v(x,z) + 0(6v2(x,z)) (B.3) v v '' OV v-v0 Oil v=v0 OV v=vo Given that the travel time for the ray path is stationary as a first order approximation according to Fermat's Principle (i.e. If =0) and neglecting higher order terms, the ° a l j V=VQ perturbations in travel times are 173 Appendix B. Inverse Problem Mathematics 174 St = t(v(x,z))-t(v0(x,z)) = — I Sv(x,z). (BA) Qy \v=V(, It is important to note that the use of Fermat's principle can introduce second order errors in the traveltimes that can potentially alter the calculated velocity models significantly. Once equation B.4 is differentiated, the following linear relationship between 6t and 8v(x,z) is obtained « = -/ *?N*- (B-») JU v£(x, z) It should be noted that the problem is solved in R A Y I N V R with respect to velocity as opposed to slowness since it is simpler to parameterize the model using velocity gradient cells than slowness gradient cells (White and Clowes, 1990). Discretization of the Velocity Model The parameterization for R A Y I N V R ensures that the velocity field is linearly inter-polated for all areas of the model. The velocity field, v(x,z), can also be discretized in terms of basis vectors, 8j(x, z), that are related to the velocity points in the model such that N v(x>z) = 12Pi(x>z)vi' (B-6) where N is the number of velocity points specified in the velocity model and Vj are the ve-locities at those points. Given this reasoning, it can also be stated that the perturbations of the velocity field can be discretized as 6v(x,z) = f^3j(x,z)6vj (B.7) j=i which can be substituted into eq. B.5 to yield Appendix B. Inverse Problem Mathematics 175 st = -£ f (3j(x,z) L0 vl(x,z) dl SVJ. (B.8) This can also be written as ~N dt (B.9) which is equivalent to St = ASv (B.10) where A is the Jacobian matrix of sensitivities with elements dU/dvj which characterize traveltime and v3- is the j th model parameter such that A has the dimensions of (number of observations) x (number of model parameters). By inspection of eq. B.8, the elements of the matrix A are equivalent to Both St and A are computed analytically during the ray tracing step of the iteration. Solving the Inverse Problem The general inverse problem attempts to minimize a specific global objective functic 4>(m) subject to the constraint that St — A Sm. Since the model parameterization performed with respect to velocity, the term Sv and Sm are equivalent; thus, Sm will used for the remainder of this appendix to represent the perturbations in the model. The global objective function is the relationships between each data point and the model parameters. U is the i th observed ( B . l l ) 4>(m) =. <j)d(m) + f3(/>m(m) (B.12) Appendix B. Inverse Problem Mathematics 176 where <$>d and <f>m are the misfit and model objective functions respectively which are equal to • M™) = \\\W4tob'--F(m))\\2 and <pm(m) = ±||Wm(m - m r e / ) | | 2 (B.13) where Wd = diag(l/ai) and W m = diag(l/aj) are weighting matrices designed to make up for differences in standard deviations between the traveltimes and the model parameters. The global objective function for a perturbed model can be Taylor expanded such that (j>(m + 8m) = (j)(m) + (V<p(m))TSm + ^8mTWT<p(m)8m + R(m, 6m) (B.14) By neglecting the remainder term, R(m,8m), and setting 'Vsm</> —> 0, we are left with Newton's equation H(m)8m = -g{m) (B.15) where H(m) = WT<p(m) is the Hessian matrix and 5(m) = V<£(Tn) is the gradient matrix. By differentiating the components of equation B.13 with respect to m, the following relationships can be deduced g(m) = V<f>(m) = V<j>d{m) + 8V<f>m(m) = -ATWjWd8t + 3W^Wm(m - mref) (B.16) and using a Gauss-Newton step, H(m) = VV0(m) = V V ^ + BVV<t>m = ATWjWdA + 8WlWm (B.17) Appendix B. Inverse Problem Mathematics 177 which can be substituted into equation B.15 to yield ( A T W j W d A + 8W^Wm)Sm = ATWj\VdSt - 8W^Wm(m - mref). (B.18) Since WjWd = diag(lfai) and W^Wm = diag(l/<jj), the above equation can be rewrit-ten such that Sm = (ATCt1A + DC-1)-\ATCr1St-aC-1(m-mref) (B.19) where C t and Cm are the estimated data and model covariance matrices given by Ct = diag{of) and Cm = diag(a)). (B.20) D, an overall damping parameter which is equivalent to 8, is usually set to 1 based on previous trial and error experiments. This fixing of 8 is very significant as it forces an equal contribution to the global objective function from the misfit and the model objective functions. R A Y I N V R actually does not solve equation B.19. R A Y I N V R solves a simplified version of equation B.19 which is derived by calculating the least-squares solution of St = A Sm directly instead of using objective functions. This technique yields the solution Sm = (ATCr1A + DC-1)-1ATCr1St (B.21) where the last term of eq. B.19 is ignored. Equation B.21 can be determined using the procedure used to derive eq. B.19 by solving for a model objective function, <f>m equal to ' ^ = ^ ! ! ^ n ( M H 2 (B.22) Appendix B. Inverse Problem Mathematics 178 which will disappear during the differentiation. Equation B.21 only manages to minimize the perturbations in the model without constraining the complexity of the calculated model. This drawback can lead to crazy, unrealistic models after many iterations and does not provide as robust a solution as equation B.19. Despite these obvious sources of error, equation B.21 which has been the classic approach to solving seismic inversion for many years (Zelt and Smith, 1992; Lutter et al., 1990), continues to be the popular choice. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items