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

LITHOSPHERIC 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) McGill University  A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T THE  REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF  SCIENCE  in THE  FACULTY OF GRADUATE  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  THE  U N I V E R S I T Y O F BRITISH  COLUMBIA  December 1999 © Joanna K i m Welford, 1999  OF  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 V 6 T 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/wideangle 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 refraction data set shows good signal-to-noise ratios; primary phases (Pg, P m P 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 forward 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 modeling 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 1  INTRODUCTION  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  1.3  1.4 2  xv  12  Previous work in northeastern British Columbia  16  1.3.1  Geological Mapping  16  1.3.2  Geophysics in the Northern Canadian Cordillera  16  Experimental Objectives  21  DATA ACQUISITION A N D PROCESSING  22  2.1  The 1997 Slave-Northern Cordillera Refraction Experiment . .  22  2.1.1  23  2.2  Location  Data Acquisition  . . .  23 iv  2.3  2.4  3  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  Preparation of Data  31  2.3.1  Processing Routines  31  2.3.2  Gain  31  2.3.3  Updating of S E G - Y Headers  32  2.3.4  Merging Data and Creating Final Data Sets  32  Data Processing  33  2.4.1  True Relative Amplitude Scaling  33  2.4.2  Filtering . .  33  DATA DESCRIPTION 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 P m P Phase  38  3.1.4  Observations of P n 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 OF T H E VELOCITY MODEL  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  5  6  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  POTENTIAL FIELD DATA MODELING  84  5.1  Gravity Data and Modeling  84  5.1.1  3-d Inversion Modeling using G R A V 3 D  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  DISCUSSION  110  6.1  110  6.2  Comparison with other Geophysical Studies in the Region 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  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  vi  6.2.3  Interpretation of Variations along-strike of the Cordillera-Craton Transition Zone  6.2.4 6.3  129  Interpreted Lithospheric Structure for Northeastern British Columbia 132  Conclusions  134  References  137  Appendices  148  A Data and Related Results  148  B Inverse Problem Mathematics  173  vii  L i s 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  Modelfittingstatistics 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 P m P phase.  39  3.4  Trace normalized plot for shot 2102 showing an example of the P n 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  4.7  Ray coverage for P n phase  69  4.8  Comparison of observed and calculated traveltimes from shot 2105 . . . .  70  4.9  Trace normalized plot for shot 2101 showing P m P phase masked by earlier R l arrivals.  .  .  68  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 G R A V 3 D inversion  92  5.8  Depth Slices through the density anomaly model generated by G R A V 3 D  94  5.9  Slice through G R A V 3 D 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  6.4  118  Comparison between interpreted reflection results from Line 11 and velocity model for Line 21  6.5  123  Interpreted Foreland-craton transitions for the southern and northern Cordilleras  130  6.6  Interpreted evolution of Foreland-craton transition in the northern Cordilleral31  6.7  Interpreted cross-section for the Cordilleran-North American craton transition  133  A.1 Shot 2101. a) Trace normalized plot. Expected location of P m P 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 observed and synthetic plots A.4 Shot 2102.  151  a) Trace normalized plot; b) Comparison of observed and  calculated traveltimes . A.5 True relative amplitude plot for shot 2102  • .  152  .  153  A.6 Shot 2102. a) Synthetic seismogram record section; b) A V O plot for observed 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  -. . . .  A.9  Shot 2104. a) Synthetic seismogram record section; b) A V O plot for observed and synthetic plots  A. 10 Shot 2105.  156  157  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 observed and synthetic plots A.13 Shot 2106.  160  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) A V O plot for observed and synthetic plots A. 16 Shot 2107.  163  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) A V O plot for observed and synthetic plots A. 19 Shot 2108.  166  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) A V O plot for observed and synthetic plots  169  xii  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 observed 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  INTRODUCTION  The dramatic topographic transition from the low plains of the Western Canada Sedimentary 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 detailed 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 (SlaveNorthern 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 interpretation 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.  1.1  INTRODUCTION  2  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 . T h e numbers 11, 21, 22 a n d 31 correspond to the four refraction lines collected d u r i n g the S N o R E '97 E x p e r i m e n t . A l l 36 shots of the experiment are shown as stars. T h e dashed b o x outlines the p r o j e c t i o n used for a l l subsequent regional m a p s .  Inset m a p shows l o c a t i o n of the S N o r C L E Transect w i t h respect to N o r t h  America.  (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 ( 1 9 8 9 ) , 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 BeltPurcell 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 northeastern 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 reveal 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 ancestral 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 Cordilleran 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 mountainous 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 tectonism 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 separates rocks with  8 7  8 7  S r / S r line. This line 8 6  S r / S r ratios less than 0.706, which are associated with oceanic 8 6  Chapter 1.  crust, from rocks with greater crust.  8 7  7  INTRODUCTION  8 7  S r / S r ratios, which are associated with continental 8 6  S r / S r ratios act as indicators of basement composition since magma will take 8 6  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 affecting 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 extending 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 Mountain 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 framework 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), allows 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 ( W C S B ) , across the deformation front of the northern Cordillera into deformed North American rocks and across the Tintina Fault, the eastern boundary of accreted terranes. W i t h 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.  9  INTRODUCTION -130°  -125°  -120°  Insular Superterrane Alexander  E£H  Coast Belt Undivided metamorphics and intrusives Accreted Terranes N  Taku Stikinia  | H  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 *, ""•t.  Ancestral North America Morphogeological belts Faults Deformation Front  *fV_  SNoRE Line 21 Shotpoints and Seismographs (inline)  ""N^  Seismographs along Line 22 (broadside)  v  -130°  -125°  -120°  ^  Political Boundaries  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 characteristics. 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 MidProterozoic 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 Y u k o n / B C border 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 Columbia 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 extent 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 miogeocline 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 approximately 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 deformation 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.  13  INTRODUCTION  Geologic Assemblages  Plutonic and Ultramafic Rocks  Cenozoic sedimentary rocks Cenozoic volcanic rocks Mesozoic-Cenozoic sedimentary rocks Mesozoic sedimentary and volcanic rocks sedimentary rocks • Mesozoic Mesozoic volcanic rocks Paleozoic-Mesozoic sedimentary rocks Paleozoic-Mesozoic sedimentary and volcanic rocks Paleozoic-Mesozoic volcanic rocks Paleozoic sedimentary rocks I I Paleozoic metamorphic rocks : I 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  • • • •  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.  14  INTRODUCTION  -130°  -120°  -125°  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,  Cordilleran Deformation Front  1000  Level  -1000  -2000  -3000  -4000  -5000  .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 1.3.1  Previous work in northeastern British Columbia Geological Mapping  The Geological Survey of Canada has been collecting geological data from the Canadian 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, aeromagnetic and seismic surveys which have been concentrated primarily in the southern Cordillera with some studies near the Yukon/Alaska border and in the Northwest Territories (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 surrounding the Line 21 area. The Peace River Arch Seismic Experiment ( P R A S E ) 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 southeast 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 L I T H O 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 reflection 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 Y u k o n / B C 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 conducted 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 anomalies with shorter wavelengths which are believed to reflect shallower depths to magnetic  20  Chapter 1. INTRODUCTION  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 resistivity 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 approved 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 and heat flow averages 2  of 103.3 ± 2 1 . 0 m W / m for the northern Cordillera (north of 59° latitude). Given that 2  the heat flow measurements for Yellowknife, Ft. Providence and the Nahanni terrane are 54 m W / m , 100 m W / m and 118-136 m W / m 2  2  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.  1.4  INTRODUCTION  21  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 ProterozoicPaleozoic 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  DATA ACQUISITION A N D PROCESSING  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 designed 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 crosssection 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 transect. 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 Columbia, Victoria, Saskatchewan, Western Ontario, Manitoba, Alberta, Calgary, Stanford and Wisconsin (Eau Claire) took part in the experiment.  22  Chapter 2. DATA ACQUISITION AND PROCESSING  2.1.1  23  Location  The SNoRE '97 refraction/wide-angle reflection survey comprised four refraction/wideangle 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  th  and 15 , 1997, and recorded on August 1 5 , th  th  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 No. 2101 2102 2104 2105 2106 2107 2108 2109  Latitude  Longitude  58.2000 58.7917 58.8362 59.5761 59.7000 59.9750 60.0750 60.2000  122.7750 123.5333 125.0253 126.5400 127.2833 128.2000 128.9500 129.5167  Julian Day 228 228 228 229 229 228 228 228  Shot Time 6:10:00 9:00:00 6:00:00 6:20:00 6:30:00 6:30:00 9:20:00 6:40:00  Elevation (m) 563 423 670 548 670 623 630 788  Charge (kg) 3000 2400 1200 2000 2000 2000 2400 3000  No. of holes 3 2 1 2 2 2 2 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 G P S 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 G P S software can reduce these uncertainties to 5 m horizontally and 10 m vertically.  Chapter 2.  DATA ACQUISITION  AND  PROCESSING  26  A l l G P S 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 G P S 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 G P S . 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 G P S was performed after the experiment. Thirty-six sites for which precise locations were not obtainable due to insufficient satellite information available 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 G P S 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 deployers 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 P R S l s ; 20 Scintrex-EDA model PRS4s; 212 RefTeks and 186 SGRs. Both the P R S l s and the SGRs are one-component  Chapter 2. DATA ACQUISITION  AND  PROCESSING  seismographs while the RefTeks and PRS4s are three-component seismographs.  28  The  SGRs were predominantly used for broadside coverage throughout the experiment while the P R S 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 G S C 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 S E G - 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 instruments. 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  29  PROCESSING  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 PRS4 RefTek SGR  178 20 212 186  GSC GSC IRIS-PASSCAL USGS  0.1 to 100 Hz 0.1 to 100 Hz 1 to 125 Hz 2 to 200 Hz  8.3333 ms 8.3333 ms 8 ms 2 ms  Table 2.2: Seismic Instrument Information.  Geophone Type L-4 L-22 L-28  Made by Mark Products Mark Products Mark Products  Band Width 1 to 20 Hz 2 to 40 Hz 4.5 to 80 Hz  Sensitivity 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 evennumbered 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 P R S 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 P R S and the S G R 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 P R S and S G R 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 S E G - Y . To enable different data formats from different instrument types to be combined into one standardized S E G - 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  32  PROCESSING  examination of the amplitudes of traces from neighboring instrument types using the PRS Is as the standard. Instrument Type PRS1 PRS4 RefTek SGR SGR  Geophone Type L-4 L-4 L-28 L-4 L-22  Gain Applied 1 0.5 -3.289 x 10 0.584 x 10 -1.136 x 10 7  7  7  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 plotsec_wsegy routine.  format to S E G - Y format using the  Chapter 2. DATA ACQUISITION  2.4  AND  PROCESSING  33  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 spreading 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  DATA DESCRIPTION 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, P m P 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 P g 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, P g 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 P m P arrivals from shots 2105, 2106 and 2107 and P n 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  37  ANALYSIS  Muted Pg amplitudes?'  m\fx~m\t^ M l Distance (km) Offset (km)  D  fu nf sieSt ("K,mf )f  100 100  0 0  ° ™ -J50  5  0  ->nn -JOU  ,  0  0  150 150  , -dx s  n  1 5  200 200  °,nn -ill"  2  0  0  250 250  , ™ -1MJ  2 5 0  300 300  ,on -TO 3  0  0  350 350  4| ou  3  5  0  n  400 400  4  u  0  0  450 450  4  su  5  0  0  Distance (km) 0 offset (km)-200  Distance (km) 0 Offset (km)-150  50 -150  50 100  100 150 -100 - 50  100 -50  150 0  200 0  200 50  50  250  300 100  350 400 450 1 50 200 250  250 100  300 1 50  350 200  400 250  450 300  Figure 3.2: Example data plots where Pg amplitudes appear muted. Muted regions are indicated 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 P m P due to its high amplitudes,  Chapter 3.  DATA DESCRIPTION  AND  38  ANALYSIS  modeling results to be presented in the next chapter reveal that it is in fact a lower crustal reflector that masks the later P m P 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 P m P Phase  P m P 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. P m P reflections are evident for all shots along Line 21 except shots 2101 and 2102 where the P m P 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 P m P arrivals with strong amplitudes and distinct codas, shots 2104 and those to the southeast exhibit more complicated arrivals. Initially, the P m P 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 P m P 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 P n phase.  3.1.4  Observations o f 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 P n 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 P n 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  S u m m a r y of phase observations  The good quality data from Line 21 had a high enough S / N ratio such that the main phases (Pg, P m P 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 corresponding 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 P n 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 procedures 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 ( R l ) masked the P m P 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  Phase Pg  PmP  Pn  R  R l  AND  Shot Point 2101 2102 2104 2105 2106 2107 2108 2109 2104 2105 2106 2107 2108 2109 2101 2102 2105 2106 2107 2108 2109 2102 2104 2105 2106 2107 2108 2109 2101 2102 2104  44  ANALYSIS  NW 160 - 40 135 - 15 70 - 0 120 - 0 150 - 0 100-0 55 - 5 20 - 10 140 - 80 200 - 115 150 - 105 100 - 75  SE 5 - 60 0 - 130 0-100 0 - 135 0 - 130 0 - 145 5 - 170 85 85 70 95 80  -  155 165 170 175 170  470 - 160 390 - 155 195 190 150 190 195  -  255 300 360 405 435  150 - 75 50 - 30 75 - 15 75 - 30 45 - 130 35 - 55 60 - 140 225 - 80 175 - 90 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  Shot No. 2101 2102 2104 2105 2106 2107 2108 2109 Total  AND  91 (119) 150 (76) 127 (104) 195 (84) 263 (75) 210 (65) 182 (43) 142 (52)  R 0 (0) 97 (105) 31 (106) 47 (141) 51 (110) 91 (133) 16 (30) 28 (132)  1360 (81)  361 (103)  P  g  45  ANALYSIS  Rl 180 (116) 109 (120) 18 (151) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 307 (120)  PmP 0 (0) 0 (0) 67 (135) 151 (148) 123 (133) 128 (111) 95 (104) 96 (116) 660 (130)  Pn 250 (142) 204 (132) 0 (0) 18 (142) 68 (136) 156 (148) 166 (125) 161 (135) 1023 (137)  Total 521 (129) 560 (82) 243 (117) 411 (116) 505 (98) 585 (181) 459 (61) 427 (69) 3711 (112)  Table 3.2: Number of observations (from traveltime picks) for each shot point with corresponding 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.  3.3  DATA DESCRIPTION  AND  ANALYSIS  46  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 developed 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 modeling 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 original 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  49  ANALYSIS  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  t = JL  ANALYSIS  v(x, z)  50  (3.2)  where L is the ray path. This equation can be discretized for practical applications such that:  (3.3) 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  Chapter 3.  DATA DESCRIPTION  AND  51  ANALYSIS  Distance -V  V  Figure 3.7: Schematic diagram of RAYINVR parameterization. Solid dots represent boundary 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 automatically 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, v is the partial derivative of v in the x direction and v x  z  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 (rootmean-square) residual traveltimes and the % factor. If these have not reached the desired 2  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 R M S 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 observations, an acceptable R M S residual traveltime should be on the order of the uncertainty of the picks. Meanwhile, the statistical % approach evaluates the differences between 2  Chapter 3. DATA DESCRIPTION  AND  ANALYSIS  55  the observed and calculated traveltimes (the misfit) with respect to their uncertainties, according to  X = 2  E  [i  *'  }i  (3.5)  "i  i=l  where o~i is the standard deviation in the traveltimes. This x factor which is subsequently 2  normalized by the number of observations, should equal 1 for a model that reproduces the observed traveltimes within the uncertainty bounds of the picks, x values less than 2  1 indicate that the data have been over-fit and that the model contains structure that is not required by the data. Conversely, x values greater than 1 normally mean that 2  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 x values greater than 2  cases, a  x  3.3.3  Amplitude Modeling  2  factor as close to  1  1  (Zelt and Forsyth, 1 9 9 4 ) . In such  as possible is sought (Zelt,  1999).  When modeling refraction/wide-angle reflection data, it is important that the final velocity 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 experiments and the sensitivity of amplitudes to small-scale heterogeneity (Zelt and Forsyth, 1994),  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  56  Chapter 3. DATA DESCRIPTION AND ANALYSIS  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 parameterization 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  DEVELOPMENT OF T H E VELOCITY MODEL  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  X  Proceed to next layer  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 P m P and P n 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 intervals 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 A V O 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 allowable 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  70  J  1 0  ; ! —t 50  OF THE VELOCITY  : ! ! • ! ! ! ' ! • !  !  —r*  100  '—i—'  150  1—  200  !  !!  f  250  61  MODEL  Vertical Exaggeration 4.25:1 '[ ' 1 l I 1 i 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.  63  Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL  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.  0  50  100  150  i  200  1  r  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. A t 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 provided 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  66  MODEL  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.  70  J  1  1  0  50  — i 100  1  150  1  200  1  250  1  1  1  <—  300  350  400  450  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 C r u s t a l Layers and the M o h o  Layer 5 is approximately 15 km thick and extends to the base of the Moho from the location 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 constraints such as lower crustal reflections, the boundary between layers 4 and 5 was again modeled as a continuous velocity boundary and P m P 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, P m P arrivals (ray coverage shown in Fig. 4.5) were used to invert both velocity and  Chapter 4. DEVELOPMENT  60 4 70  67  OF THE VELOCITY MODEL  ® 0  50  100  150 200 250 Distance (km)  300 350  400  450  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 P m P 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 P n arrivals. While it is common practice to invert for Moho depths using both P m P and P n arrivals, only P m P 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 P n 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  70  J  1  1  0  50  - 1  100  68  OF THE VELOCITY MODEL  1  1  '  150 200 250 Distance (km)  1  1  300 350  1  400 450  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  1 0  1  50  OF THE  i  i  100  VELOCITY  150  200  i  250  69  MODEL  i  300  1  350  i  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 P n 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 easternmost portion of the line to attempt to obtain the P n 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 P n phase (Fig. 4.8) since they still had to travel through the higher upper mantle velocities 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,  70  Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL  14 12  10J  Pg  2H Distance (km) 0  ^  mmiugajia**"  Q  Offset (km).200  Pn  PmP  co8  50  -150  100  -100  150  -50  0  200  250  50  300  100  350  150  400  450  200 250  Figure 4.8: Comparison of observed (vertical bars with length equal to twice the pick uncertainty) and calculated (thin lines) traveltimes from shot 2105. All 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 ( R l ) that existed a few kilometers above the Moho and that masked the later P m P 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 P m P 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.  4.1.2  DEVELOPMENT  OF THE VELOCITY  MODEL  72  Spatial Resolution and Absolute Parameter Uncertainty  The velocity model that was developed in this thesis was arrived at using a minimumparameter/prior-structure model approach (Zelt, 1999) which signified that a priori information (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 important 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 parameterization 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 inversion 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 calculated. These calculated traveltimes are then interpolated to the receiver locations from  Chapter 4. DEVELOPMENT  OF THE VELOCITY  73  MODEL  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 traveltimes 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 x values and 2  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 complicated 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 P n arrivals indicated 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  V V V V V B B B V V  Surface Top of Upper Crust Within Upper Crust Top of Middle Crust Base of Lower Crust Moho west of 2104 R l Reflector Moho east of 2104 Top of Upper Mantle Base of model  1(T) 2(T) 3(T) 4(T) 5(B) 5(B) 5(B) 6(B) 7(T) 7(B)  Approximate lateral resolution (km) <15(40)-110 <15(40)-110 75 60 80 (120) 60 100 100 90 140  Approximate absolute uncertainty 0.1-1.0 (0.1-1.6) km/s 0.1-0.5 km/s 0.1-0.2 km/s 0.1-0.2 km/s 0.3 km/s 1.5-2.5 km 1.5-2.5 km 3.5 km 0.2-0.25 (0.25-0.4) km/s 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 P g ray coverage. This improvement in lateral resolution was accompanied by a decrease in absolute parameter 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 P m P traveltimes. Secondly, while the P m P arrivals provided a lateral resolution of 60 km on the boundary 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.52.5 km across the layer. East of shot location 2104, the lack of P m P arrivals required the depth of the Moho to be constrained solely by P n 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 interactive 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  77  MODEL  model (Fig. 4 . 1 0 ) bears little resemblance to the model parameterization (Fig. 4 . 2 ) , signifying 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 R M S residual traveltime of 120 ms, its normalized x value of 2 . 9 9 1 and its ability 2  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 % value of 2 . 9 9 1 was the lowest achieved for Line 2 1 2  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  (1998)  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 120 ms R M S residual traveltime  was so close to the overall pick uncertainty of 112 ms that the introduction of more complexity would have incorporated unnecessary structure into the model. Pick Uncertainties (ms)  R M S Travel Time Residual (ms)  Normalized  4.801  Phase  Picked Traveltimes  . #of Calculated Traveltimes  Pg PmP Pn R Rl Total  1360  1336  81  106  660  658  130  105  1.217  #of  x Misfit 2  1023  1018  137  152  2.494  361  207  103  49  1.075  307  302  120  119  2.253  3711  3521  112  120  2.991  Table 4 . 2 : Model fitting statistics for each phase.  Chapter 4. DEVELOPMENT  50  0  OF THE VELOCITY  150  100  200  250  78  MODEL  Distance (km)  400  300 350  450  b) ESE  WNW  Foreland Belt  Omineca Belt 9  8 |  7  6  5  Ancestral N.A. 4 I  2  0  E  -20  g-40 -60 H 150  200 250 300 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 4.10: Final velocity structural model for Line 21, a) plotted with no vertical exaggeration; 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.  79  Chapter 4. DEVELOPMENT OF THE VELOCITY MODEL  Shot No. 2101 2102 2104 2105 2106 2107 2108 2109  #of Picked Traveltimes 521 560 243 411 505 585 459 427  #of Calculated Traveltimes 511 482 228 402 484 543 449 422  Pick Uncertainties (ms) 129 82 117 116 98 181 61 69  R M S Travel Time Residual (ms) 112 165 104 137 108 98 114 103  Normalized x Misfit 1.354 4.505 2.271 3.546 3.439 2.794 3.709 2.130 2  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  80  OF THE VELOCITY MODEL  ESE  WNW  Foreland Belt  Omineca Belt 9  Ancestral N.A.  6  8  5  0  7  -20  J--40  -60 A 0  50  100  lntensity=function(ray density) Full Colour=64 rays White=0 rays  150  200  250  300  350  400  450  Distance (km)  to Q64  £  (km/s) l!9 3!o 4!o 5!8 6.0 6l2 6l4 8.8 7!0 7l5 8.\ 8A 8.9  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. :  P m P 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 P m P 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  83  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 FIELD DATA 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 improve 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/cm  3  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 gridding 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  85  relevance to this study are: 1) the Tintina Fault which is defined by a gravity expression 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 approach 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 interpretation 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 G R A V 3 D 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 distancesquared 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  88  MODELING  500 400 ^  300  200 100 100  200  300  500  400  Eastings (km)  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 refraction 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, G R A V 3 D 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  89  MODELING  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),  200  300  400  I -2200  0  400  700  |  500  |  600  q (meters)  900 1100 1300 1500 1800 2600  0  100  200  300  500  600  I (blocks removed)  ITT" 0  400  *  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 G R A V 3 D 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  90  MODELING  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/cm  z  with respect to a reducing density of 2.67 g/cm  3  (i.e. the upper mantle  was assigned a density of 3.20 g/cm ). Within the crust, the assigned density anomalies 3  increased linearly with depth from -0.37 g/cm  3  at the surface up to 0.33 g/cm  base of the crust corresponding to a density range of 2.30-3.00 g/cm  3  3  at the  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  50  FIELD  100  150  DATA  200  250  91  MODELING  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  0.2  0.3  0.4  0.5  0.6  Distance (km) (g/cm**3) -0.7  -0.5  -0.4  -0.3  -0.2  -0.1  0.0  0.1  Figure 5.6: Reference density anomaly model used in the inversion.  Chapter 5. POTENTIAL  FIELD DATA  92  MODELING  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 x factor of the misfit was 2  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 % factor was readjusted, the inversion was 2  rerun. With the increased % factor, the algorithm converged after 32 iterations. The pre2  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 constraints 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 refraction line (Fig. 5.9) was taken to examine whether the resolved velocity blocks corresponded with distinct regions of anomalous densities. Comparison of the final velocity structural model (Fig. 4.10) with the G R A V 3 D 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 G R A V 3 D 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  94  MODELING  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  0  100  DATA  200  95  MODELING  300  400  500  600  Eastings (km)  P P W H P P P - T - ' J ^ -160 -90 -40 -30 -20 -10 0  10 20 30 40  (mgals) 120  WNW  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.  5.1.2  POTENTIAL  FIELD DATA  MODELING  96  2.5-d Forward Modeling using S A K I  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 structural constraints derived from the seismic refraction model and the 3-d gravity inversion, 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 (Barton, 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 P m P wide-angle reflections were included in the model as constraints. The density model parameterization used in S A K I involves a polygonal parameterization 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 velocity. 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 experimentally 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  -160.0 ' -40.0  FIELD DATA MODELING  1  60.0  1  ' • ' 160.0 260.0 Model Distance (km) 1  98  ! 360.0  '  '— 460.0 1  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  * Observed - Calculated  ii ii n ii H II in  -80 -90 |-100 E  -110 -120 100  Model Distance (km) 200 300  400  ^ 20 E o. 40 5  •  60  B l 2.73  275 2.80 2.82 2.85  j 2.90 2.95  3.15 3.00 3.10  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 G R A V 3 D inversion, successfully reproduces the increase in Bouguer gravity values at the eastern end of the line. A t 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  101  MODELING  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  cr = pgz  (5.1)  zz  where o~ is the normal stress acting at the base of the column due to the overlying rock, zz  p is the column's density, g is the gravitational acceleration (9.81 m/s ) and z is the 2  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/cm . 3  Next the depths of  the five main boundaries, z\ to z are determined; see, Fig. 5.12 for definitions. Finally, 5  the equation  {zi + z )pb + (z4 - z )p •+ (z - z )p 2  2  c  5  4  m  = k  (5.2)  Chapter 5. POTENTIAL  FIELD DATA  102  MODELING  is solved, where k is the mass per unit area in kilograms per square meter at the base of each column. -1  Pc Z| 4  •m  Figure 5.12: Model for isostatic balance calculations. p , p and p represent the main blocks making up each column. z is the height of topography, z is the depth below which all densities are greater than or equal to the Bouguer reducing density of 2.67 g/cm , z is the depth to the base of the crust if z = z , z is the depth to the base of the crust including the crustal root, and z is the maximum depth for lateral changes in density. b  x  c  m  2  3  z  x  0  4  5  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 compensation. 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. A t 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 equilibrium; 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  104  MODELING  w 2.06x10'  400  200  cd  S  Distance along the line (km)  3  T  TfT  3  2.9  b)  S  2  8  ft F-  ' i . ' '  11 l|  ii  CP  1 L r-'  «  1 _ l  I | r  I - ,  I  f  ,1  •  'Ll I._ J'  2.7  -  ,  II r 1  If  j  K  2.6  3  400  200  Distance along the line (km) -1  o QfJ  )1  c  3.4 3.3  CO  C  CD  Q jU C  (0  S  3.2  "  1  .  ~  •I  .'. ^  f 1  .i  !!  3.1 - >  1  ;  J \  |i  "  »  ' i n  •  -  1 1  •  J  1  r' 1  -1 J  r  1  ~  1  Jl  r - - I  1  1  1  i  200  i  i  1  'I'I  i  ' r.; F  h  11 II  1  1 —  r  II  : -  i  1  1  i  r ll  i i  i  i  " i  i  i  -  '  r  -  Z  1' 1 i |  1 ^ 1  1  r  400  Distance along the line (km)  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  5.3  FIELD DATA MODELING  105  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 magnetic 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 magnetization 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  -760 -200 -100 -50  0  60  MODELING  106  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 (corresponding 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 declination, 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  Original Data  108  MODELING  Reduced to pole  Figure 5.15: Comparison of aeromagnetic data before and after reduction to pole. The refraction line has been overlain on both plots for reference.  This positive anomaly may be the result of realignment of magnetic minerals from sediments filling in the fault scar due to metamorphic remelting. At the Yukon/British Columbia border where the fault has no surface expression, an enlargement of the reducedto-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  -760 -200 -100 -50  0  MODELING  60  109  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 geophysical 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 i n the College F j o r d of Alaska and another 13 shots detonated near Skagway, Alaska. T h e 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 B C / s o u t h e r n Y u k o n area. Their results revealed that the crust i n the region is relatively uniform w i t h an average upper crustal P-wave velocity of 6.1 k m / 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. W h i l e M o h o reflections were detected from their survey, no depth estimates were suggested and they did not observe a significant increase i n velocity w i t h depth i n the lower crust.  The P R A S E 85 Experiment Southeast of Line 21 i n northeastern B r i t i s h C o l u m b i a and northwestern A l b e r t a between 54°N and 59°N and 122°W and 115°W, the Peace River A r c h Seismic E x p e r i m e n t ( 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 R i v e r A r c h , a regional geological structure with anomalous vertical displacements compared to the rest of the W C S B (Zelt and Ellis, 1989). W h i l e the P R A S E experiment's proximity to Line 21 provides some hope that direct 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 basement 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 S N o R E '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 P m P 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  B r i t i s h Columbia. T h e projected model ( F i g . 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 k m . 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 ( F i g . 1.1) which lies to the northwest of Line 21 and extends 300 k m N E to S W i n the central Y u k o n , crosses, from east to west, the deformed miogeoclinal rocks of the O m i n e c a Belt, the T i n t i n a Fault, the displaced continental margin and the accreted terranes of the Intermontane Belt. Despite the T i n t i n a 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 T i n t i n a Fault along b o t h lines. Preliminary results from Creaser and Spence (1999) indicate a general increase i n upper and middle crustal velocities westward across the T i n t i n a Fault and a corresponding increase i n M o h o depths from 32 k m to 37 k m . T h e low upper crustal velocities seen near the surface of the Omineca Belt along Line 31 agree w i t h the overall low velocities at the surface at the western edge of Line 21. However, the modeled increase i n M o h o depths across the T i n t i n a 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. D u r i n g 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 i n t i n a Fault. T h e analysis and modeling of these broadside data will begin once the inline profile analyses for lines 21 and 22 are complete.  aptei6.  DISCUSSION  (iu>l) qjdsa  Chapter 6.  118  DISCUSSION -130  -  -120"  125'  a)  -130  b)  OMINECA  INTERMONTANE I  ~  Stikinia  100  125"  ANCESTRAL N.A.  FORELAND  |CC CA Cassia SM  200  300  j  FS  500  400  Distance (km) 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, C A , 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  6.1.2  119  Comparison with Refraction Studies in the Southern Cordillera  Several extensive seismic refraction surveys have been conducted in the southern Cordillera 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 summers of 1962-1969 from Greenbush Lake, British Columbia, to Swift Current, Saskatchewan. 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 achieved 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 i n velocities eastward. A deepening of the Moho to ~ 45 km beneath the Foreland 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 velocities 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. P n 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 S C o R E '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 Refraction 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.  6.1.3  DISCUSSION  122  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.  123  Chapter 6. DISCUSSION  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, R l 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.  124  DISCUSSION  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 while in the Cordillera they range between 2  80-100 m W / m (Hyndman and Lewis, 1999). 2  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 ; west of Yellowknife, N W T , situated on the Slave Craton, the heat flow is 53 2  mW/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.  125  DISCUSSION  have recorded heat flow values of 89 ± 17 m W / m  2  across the terranes and up to 118-  136 m W / m in the Nahanni terrane itself (Lewis and Hyndman, 1999). These values 2  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 (Lewis and Hyndman, 1999), an extremely high 2  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 . This value suggests that 45 3  m W / m of the heat flow in the region originates within the upper 10 km of the crust. 2  These results offer some explanation as to how the upper mantle beneath the Foreland Belt could have remained cold (inferred from the high P n 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 upper 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 6.2.1  Final Velocity Model within the Context of the Regional Tectonics 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 B C / Y u k o n 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 P n 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 anomalously 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 material, 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 sediments 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, relatively 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 Transition 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.  130  DISCUSSION  crust and upper mantle are inferred to be colder and thus, stronger. Since the westfacing 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  / V  js  \  WCSB  North American Craton  Blfes*-*-  Northern Cordillera  -—  J  WCSB^^-^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.  131  Chapter 6. DISCUSSION Proterozoic  Low Velocity Upper Mantle High Velocity Upper Mantle I  Middle Jurassic  Passive— b - - i - * . Low V e l o ^ i i y - ^ - , ~-. Upper Mantle  Cratonic North America  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 transition suggest a hotter and weaker lower crust and upper mantle. As the miogeoclinal sediments 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 mantle 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  Ancestral North America  Foreland Belt  Fort Nahanni : Terrane i  Simpson Terrane  Proterozoic Mantle  North American Lower Lithosphere  F i g u r e 6.7: Interpreted cross-section for the C o r d i l l e r a n - N o r t h A m e r i c a n craton t r a n s i t i o n . L o c a t i o n s o f S N o R E L i n e 21 a n d S N o r C L E reflection line 11 are shown across the top o f the d i a g r a m . T h e l a t e r a l extents of cratonic N o r t h A m e r i c a , the F o r e l a n d B e l t a n d the O m i n e c a B e l t ( O . B . ) are shown across the top of the d i a g r a m . T h e inferred dip of the T i n t i n a F a u l t is i n d i c a t e d w i t h a t h i c k dashed line at the western end of L i n e 21. D u r i n g t h e c o m p r e s s i o n a l events t h a t r e s u l t e d f r o m t h e a c c r e t i o n of t h e w e s t e r n C o r d i l l e r a n b e l t s d u r i n g t h e M i d d l e J u r a s s i c ( c a . 170 M a ) , these passive m a r g i n m e t a s e d i 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 u p t h r u s t onto the c r a t o n i c m a r g i n to form the Foreland B e l t . However, i n contrast to the southern C o r d i l l e r a where t h e r a m p e d edge o f c r a t o n i c N o r t h A m e r i c a p l a y e d a d i r e c t r o l e i n t h e u p t h r u s t o f p a s s i v e m a r g i n s e d i m e n t s , t h e r a m p e d edge o f c r a t o n i c N o r t h A m e r i c a l a y east o f t h e n o r t h e r n R o c k y M o u n t a i n s o f t h e F o r e l a n d B e l t a n d as s u c h , s e v e r a l h u n d r e d s o f k i l o m e t e r s o f t h e F o r t S i m p s o n B a s i n w e r e left u n d e f o r m e d . T h e u p t h r u s t i n g of passive m a r g i n s e d i m e n t s t o f o r m t h e t h i c k e r c r u s t o f t h e F o r e l a n d B e l t resulted i n a depression i n the M o h o due to a weaker crust and upper m a n 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 v e l o c i t y u p p e r m a n t l e m a t e r i a l t o e i t h e r side o f t h e F o r e l a n d B e l t a n d t h e u p w e l l i n g o f  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 amplitude modeling, a velocity structural model for Line 21 of the SNoRE '97 refraction/wideangle reflection experiment was constructed 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. 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 interpretation 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 directly 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 response 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 remained 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 margin 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 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.99124. 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 Astronomical Society, Vol.87, p. 195-208. Berry, M . J . , W.R. Jacoby, E.R. Niblett, and R . A . Stacey, 1971. A review of Geophysical 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, Cambridge 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 Canadian 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, Geophysical 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 structural 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 Cordillera 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 Northern 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 Wopmay 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.  140  References  Firbas, P., 1981. Inversion of travel-time data for laterally heterogeneous velocity structure - linearization approach, Geophysical Journal of the Royal Astronomical Society, Vol.67, p. 189-198. Gabrielse, H . , 1967. Tectonic evolution of the northern Canadian Cordillera, Canadian Journal of Earth Sciences, Vol.4, p.271-298. Gabrielse, H . , 1985. Major dextral transcurrent displacements along the Northern 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 Orogen 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 framework. 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.1628. 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 Research 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 Plutonic 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 southern 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 Cordilleracraton 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 P A S S C A L 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 introduction, 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 geochronology of the Mesoproterozoic Belt Supergroup (northwestern United States): implications 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 Processing 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 seismograms 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  145  References  Framework. Part C: Crustal Geophysics, 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.39-58. Tatel, H . E . and M . A . Tuve, 1956. Seismic crustal measurements in Alaska, E O S , 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 Tectonics, 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. S A K I : A Fortran program for generalized linear inversion of gravity and magnetic profiles, United States Department of the Interior, Geological 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 velocity 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 southeastern Canadian Cordillera, Journal of Geophysical Research, Vol.100, p.24,25524,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.209221.  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 S E 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  2  149  i  Distance (km) 0 ~~50 Offset (km) -450 -400  100 -350  150 -300  200 -250  250 -200  300 -150  350 -100  400 -50  450 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  Figure A.2: True relative amplitude plot for shot 2101.  150  151  Appendix A. Data and Related Results  Pg ( c a l c u l a t e d )  R l (calculated)  P n (calculated)  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  NW  o  J  Distance (km) Offset (km)  152  SNoRE 2102  1  rj  1  50 -350  1  100 -300  -250  150 -200  SE  1  200 -150  250 -100  1  300  350 -50  400 0  450 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  Figure A.5: True relative amplitude plot for shot 2102.  153  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  NW  Distance (km) o Offset (km) .300  155  SNoRE 2104  50 -250  100 -200  150 -150  200 -100  250 -50  SE  300 0  350 50  400 100  450 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.  158  Appendix A. Data and Related Results  SNoRE 2105  Distance (km) 0 Offset (km)-200  Distance (km) 0 Offset (km)-200  50 100 150 -150 -100 -50  50 100 150 -150 -100 -50  200 0  250 50  200 0  250 50  300 100  300 100  350 400 450 150 200 250  350 400 450 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  Figure A.11: True relative amplitude plot for shot 2105.  159  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.  164  Appendix A. Data and Related Results  12  PmP  PmP  R  ~*  l ! |  Pn  JiJl  Wi*** ^"" 1  400 450 350 300 250 150 Distance (km) 0 50 100 100 150 200 200 150 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.  Appendix A. Data and Related Results  Figure A. 17: True relative amplitude plot for shot 2107.  165  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 Offset (km) -50  50 0  100 50  150 100  200 150  250 200  300 250  350 300  400 350  450 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  Figure A.20: True relative amplitude plot for shot 2108.  168  169  Appendix A. Data and Related Results  0  50  100  150  200 250 DISTANCE ( k m )  300  350  400  450  b) Pg (observed)  P m P (observed)  P n (observed)  Pg (calculated)  P m P (calculated)  P n (calculated)  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.  170  Appendix A. Data and Related Results  12-  10  co  Q  2H \  0 Distance (km) Offset (km)  0 0  50 50  100 100  150 150  200 200  250 250  300 300  350 350  400 400  450 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) Offset (km)  0 0  50 50  100 100  150 150  200 200  250 250  300 300  350 350  400 400  Figure A.23: True relative amplitude plot for shot 2109.  450 450  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 (t ) of a ray traveling through a continuous 0  reference velocity field (v ) is considered such that 0  f  t (v (x,z))= 0  0  JLo  j~~ rdl  (B.2)  V (X,Z) 0  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, v (x, z), such that v(x, z) = v (x, z) + 6v(x, z), the travel 0  0  time t(v(x, z)) can be written as a Taylor series expansion about v (x, z): 0  t(v(x,z)) = t(v (x,z)) + —  6v(x,z)+—  0  v  v  ''  OV v-v  Oil  0  —  6v(x,z) + 0(6v (x,z)) 2  (B.3)  v=v OV v=vo 0  Given that the travel time for the ray path is stationary as a first order approximation according to Fermat's Principle (i.e. If °  a  l  j  =0) and neglecting higher order terms, the V=VQ  perturbations in travel times are  173  174  Appendix B. Inverse Problem Mathematics  = — I  St = t(v(x,z))-t(v (x,z)) 0  Qy  Sv(x,z).  (BA)  \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 * JU  (B-»)  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 interpolated 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( > ) = 12Pi( > ) i' x  z  x  z v  (-) B 6  where N is the number of velocity points specified in the velocity model and Vj are the velocities 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^3 (x,z)6v j  j=i  which can be substituted into eq. B.5 to yield  j  (B.7)  175  Appendix B. Inverse Problem Mathematics  st = -£  f (3j(x,z) L  vl(x,z)  ~  dt  0  dl  SVJ.  (B.8)  This can also be written as N  (B.9)  which is equivalent to  (B.10)  St = ASv  where A is the Jacobian matrix of sensitivities with elements dU/dvj which characterize the relationships between each data point and the model parameters. U is the i t h observed traveltime and v - is the j t h model parameter such that A has the dimensions of (number 3  of observations) x (number of model parameters). B y inspection of eq. B.8, the elements of the matrix A are equivalent to  (B.ll) 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  4>(m) =. <j)d(m) + f3(/> (m) m  (B.12)  176  Appendix B. Inverse Problem Mathematics  where <$> and <f> are the misfit and model objective functions respectively which are d  m  equal to  • M™) = \\\W4t '--F(m))\\ ob  and  2  where W = diag(l/ai) and W d  m  <p (m) = ±||W (m m  m  - m  r e /  )||  2  (B.13)  = 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)) Sm + ^8m W <p(m)8m + R(m, 6m) T  T  (B.14)  T  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) = W <p(m) is the Hessian matrix and 5(m) = V<£(Tn) is the gradient T  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> {m) + 8V<f> (m) = -A WjW 8t T  d  m  d  + 3W^W (m - m ) m  ref  and using a Gauss-Newton step,  H(m) = V V 0 ( m ) = V V ^ + BVV<t> = A WjW A T  m  d  + 8WlW  m  (B.17)  (B.16)  177  Appendix B. Inverse Problem Mathematics  which can be substituted into equation B.15 to yield  ( A W j W A + 8W^W )Sm = A Wj\V St T  - 8W^W (m - m ).  T  d  m  Since WjWd = diag(lfai)  d  and W^W  m  (B.18)  ref  = diag(l/<jj), the above equation can be rewrit-  m  ten such that  Sm = (A Ct A T  where C and C t  m  + DC- )-\A Cr St-aC- (m-m )  1  1  T  1  (B.19)  1  ref  are the estimated data and model covariance matrices given by  C = diag{of) t  C  and  = diag(a)).  m  (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 version of equation  B.19. R A Y I N V R solves a simplified  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 = (A Cr A T  1  + DC- )- A Cr St 1  1  T  (B.21)  1  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> equal to m  '  ^  = ^!!^n(MH  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