Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Crustal velocity structure in the southern coast belt, British Columbia, from a seismic refraction survey O'Leary, Deirdre 1992

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

Item Metadata


831-ubc_1993_spring_oleary_deirdre.pdf [ 6.2MB ]
JSON: 831-1.0053004.json
JSON-LD: 831-1.0053004-ld.json
RDF/XML (Pretty): 831-1.0053004-rdf.xml
RDF/JSON: 831-1.0053004-rdf.json
Turtle: 831-1.0053004-turtle.txt
N-Triples: 831-1.0053004-rdf-ntriples.txt
Original Record: 831-1.0053004-source.json
Full Text

Full Text

CRUSTAL VELOCITY STRUCTURE IN THE SOUTHERN COAST BELT, BRITISH COLUMBIA, FROM A SEISMIC REFRACTION SURVEY  by Deirdre O'Leary B.Sc. (Geophysics), McGill University, 1990  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOPHYSICS AND ASTRONOMY  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA December 1992 ©Deirdre O'Leary, 1992  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.  (Signatur  Department of  G el:v(1y g 1%cs an t A c÷ron 0 r-t  The University of British Columbia Vancouver, Canada  Date  DE-6 (2/88)  Dec. 24 ) Iqq2  e  ..  ABSTRACT The Coast Belt is one of five morphogeological belts of the Canadian Cordillera. It was created through the complex tectonic processes that accommodated the Mesozoic accretion of the allochthonous Insular superterrane to the Intermontane superterrane, then the leading edge of western North America, and the subsequent overprinting of the suture zone by granitic intrusions associated with east-dipping subduction. As part of the LITHOPROBE Southern Cordillera transect, seismic refraction data were recorded along a 350 km long strike profile in the Coast Belt. An iterative combination of two-dimensional travel time inversion and amplitude forward modelling was used to interpret crust and upper mantle P-wave velocity structure. The structural model features a thin (0.5 to 3.0 kin) near-surface layer with an average velocity of 4.40 km/s, a significant vertical velocity gradient and lateral variations. This uppermost stratum overlies a crustal velocity structure with three layers of approximately 10.0 km each. The upper and middle crust have average velocities of 6.20 and 6.30 km/s, respectively. The lower crust consists of an upper and lower unit with average velocities of 6.50 and 6.65 km/s, respectively. Beneath the lower crust lies Moho, the crustmantle boundary, which is modelled as an approximately 2 km thick transitional layer with an average depth to its upper boundary of 34.5 km and a maximum depth of 35.5 km in the southeast. The transitional layer exhibits lateral velocity variations (7.40 to 7.80 km/s) and overlies an upper mantle which is poorly defined in terms of velocity, although the data indicate relatively high values (> 8.05 km/s). Interpretations of LITHOPROBE reflection data and other refraction lines which cross, or nearly cross, the Coast Belt profile suggest that this profile is within the region of the collision zone between the Insular and Intermontane superterranes. These interpretations indicate two models for the collision zone: (i) a crustal delamination model in which the Insular superterrane was displaced along east-vergent faults over the terranes below and (ii) a crustal imbrication model in which imbrication of Insular and Intermontane rocks occurs ii  throughout the crust. The latter model involves the presence of thick layers of Insular material beneath the Coast Belt refraction profile. However, the velocity model indicates predominantly Intermontane material, thereby favoring the crustal delamination model and the eastward extent of the Insular Belt being west of the refraction profile. From comparisons of the refraction velocity model with the reflection data, the top of the Moho transition zone corresponds to the top of a prominent band of reflectivity which constitutes the reflection Moho. When combined with the seismic reflection as well as other geophysical studies, three likely sources for wide-angle reflections from the upper, middle, and lower crust, as observed in the refraction data, are structural features (e.g. fault zones), Ethological contrasts (e.g. batholithic granites to accreted volcanics) and the transition zones at the top and bottom of a region of layered porosity in the crust.  Table of Contents ABSTRACT ^  ii  List of Tables ^  vii  List of Figures ^ Acknowledgements ^ Chapter I^INTRODUCTION ^ 1.1 Tectonic History — Terrane Accretion ^  viii ix 1 1  1.2 Geological and Geophysical Background ^  6  1.2.1 Geology of Southern Coast Belt ^  6  1.2.2 Previous Geophysical Studies ^  8  1.3 Experimental Objectives ^  10  Chapter II^DATA ACQUISITION AND PROCESSING ^ 11 2.1 The Refraction Experiment ^  11  2.1.1 LITHOPROBE Southern Cordillera Transect ^  11  2.1.2 SCoRE '89 & SCoRE '90 ^  11  2.1.3 Line 10 — Acquisition Geometry of the Field Data ^ 12 2.1.4 Shot and Receiver Site Positioning ^  13  2.2 Instrumentation ^  13  2.3 Processing ^  14  2.3.1 Initial Data Reduction ^  14  2.3.2 Filtering ^  14  2.3.3 Amplitude Corrections ^  15  iv  Chapter III^DATA ANALYSIS AND INTERPRETATION PROCEDURES . . 17 3.1 Initial Observations and Interpretations ^  17  3.1.1 Data Characteristics and Quality ^  17  3.1.2 Travel Time Picking ^  36  3.2 Interpretation Methods ^  37  3.2.1 The Modelling Algorithms ^  37  3.2.2 Inversion mathematics ^  40  3.2.3 The Modelling Techniques ^  44  3.2.4 Modelling the Refraction Data ^  45  3.2.4a The Upper Crust ^  45  3.2.4b The Middle Crust ^  53  3.2.4c The Lower Crust and Upper Mantle ^  55  Chapter IV^INTERPRETATION ^ 4.1 The Final Model ^  61 61  4.1.1 Primary Features of the Velocity Model ^  61  4.1.2 Non-uniqueness of the Model ^  63  4.2 Spatial Resolution and Absolute Parameter Uncertainty ^ 63 Chapter V^DISCUSSION AND CONCLUSIONS ^ 66 5.1 Comparison with LITHOPROBE Refraction Interpretations ^ 66 5.2 Comparison with Other Geophysical Studies — Reflection, Heat Flow, Gravity, and Geomagnetic ^  71  5.2.1 Comparison with LITHOPROBE Reflection Data and Interpretations . . . ^ 71 5.2.2 Heat Flow Data and Interpretations ^  83  5.2.3 Gravity Data and Interpretations ^  84  5.2.4 Geomagnetic Data and Interpretations ^ 5.3 Conclusions ^ Appendix A^Terrane Descriptions ^  vi  87 88 98  ^  List of Tables 1. Noisy traces deleted from true amplitude data sections ^  16  2. Number of observations (with corresponding uncertainties) for each shotpoint . . ^ 31 3. Source-receiver offsets of identified phases ^  34  4. Source-receiver offsets and reduced arrival times of critical points of PmP ^ 36 5. Q structure used in amplitude modelling ^  39  6. Inversion results for the final velocity model ^  49  7. Estimated lateral resolution and absolute parameter uncertainty ^ 65 8. Comparison of interpreted velocity structure for Line 10 and Line 2 of SCoRE'89 ^ 67 9. Comparison of interpreted velocity structure for Line 10 and Line 3 of SCoRE'89 ^ 68 10. Comparison of interpreted velocity structure for Line 10, Line 1 of ScoRE'89 and Line I of VISP '80  vii  70  List of Figures 1. Map of the Canadian Cordillera showing accreted terranes ^ 2 2. Map of study area showing Line 10 and other regional seismic profiles ^ 4 3. Collisional model of the creation of the Coast Belt ^  5  4. Location map of Line 10 showing the geological features of the surrounding area ^ 7 5. Comparisons of observed and calculated data for SP-46 ^ 18 6. Comparisons of observed and calculated data for SP-47 ^ 21 7. Comparisons of observed and calculated data for SP-48 ^ 23 8. Comparisons of observed and calculated data for SP-49 ^ 25 9. Comparisons of observed and calculated data for SP-50 ^ 27 10. Comparisons of observed and calculated data for SP-51 ^ 29 11. Final velocity model and schematic interpretation ^  46  12. Location of nodes used in inversion process ^  47  13. Ray tracing diagrams ^  48  14. Refraction comparisons — 1–D velocity versus depth profiles ^ 69 15. Reflection data from LITHOPROBE profiles 88-13,88-14, and 88-17 ^ 73 16. Comparison of refraction model with results of 88-13 interpretation ^ 76 17. Comparison of refraction model with results of 88-14/88-17 interpretation ^ 77 18. Reflection interpretations ^  81  19. Bouguer gravity anomaly map ^  85  vu  '  Acknowledgements I wish to thank my supervisor Dr. Ron Clowes for his deft "coaching" abilities. I would also like to thank my co-supervisor Dr. Robert M. Ellis and Dr. Garry Clarke for critically reviewing this thesis. To my fellow linemates Barry Clint Zelt (#3), Mike Perz (#18), and Brad Isbell (#X): you made the GPRF room a place of final solutions, insane humour, and impeccable safety. My thanks also go to manager John Amor, for managing both the system and my nerves. The other departmental team members have made this a great place to study — I am grateful to all of them for their friendship. And, mesdames et messieurs, for their constant encouragement, the MOLSON/GRPF three star selection: 1. Les premiere etoiles, the first stars, my MOM and DAD! 2. La deuxieme etoile, the second star, DAVE HADLEY! 3. La troisieme etoile, the third star, CLAIRE DAT! SCoRE '90 was carried out by about 40 participants under the co-leadership of Dr. E. R. Kanasewich (U. of Alberta) and Dr. R. M. Ellis (U.B.C.). I thank all of them for their efforts in the field. Dr. Colin Zelt provided the ray-tracing inversion code used in this study. SCoRE '90 data were acquired primarily with funds from the NSERC Collaborative Special Project and Program grant in support of LITHOPROBE, with additional funds from the Geological Survey of Canada and NSERC Research grants of E. R. Kanasewich, R. M. Ellis and R. M. Clowes. Financial support for my student stipend and analysis costs derived from NSERC Research and Infrastructure grants to R. M. Clowes. Additional financial support was provided by Petro-Canada Inc.  ix  1 I. INTRODUCTION 1.1 Tectonic History — Terrane Accretion The Cordillera of North America is one of Earth's great orogenic belts. Until the 1960's, it was believed that Cordilleran orogenesis was the result of geosynclinal deposition (Monger 1992). The advent of the plate tectonic hypothesis coupled with the completion of regional stratigraphic analysis of supracrustal rocks led to a new theory of Cordilleran evolution based upon terrane accretion. This theory asserts that the Cordillera has grown 500 km westward from the ancestral North American margin, which was formed by extension, rifting, and sea floor spreading in Early to Late Proterozoic time (1600-570 Ma). The scope of the processes involved in this tectonic growth include (i) the deposition of substantial thicknesses (>10 km) of sedimentary and volcanic rocks in a passive margin setting, beginning in the Middle Proterozoic, (ii) the amalgamation and accretion of volcanic island arc and oceanic terranes, and associated capacious magmatism through Mesozoic and Cenozoic time, and (iii) very large transcurrent displacements during Cretaceous and Cenozoic time (Gabrielse and Yorath 1989). The term "terrane" refers to a block of the Earth's crust which holds a geological record distinct from those of surrounding terranes; the term has no significance pertaining to the origin of a crustal block, or to its past or present position. By definition, the boundaries between adjacent terranes are faults. The description of terranes as "accreted" defines them as those which became a part of ancestral North America at an advanced point of their tectonic development. "Amalgamated" terranes, or "superterranes", consist of at least two terranes which were joined prior to their accretion to North America (Gabrielse and Yorath 1989). The map shown in Figure 1 delineates the terranes which make up the Canadian portion of the Cordillera. This present structure is considered to be the result of the accretion of two allochthonous superterranes — the Insular superterrane and the Intermontane superterrane —  FIGURE 1: Map of the Canadian Cordillera showing the accreted terranes in different patterns and labelled with circled letters (from Clowes 1989). The list of terranes is shown at the bottom left, with stars indicating terranes that were formed near the west coast of the ancient craton; their relationships to North America are not certain. The other terranes are believed to have been formed on oceanic crust and then accreted to the ancient margin. The larger rectangle (solid line) shows the study area of the LITHOPROBE Southern Cordillera Transect; the inner rectangle (dotted line) shows the map area of Figure 2. Descriptions of the terranes discussed in the text are provided in Appendix A.  3  to the leading edge of the westward-travelling North American Plate prior to mid-Cretaceous time (e.g. Monger and Price 1979, Coney et al. 1980). The Insular superterrane, which was positioned adjacent to ancestral North America, comprises the Alexander and Wrangellia terranes, while the Intermontane superterrane to the east consists of the Stikine, Bridge River, Cache Creek, and Quesnel terranes. These terranes are briefly described in Appendix A. The Coast Belt is one of the five morphogeological belts of the Canadian Cordillera (inset, Figure 2); their boundaries are largely coincident with terrane boundaries (Gabrielse and Yorath, 1989). Numerous detailed explanations of the tectonic processes involved in creating the Coast Belt through the accretion of geologic terranes have been provided (e.g. Davis et al. 1978, Monger et al. 1982, Gehrels and Saleeby 1985, and van der Heyden 1992). In a currently widely-accepted model, Monger et al. (1982) described the metamorphic and structural character of the Coast Belt as the result of the mid-Cretaceous collision between the Insular superterrane and the Intermontane superterrane, which was by then attached to North America. According to this hypothesis, the Insular-Intermontane collision caused the closure of an intervening ocean basin, and the development of a collisional suture in its place. This suture was in turn overprinted by the granitoid intrusions of the Coast Belt, a magmatic arc linked to the east-dipping subduction beneath the new west coast of the North American continent. Most recently, van der Heyden (1992) has proposed a variation of this collisional model, whereby the Coast Belt was formed as an Andean-type magmatic arc superimposed upon a Cordilleran structure resulting from the mid-Jurassic accretion of a superterrane composed of Wrangellia, Alexander and Stikinia to ancestral North America. Although the collision model of the Coast Belt orogeny (which is illustrated in Figure 3) is generally accepted, the nature of the superterrane collision zone has been obscured by the magmatic intrusions of the Coast Belt, as well as the other regional geological complexities mentioned below.  4  FIGURE 2: Location map showing shotpoints (triangles) and receiver locations (small squares) for Line 10 of SCoRE '90 and Lines 1, 2, and 3 of SCoRE '89. Other seismic surveys in the region are shown by clashed lines; these include the refraction Line I of VISP '80, and the LITHOPROBE reflection profiles 88-13, 88-14, 88-17 and 88-18. The map area is shown by the rectangle outlined with dotted lines in Figure 1. The inset shows the morphogeological belts of the Canadian Cordillera.  5  FIGURE 3: Schematic illustration of the collision model for the evolution of the Coast Belt: Cretaceous magmatic arc, with associated subduction west of the Insular superterrane (from van der Heyden 1992).  6  1.2 Geological and Geophysical Background The Coast Belt runs for 1700 km between latitudes 49° N and 62° N, is 100-200 km wide, and reaches elevations above 4 km (although most peaks are in the range of 2-3 km) (Monger and Journeay 1992). It consists mostly of Late Jurassic to Early Tertiary granites (80-85% by area according to Roddick 1983), and its dominant structural features are Mid-Cretaceous to Early Tertiary in age. Due to the vast amount of plutonism, along with the high-grade metamorphism and complex deformation that characterize the Coast Belt (CB), the orogenesis and early tectonic past of this morphogeologic belt have proven difficult to discern.  1.2.1 Geology of Southern Coast Belt The area of interest for this study (see Line 10, Figures 2 and 4) lies in the southern CB between latitudes 49° N and 52° N, more than 200 km east of the Cascadia subduction zone which marks the east-clipping subduction currently taking place beneath the western continental margin. Based upon terrane affiliation and the distribution of dated plutonic rocks (Friedman and Armstrong 1992), the southern CB can be divided into eastern and western parts, both of which were traversed by Line 10 (Figure 4). Along the eastern margin of the Insular superterrane, there lies a Late-Cretaceous, westvergent contractional fault belt which exists in both the northern and southern CB (Monger and Journeay 1992). The Coast Belt Thrust System (CBTS) of southwestern British Columbia has been identified as a part of this contractional system (Journeay and Friedman, in press); Line 10 passed through the CBTS. On Figure 4 four major faults of the CBTS in the region of the present study have been identified. The southwestern CB (WCB, Figure 4) is characterized by low grade metamorphism and a mainly plutonic nature, and therefore probably acted as a rigid crustal block during Late Cretaceous west-directed thrust faulting (Journeay and Friedman, in press). To the east, the  7  -125'  -124'  ASCoRE '90 shot point  -123'  -122'  ;:Z5 Cadwallader T: rrane fig! Stikine Terrane Bridge River Terr ne Shuksan Terrane arrison Terrane ^Plu o• Rocks ost-Terrane Accretion Intrusions =Undivided Meta orphic Assemblages Gambier Overla Assemblage  011  52'  -121'  - 52'  Other Post-Terrane Accretion Overlap Assembla es  NTER ONTANE B T - 51'  5V  50'  - 50'  \100 ilometers 49'  49' -125'  -124'  - 123'  -122'  -121'  FIGURE 4: Location map showing shotpoints (triangles) and receiver locations (small squares) for Line 10. Heavy lines indicate the divisions between the Insular, Coast, and Intermontane belts, and lines show the positions of the major faults of the CBTS (Coast Belt Thrust System) in the region (CCBD denotes the Central Coast Belt Detachment). The Harrison Fault divides the Eastern Coast Belt (ECB) from the Western Coast Belt (WCB) (R. Friedman, pers. comm.). The regional geology is shown in the area surrounded by the dashed line (Wheeler et al. 1991).  8  WCB is thought to have been underplated beneath the southeastern CB (ECB) along westvergent thrust faults of the CBTS (Monger and Journeay 1992) during tectonic adjustments linked with the impingement of the Insular superterrane in Late Cretaceous time (Journeay 1990). The ECB (Figure 4) is characterized by the presence of oceanic Bridge River terrane rocks. The rocks range widely in age from at least 160 to 350 million years, and have therefore been identified as a former large ocean basin (Monger and Journeay 1992). The Line 10 profile traversed mainly granitic plutons created during post-terrane accretion intrusion, as well as slivers of terrane materials and some rocks excluded from terrane classification, as shown in the region in Figure 4 enclosed with a dashed line. In the northwest, the Stikine and Cadwallader terranes were crossed, as were the Gambier overlap assemblage (comprising arc volcanics and local rift volcanics) and another undifferentiated overlap assemblage (characterized by foredeep and southwesterly-derived clastic wedge material, foredeep marine shales, and arc volcanics). In the center and southeast, where it runs along the Harrison fault, Line 10 passed over pods of the Cadwallader terrane, the Gambier assemblage, and undivided metamorphic assemblages (Figure 4). The terranes of the CB are described briefly in Appendix A.  1.2.2 Previous Geophysical Studies In comparison with the other belts of the Canadian Cordillera, the Coast Belt (CB) has until recently been a frontier area in terms of geophysical studies. The pioneer seismic research effort in the CB was the interpretation of Berry and Forsyth (1975) of refraction data recorded prior to 1971. However, as the data were acquired with large shot and receiver spacings, the resolution of both lateral velocity variations and changes in velocity with depth were poor. The velocity model in the vicinity of Line 10 that resulted from this study has a crustal velocity of 6.0 km/s to a depth of 20 km, a lower crust with velocities of 6.1 km/s at the top and 6.9 km/s at the bottom, an upper mantle velocity of 7.8 km/s, and a crustmantle boundary at 33 km depth. A subsequent surface-wave study by Wickens (1977) in  9  the southern Canadian Cordillera that was generally consistent with that of Berry and Forsyth (1975) lacked data in the center of the southern CB. Although the Berry and Forsyth (1975) experiment was the only seismic survey coincident with the present area, previous refraction measurements for the north (Johnson et al. 1972) as well as the east (Cumming et al. 1979) and west (McMechan and Spence 1983, Spence et al. 1985, and Drew and Clowes 1990) have also been interpreted. The two former studies were not important to the current one as the Johnson et al. (1972) effort was restricted by the limitations described for the Berry and Forsyth (1975) case, and profiles recorded similarly to Line 10 have yielded updated results for the region to the east (Zelt et al. 1992a, 1992b). In the context of linking known features of the subsurface region to the west with the structure determined from Line 10, it is important to consider the Spence et al. (1985) interpretation of the refraction data from a refraction line (VISP '80 Line I, see Figure 2) which ran from Vancouver Island to the western CB. Principal features of the model of Spence et al. (1985) include upper crustal velocities of 6.4 km/s to 6.7 km/s (to a depth of 18 km), and a Moho depth of 38 km; however, the model is considered well-constrained only to 18 km depth. Other types of geophysical studies performed in the region of interest include geomagnetic (Caner 1970), magnetic (Coles and Currie 1977), and gravity (Riddihough 1979) surveys. The resolution of the structural information provided by these potential field data is limited in comparison with the detail which is of interest for the present study. The SCoRE '89 and '90 (Southern Cordillera Refraction Experiment) projects which encompassed Line 10 were run as part of the LITHOPROBE Southern Cordillera Transect (Clowes 1989, Clowes et al. 1992). LITHOPROBE is a multidisciplinary program designed to supplement the existing geological and geophysical data, in order to elucidate the processes involved in the westward growth of the North American continent. The profiles which  io  coincide with the present study are reflection profiles 88-13, 88-14, and 88-17, and refraction Lines 1, 2, and 3 of SCoRE '89 (see Figure 2). Heat flow (Lewis et al. 1988, 1992) and magnetotelluric (Jones et al. 1992) studies which formed part of the LITHOPROBE Southern Cordillera Transect activities were also undertaken. The elements of these and the aforementioned surveys which are of particular relevance will be discussed within the framework of the results determined in this thesis. 1.3 Experimental Objectives SCoRE '89 and '90 were designed to complement LITHOPROBE seismic reflection data acquired in 1985 and 1988 (Cook et al. 1988, 1991, 1992; Varsek et al. 1992) as they offered a means of extending the reflection results beyond 2—D interpretations. This would fulfill the general objective of providing structural models for the southern Canadian Cordillera to tie the models which exist for the convergent margin (e.g. McMechan and Spence 1983, Spence et al. 1985, Drew and Clowes 1990) with the interior of the continent. More specifically, the ultimate objectives put forth for the SCoRE '89 and '90 surveys were to study the velocity field, crust and upper mantle structure, spatial range of terranes, fault geometry and westernmost extent of cratonic basement in the southern Canadian Cordillera. This was to be achieved by the integrated interpretation with the seismic reflection data and other geophysical and geological studies. The specific objective for the present study is to analyse the refraction data of the SCoRE '90 Line 10 profile in order to map crust-mantle transition and intracrustal boundaries, and to integrate this interpretation with coincident seismic reflection and refraction data, gravity and heat flow profiles, and magnetotelluric information in order to unravel some of the mysteries surrounding the formation of the Coast Belt.  11 II. DATA ACQUISITION AND PROCESSING 2.1 The Refraction Experiment 2.1.1 LITHOPROBE Southern Cordillera Transect As part of LITHOPROBE's multidisciplinary scientific studies in the Southern Cordillera Transect, the extensive seismic refraction program referred to as SCoRE was can -ied out during the summers of 1989 and 1990. Line 10, the focus of this study, was recorded as one part of SCoRE '90.  2.1.2 SCoRE '89 & SCoRE '90 The first phase of the refraction experiment, SCoRE '89 (Figure 2), focussed on the Intermontane Belt and the Coast Belt, with an incursion into the eastern Insular Belt to provide a closer tie with the VISP '80 experiment (Zelt et al. 1992a). Two shots also were detonated west of Vancouver Island to extend the VISP refraction Line I across the CB. SCoRE '89 was an international collaboration with participants from Canadian universities, the Geological Survey of Canada, the United States Geological Survey, and Japan. It represented the first thorough application of the spatial seismic refraction recording (S 2 R2 ) method (Kanasewich and Chiu 1985). By recording shots from each apex of the triangular array, the center of each arm, and the center of the array, a set of inline and broadside data was obtained. Inline profiles are interpreted as 2—D crustal sections, whereas the suite of broadside data provide the basis for a tomographic study of the interior of the triangle, in this case focussed on the Fraser River fault system, the boundary zone between the Coast Belt and the Intermontane Belt. The second phase of the experiment, SCoRE '90, completed studies in the CB including Line 10, and extended coverage eastward. This also was a cooperative, international effort, involving the University of Alberta, the University of British Columbia, the Geological Survey of Canada, the U.S. Geological Survey , the University of Victoria, and the University  12  of Saskatchewan. It spanned a major portion of the Canadian Cordillera south of 52° N latitude. Line 10, which was approximately 350 km long and ran entirely within the CB, was the westernmost profile. The remainder of the project included east-west Lines 7 (400 km long) and 9 (360 km long) which spanned the Intermontane, Omineca Crystalline, and Foreland Belts of the Cordillera; north-south Line 8 (340 km long) which formed an alongstrike profile in the Omineca Belt; and Line 6 (230 km long) located just north of the US border and deployed to obtain broadside recording of shots along Line 8 (Burianyk et al. 1992).  2.1.3 Line 10 — Acquisition Geometry of the Field Data Line 10 extended northwest from Harrison Lake in the western CB to north of Chilco Lake in the eastern CB (Figure 2). Based on the experimental plans of shot points at - 50 ,  km intervals and the limitations of access imposed by the available road system, Line 10 included 6 shot points (SP) (Figures 2 and 4). Large shots were detonated at the endpoints: 2500 kg of explosives were used at SP-46, and 1800 kg at SP-51. An 800 kg shot was fired at SP-49 near the center of the line, while smaller sources were used at the intervening shotpoint locations (400 kg at SP-50, 200 kg at SP-48 and SP-47) in order to provide detailed information on shallow crustal structure. Three hundred portable seismographs were placed at intervals varying from 1 to 3 km along the profile. The receiver locations were determined by marking off the intervals on a straight line drawn between the shot points on a map, and then projecting these sites onto the available road system for ready access (Burianyk et al. 1992). Line 10 represented the only possible strike line within the CB along which such a survey could be carried out at reasonable cost. Nevertheless, many stations northwest of SP-49 and SP-50 had to be deployed by helicopter due to a lack of roads, and those along Chilco Lake were deployed by boat.  13  2.1.4 Shot and Receiver Site Positioning All of the shot points and 94% of the receiver locations were located using two types of GPS (Global Positioning System) navigational instruments: the Trimble Advanced Navigation Sensor GPS, and the Magellan GPS NAV 100 PRO. Two-dimensional fixes (i.e. 3 satellites used in the positioning solution) were used to position most coordinates, as this is the most accurate approach where altitude is known. Altitudes were read from topographic maps with a scale of 1:50 000, and the resulting measurement errors were approximately 100 m for latitude and longitude, and 20 m for altitude. When poor satellite visibility made GPS positioning difficult, the alternative was to read the locations from the 1:50 000 topographic maps, which resulted in errors of 100-150 for latitude and longitude (B. Isbell, personal communication). Due to a limited supply of instruments, Loran C positioning was used for a number of helicopter deployments at the northeastern portion of the line; these were then corrected using redundant GPS measurements so that accuracy levels were comparable to that of GPS (Burianyk et al. 1992).  2.2 Instrumentation The recording instruments consisted of 180 GSC/LITHOPROBE EDA model PRS1 digital refraction seismographs and 120 USGS analog cassette recorders, which were connected to Mark Products L-4A 2 Hz vertical seismometers. The sampling rate for the PRS-1 s was 120 samples/s, while the analog USGS data were digitized at 200 samples/s (Burianyk et al. 1992). For convenience of discussion, the data sets recorded on the PRS-ls and the analog systems will be referred to as the GSC and USGS data, respectively. The relative velocity responses of these instruments have been given by Zelt et al. (1990). Each shot point had one to three drill holes holding charges of Nitropel®, a nitroglycerinebased pellet explosive, at depths of 30-42 m. The detonation time of the charges was controlled by a clock device which uses the same satellite-based timing system (GOES) as the  14  recording instruments. The precision involved is such that instrument timing errors are considered negligible (I. Asudeh, personal communication); in other words, the beginning of each recorded trace is assumed to be the exact moment of shot detonation.  2.3 Processing 2.3.1 Initial Data Reduction Digital field data were initially available on 9—track tape in SEGY-LDS format (Spencer et al. 1989), with the GSC and the USGS data forming separate data sets. The first step in processing was conversion of the data recorded by the USGS instruments to the standardized format of the GSC data. The data were redigitized to 120 Hz from 200 Hz, gains were corrected to units of nanometers/s, and reformatting from SEGY format 3 (i.e. 2 byte fixed point) to SEGY format 1 (4 byte floating point) was carried out. Finally, 20 seconds of data (sample values of zero) were added to make the USGS traces equal in length to the GSC traces. At this stage, the headers of both the GSC and USGS traces were edited to update shot and receiver site locations. Corrections were therefore necessary to dependent header values such as source-receiver offset. This initial data reduction was done by the University of Alberta Seismology Laboratory (Burianyk et al. 1992). Additional processing was designed to extract the most information possible from the data.  2.3.2 Filtering Spectral analyses were performed on representative signal and noise recordings from the refraction data. This testing revealed that the power spectrum of the desired signal was within the 2-12 Hz range, whereas that of the noise was mainly outside of this range. An eightpole Butterworth bandpass filter (Kanasewich 1981) was applied to the data in an attempt  15  to improve the signal-to-noise ratio. These filtered sections were used concurrently with the respective unfiltered sections to make travel time picks.  2.3.3 Amplitude Corrections Amplitude corrections were applied to the data to account for spherical spreading Each trace was scaled by a factor proportional to d 15 , where d is the shot-receiver distance. This process enhanced the amplitude of the far-offset traces with respect to those near the source. There was one exception to this procedure: since the amount of energy released at SP-46 was so large, it was necessary to scale the corresponding traces proportionally to d 1 •. Otherwise, the amplitudes of the near-offset traces would be so large that separate phases could not be distinguished. There was no amplitude variation observed between traces recorded by the two types of seismographs, and since the velocity responses of each were nominally the same, no further amplitude corrections were necessary. The final step in preparing the data for interpretation was singling out traces which were so noisy that they obscured the signal visible on adjacent traces when the data were displayed with the amplitude corrections described above. For the purpose of these displays (Figures 5a-10a), these traces were deleted (see Table 1); this was not necessary when the data were plotted in common maximum amplitude format (Figures 5d-10d).  16  Shot Point SP-46 SP-47 SP-48 SP-49 SP-50 SP-51  Trace Numbers of Deleted Traces 147-151,161 55,58,69,75,78,97,147-151,161,193 55,58,69,91,97,147-151,161 198,247 147-151,161,178 135,147-151,161,178,189  Table 1. The trace numbers of those traces which had to be deleted (due to their excessive noisiness) in Figures 5a — 10a are specified for each shot point. Note that the traces are numbered 1-290 from left to right for each shot point in the data sections shown in the figures.  17  III. DATA ANALYSIS AND INTERPRETATION PROCEDURES 3.1 Initial Observations and Interpretations 3.1.1 Data Characteristics and Quality Acquisition of the field data resulted in a data set of 1734 traces. The six record sections are displayed in Figures 5a-10a. Of these 1734 traces, approximately 1100 contain primary and secondary arrivals which could be identified and timed. Distances in Figures 5-10 refer to model distances, with zero distance coinciding with SP-51. References to shot points accompanied by a "NW" or "SE" refer to data recorded to the northwest or southeast of the shot point. In the text, offset distances will be mentioned; these refer to distances from the shots to the receivers. The branches of the phases that were identified and used in the modelling procedure are labelled in Figures 5a-10a; in addition, Figure 5e shows the two phases which are not observable on a relative amplitude section unless the traces of SP-46 are scaled proportionately to d 1.5 . In total, eight phases were identified: Ps, Pg, and Pn, which represent energy which refracted back to the surface through the near-surface, upper crust, and upper mantle, respectively; and R1, R2, R3, PmP, and R, which denote energy which reflected back to the surface. The R1 reflector is at the base of the upper crust, R2 is at the base of the mid-crust, R3 is within the lower crust, and PmP corresponds to energy reflected from the crust-mantle boundary (see Figures 5a-10a). R is a sub-Moho reflection phase visible only on the far-offset traces of the SP-46 section when the appropriate amplitude scaling is applied (see Figures 5d and 5e). Due to low signal-to-noise ratio, inadequate shot-receiver offset, or nonexistence of a particular phase, not all phases are observed on all record sections. The quality of the data is usually good, with the following exceptions: SP-47 NW (50-110 km offset, Figure 6a), and SP-50 SE (120-180 km and 215-255 km offset, Figure 9a). Table 2 gives a synopsis of  18  SE  DISTANCE (k g/1111■1/1111M, AI^:11.111,  0  50  100  150^200  250  300  350  it.;,!. ar.-:.. L„  2.: I 1 .^g■^I , PrilLIIIIIMI 11011131111 ilia..mairimmatimmolun min  I  1E'1 pirgiliiti,  1•011till 1.  1116111411111114101111111M11111113IGNIGANNAIMIRIMIC15:111111MINIffeffloimiliminummtorast  12  10  (f) 0 0  50  100  150  200  250  300  350  21  FIGURE 6:(a) Observed record section for SP-47 (see Figure 5 caption for plotting parameters and general description). (b) Synthetic section for SP-47 (see Figure 5 caption for scaling and wavelet).  0^50^100^150^200^250^300^350 DISTANCE (km) FIGURE 6:(c) Comparison of observed and calculated travel times (see Figure 5 caption for plotting style). (d) Observed record sections scaled to common maximum amplitude.  '11  0^50^100^150^200^250^300^350 DISTANCE (km) FIGURE 7:(a) Observed record section for SP-48 (see Figure 5 caption for plotting parameters and general description). (b) Synthetic section for SP-48 (see Figure 5 caption for scaling and wavelet).  FIGURE 7:(c) Comparison of observed and calculated travel times (see Figure 5 caption for plotting style). (d) Observed record sections scaled to common maximum amplitude.  L/I \J I rlisIN-11-- ‘1\111/  FIGURE 8:(a) Observed record section for SP-49 (see Figure 5 caption for plotting parameters and general description). (b) Synthetic section for SP-49 (see Figure 5 caption for scaling and wavelet).  26  .^. FIGURE 8:(c) Comparison of observed and calculated travel times (see Figure 5 caption for plotting style). (d) Observed record sections scaled to common maximum amplitude.  27  FIGURE 9:(a) Observed record section for SP-50 (see Figure 5 caption for plotting parameters and general description). (b) Synthetic section for SP-50 (see Figure 5 caption for scaling and wavelet).  28  .^.  FIGURE 9: (c) Comparison of observed and calculated travel times (see Figure 5 caption for plotting  style). (d) Observed record sections scaled to common maximum amplitude.  29  .^. FIGURE 10:(a)Observed record section for SP-51 (see Figure 5 caption for plotting parameters and general description). (b) Synthetic section for SP-51 (see Figure 5 caption for scaling and wavelet).  30  FIGURE 10:(c) Comparison of observed and calculated travel times (see Figure 5 caption for plotting  style). (d) Observed record sections scaled to common maximum amplitude.  31  Shot point SP46 SE SP46 NW SP47 SE SP47 NW SP48 SE SP48 NW SP49 SE SP49 NW SP50 SE SP50 NW SP51 SE SP51 NW Total  Ps & Pg 1 (50) 101 (63) 48 (52) 42 (72) 64 (68) 60 (93) 73 (81) 129 (68) 65 (61) 54 (56) 108 (50) 2 (50) 749 (67)  R1  12 (144)  40 (100)  30 (88) 7 (62) 56 (92) 145 (96)  R2  R3  21 (150) 14 (111) 6 (153) 20 (95) 14 (129) 10 (127) 23 (114)  108 (124)  PmP  62 (143) 12 (108) 23 (104) 16 (78) 41 (120)  Pn  R  33 (92)  14 (93)  33 (92)  14 (93)  35 (115) 18 (120) 71 (109) 38 (129) 55 (100) 28 (96)  15 (80)  43 (97)  107 (103)  350 (114)  Total 1 (50) 243 (100) 48 (52) 105 (96) 128 (89) 174 (101) 141 (98) 235 (87) 146 (82) 61 (57) 222 (78) 2 (50) 1506 (88)  Table 2. The number of observations (travel time picks) for each phase and each shot point. Corresponding average uncertainties in ms are given in brackets. The right column gives the total number of observations (and corresponding average uncertainties) for each shot point; the bottom row lists the total number of observations (and corresponding average uncertainties) for each phase.  32  Ps is visible on all record sections, except on SP-50 NW (Figure 9a). This nearoffset portion of SP-50 NW was the least successfully modelled feature of the data, mainly because it was difficult to concurrently model the relatively prominent Ps phase on SP-50 SE. From section to section, the Ps phase shows large variations of approximately 2 to 5 km/s in apparent velocity (i.e. the inverse of the slope of the travel time-distance curve). The amplitude of Ps is invariably high, which indicates a large velocity gradient in the near-surface layer. The Pg phase appears on all sections, after the Pg-Ps cross-over distance, as the first arrival to offsets of 110 — 170 km; its average apparent velocity is 6.2 km/s. The phase is essentially linear with no distinct cross-over points, which indicates that there are no distinct layers of different velocities and no strong vertical velocity gradients in the upper crust. As well, there are no substantial travel time offsets to indicate the presence of thick low velocity zones. The fluctuations which appear within the branch were considered to be caused by near-surface features such as topography and relatively rapid changes in velocity throughout the near-surface layer. With the exception of SP-46, the amplitude pattern of Pg varies minimally from section to section, which suggests that there are no major lateral changes in the velocity gradient within the upper crust. As shown in Figure 5a, however, the amplitude of the phase is exceptionally large for SP-46, indicating a strong gradient in this region. Within each coherent Pg branch, the amplitude decreases quite abruptly at certain offsets (ranging from 30 km on SP-50 ,  SE to -100 km on SP-46 NW; see Figures 9a and 5a). If near-surface effects can be ruled ,  out, such fluctuations indicate a sudden decrease in velocity gradient at some depth within the upper crust. The Pn phase appears as a first arrival on only one section — SP-46, where the largest shot was fired. The Pg-Pn crossover point occurs at an offset of about 180 km and the apparent velocity of Pn is 8.0 km/s. This Pn phase is also essentially linear, with fluctuations being ,  33  attributed to near-surface variations as in the Pg case. The secondary arrivals from the upper crust, identified as R1, are apparent on 5 record sections, at the source-receiver offsets shown in Table 3. The R1 arrivals on SP-51 SE, SP-50 NW and SE, and SP-48 SE are virtually parallel to the Pg first arrivals, with an average delay of 0.30 s; the R1 arrivals on SP-46 are visible at offsets of 140-200 km, where the Pg arrivals are obscured by noise. The R1 phase is generally prominent — its amplitude is greater than that of the Pg phase in all cases. The amplitudes of R1 were used to constrain the velocities at the top of the middle crust. The phase labelled R2 was identified as a wide-angle reflection from the base of the middle crust and is visible on 7 sections; see Table 3 for the source-receiver offset ranges at which R2 appears. The asymptotic apparent velocity of the R2 phase is  ,  6.40 km. The  amplitude of the R2 arrivals is invariably smaller that of the R1 arrivals. The R2 phase was the principal constraint for the velocities of the middle crust and the amplitudes of the R2 arrivals were used to determine the velocity contrast across the middle crust — lower crust interface. However, it was not possible to obtain any information from the R2 amplitudes of SP-47 NW as this phase is visible on sections plotted in common maximum amplitude format (Figure 6d), but cannot be accurately discerned on true amplitude plots (Figure 6a). The secondary arrivals from the lower crust consist of two phases: R3, identified as a wide-angle reflection from a horizon within the lower crust, and PmP, identified as a reflection from the crust-mantle boundary. The R3 reflections appear on the seismic sections from 5 shot points. The source-receiver offsets at which the respective R3 branches are observed are shown in Table 3. The R3 phase typically appears 0.25-0.40 s before the PmP arrivals (except for SP-49 NW, Figure 8a, where the R3 branch is  ,  0.8 s before PmP at offsets  between 50 and 100 km). The asymptotic velocity of the R3 phase is between  ,  6.40 km/s  34  Phase^Shot Point^Source-receiver Offset R1^SP-46 NW^140 -200 km SP-48 SE^60 -110 km SP-50 NW^60 - 110 km SP-50 SE^60 -30 km SP-51 SE^80 - 150 km R2^SP-46 NW^185 - 220 km SP-47 NW^100 -140 km SP-48 NW^100 - 115 km SP-48 SE^70 - 100 km SP-49 NW^150 - 175 km SP-49 SE^75 - 100 km SP-50 SE^80 - 105 km/180 - 210 km R3^SP-47 NW^140 - 160 km SP-48 NW^120 - 170 km SP-49 NW^150 - 175 km SP-49 SE^50 - 170 km SP-51 SE^110 -135 km PmP^SP-46 NW^155 - 275 km SP-47 NW^120 - 180 km SP-48 NW^110 - 230 km SP-48 SE^80 - 115 km SP-49 NW^90 - 170 km SP-49 SE^100 - 175 km SP-50 SE^85 - 100 km/180 210 km SP-51 SE^70 - 170 km  Table 3. The source-receiver offsets at which each phase was observed is specified for the relevant shot points.  35  and -- 6.50 km/s in each case. ,  The amplitudes of the R3 phase vary: for SP-51, it is larger than R1 and quite prominent over the small distance that it is apparent; for SP-49, the R3 phase is less prominent than its SP-51 counterpart, but still larger in amplitude than the earlier R2 reflection; and the R3 phase visible on the seismic section from SP-48 is comparable in amplitude to the R2 phase, as was assumed to be the case for R3 of SP-47. The overall trend is decreasing R3 amplitudes from NW-SE. The R3 phase was used to constrain the velocities in the upper portion of the lower crust, and its amplitudes were used to constrain the velocities below the R3 reflecting horizon. The very prominent PmP arrivals are visible on 8 record sections at the offsets shown in Table 3. The asymptotic apparent velocities average 6.7 km/s. The source-receiver offset and reduced arrival time near each critical point — values which are sensitive to the Moho depth and the overlying velocity structure — vary significantly, as shown in Table 4. These variations indicate that lateral velocity changes in the base of the lower crust, topographic variations along the crust-mantle boundary, and lateral changes in the velocity contrast across the Moho are likely to be present. In terms of amplitudes, SP-48 NW is anomalous as it is not relatively large — this may be a further indication of significant lateral complexity in the Moho. The PmP arrivals were used to constrain the velocity at the base of the lower crust and the topography of the Moho, and their amplitudes were incorporated to give information about the velocity contrast across the crust-mantle boundary. The wide-angle reflection R is present on the seismic section from SP 46, starting at 315 km shot-receiver offset and 8.0 s reduced travel time (see Figures 5d and 5e). The asymptotic apparent velocity of this feature is > 8.0 km/s, and its amplitude is bigger than Pn and smaller than the other reflections. These features of the phase led it to be identified as a sub-Moho reflection.  36  Shot Point  ^  Offset of Critical Point (km)^Reduced Arrival Time of  Critical Point (s) SP-51 SE^80^ 6.9 SP-50 SE^100^ 6.9 SP-49 NW^100^ 7.2 SP-49 SE^100^ 6.5 SP-48 NW^110^ 7.0 SP-48 SE^165^ 7.0 SP-47 NW^80^ 6.8 SP-46 NW^105^ 7.6 Table 4. The source-receiver offsets and reduced arrival times of the critical point for PmP is given for each appropriate shot point.  3.1.2 Travel Time Picking Identifying and picking travel times of the phases involved looking at plots (of various scales) as well as computer screen presentations using interactive software. The procedure followed for travel time picking was dependent upon the nature of the phase, as outlined below; the preferred trace format was wiggles in all cases. First, the Ps and Pg first arrivals were picked from large-scale presentations of unfiltered traces using the interactive software package. This was a relatively straightforward task due to the high quality of the data and the uniform nature of the waveforms. The Pn first arrivals were picked in the same way. Each travel time pick was assigned an uncertainty dependent upon the difficulty of pinpointing the first motions. Ranging between 25 and 150 ms, these uncertainty values were used in the interpretation procedure (see below). The average errors for Ps, Pg, and Pn, as well as all other phases, are shown in Table 2. The picking of most of the R1 travel times was accomplished in the same manner. In regions where the onsets of the arrivals were difficult to discern in comparison with Pg, it proved useful to make initial picks on large-scale displays of the data plotted at a reducing  37  velocity of 6 km/s, and then make fine adjustments using the interactive software. The travel times for R2 were also picked in this manner for sections from all shotpoints except SP-47, for which it was necessary to use a filtered section due to high noise levels. For the R3 picking, the method used was the same as that of R1, except that the plots were displayed using an 8 km/s reducing velocity, in order to increase the coherency of this higher-velocity phase. For the PmP arrivals, the picking procedure was complicated by variations in the waveform. Plotting the sections in common maximum amplitude format (Figures 5d-10d) was helpful for the picking of the PmP phase, as this sometimes enhanced the phase coherency. Common maximum amplitude plots were also useful in the picking of the upper mantle arrival R which appeared on SP-46 NW, as seen by comparing Figure 5d with Figures 5a and 5e.  3.2 Interpretation Methods 3.2.1 The Modelling Algorithms The modelling procedure applied to the refraction data combined two approaches: 1. the iterative modelling algorithm of Zelt and Smith (1991) was used to create a model of interfaces and velocities via the inversion of travel time data; 2. the forward amplitude modelling algorithm of Zelt and Ellis (1988) was then used to adjust the model based upon the additional constraints of amplitude information. Both the travel time algorithm and the forward amplitude modelling algorithm are 2-D methods based on the numerical solution of the ray-tracing equations of zero-order asymptotic ray theory (terven5r et al. 1977), with first-order asymptotic ray theory incorporated for the calculation of head-wave amplitudes. The numerical evaluation of the travel time along a ray path uses the Runge-Kutta method (see, for example, Sheriff and Geldart 1982) with error control as suggested by terveq et al. (1977). The application of Snell's law at each point where a ray intersects a boundary was the final aspect of the ray tracing procedure.  38  An inversion method for the travel time modelling was chosen over a forward travel time modelling procedure primarily because it is a less time-consuming, more objective approach.  In addition, the Zelt and Smith (1991) routine provides useful estimates for the resolution of model parameters (i.e. boundary depths and layer velocities). Furthermore, the chosen inversion program was written for crustal seismic data with a shot-receiver density too low to fulfill the requirements of conventional tomographic or full wavefield techniques. In this sense, the algorithm is ideally suited to the Line 10 data. The Zelt and Ellis (1988) algorithm was chosen because it efficiently allows rapid forward modelling of amplitudes; although more accurate methods (e.g. finite-difference algorithms) exist, they were not deemed computationally practical for this study. An important aspect of this algorithm is that in order to calculate synthetic seismograms, a value of the physical constant Q (the "quality factor") must be assumed in order to account for attenuation, the reduction in amplitude or energy caused by the characteristics of the transmitting media. The values of Q, which is defined as 27/(fraction of energy lost per cycle) (Sheriff and Geldart 1982), used to calculate the amplitudes of the synthetic seismograms were not determined through the modelling procedure, but were adapted from the values of Hough and Anderson (1988) determined for the Anza, California region. That study area was very similar to the CB site, as it was situated within the Southern California Batholith. The dominant geology (granitic plutons), age (Jura-Cretaceous), and tectonic origin (subduction off the west coast of North America) of the Southern California Batholith make it a suitable analog for this comparison with the CB. The Q structure used in modelling is shown in Table 5. One convenience of using the travel time inversion and the forward modelling algorithms concurrently was that the model parameterizations for the two methods were identical. The parameterization procedure is based upon the positioning of two types of nodes: boundary and velocity nodes. The boundary nodes are essentially depth nodes, as their horizontal position is user-specified, with the inversion program accounting for the adjustment of the  39  Depth From Surface^ Q (dimensionless) 100 1.5 km^ 300 8.0 km^ 1000 22.0 km^ 700 72.0 km^  Table 5. The Q structure used in modelling is shown in terms of depth from the surface. Notes: each Q represents the value assigned to a homogeneous Q layer with its bottom at the corresponding depth and its top at the previously-mentioned depth. Values were adapted from Hough and Anderson (1988). depth. The boundary nodes are used to construct the model in terms of layers; a layer boundary is represented by an arbitrary number and spacing of boundary nodes which are connected by linear interpolation. These horizons must cross the model from left to right without crossing another boundary. The velocity structure within each layer is represented by a series of velocity nodes at the top and at the bottom of each layer. Again the number and spacing of nodes is arbitrary. The velocity between the upper velocity nodes in a layer is determined through horizontal linear interpolation, as is the velocity between the lower velocity nodes. The final picture of the layer velocity is achieved by linearly interpolating vertically between the upper and lower velocities in each layer (see Zelt and Smith 1991). The resulting parameterized model is made up of layers composed of different-sized trapezoids which are divided by vertical segments that are included wherever there is a boundary or velocity node. These vertical segments within each layer do not represent velocity discontinuities — the velocity is laterally continuous within each layer. The quasihorizontal boundaries between layers can represent velocity discontinuities if required; they can also represent second-order discontinuities (i.e. a change in velocity gradient), or pseudoboundaries used for ease of modelling. A smoothing procedure is applied to reduce the effects  40  (i.e. scattering and focusing of ray paths and geometrical shadows) of this blocky model parameterization (Zelt and Smith 1991). Since the structure of the final model is strongly dependent upon the nature of the inversion algorithm, the aspects of the inversion procedure that are important in this context are discussed below. This derivation of the algorithm's essential equations is an alternative approach to the one taken by Zelt and Smith (1991), who solved equation (1) using a damped least-squares technique presented by Aid and Richards (1980), following the method described by Lutter et al. (1990).  3.2.2 Inversion mathematics If a velocity field is considered in terms of its discrete form instead of its form as a continuous velocity field v  (x,  y), the travel time t along a raypath between a source and  a receiver is t=  E li/vi  where vi is the velocity of the i th ray segment, li is the length of the i th ray segment, and n is the number of segments comprising the ray path (i.e. the number of cells the ray samples in passing through the velocity field). This travel time expression is linear in terms of slowness (reciprocal velocity), but non-linear in terms of raypath since the path of the ray depends upon the velocity structure it encounters. The following linearized expression for travel time residual (the difference between observed and calculated travel time) is obtained by expanding the velocity field about a starting velocity model in a Taylor's series expansion and neglecting the higher order terms: At = AAm^  (1)  where At is the travel time residual vector, Am is the model parameter adjustment vector,  and A is a matrix of partial derivatives. A contains the elements  ati/ami,  where ti is the  41 i th observed travel time, and mi is the j th model parameter selected for inversion (i.e. either  a velocity value or the z-coordinate of a boundary node). To solve equation (1), rays are traced through the current model during a particular iteration and At and A are computed, with partial derivatives calculated analytically using the method described in Zelt and Smith (1991). This determines the model update vector Am, which is added to the current model for further iterations or for stopping if the stopping criteria are met. Crustal refraction experiments often give rise to a system (1) which is overdetermined, underconstrained, and ill-conditioned. For such systems, the damped least squares method of solution is a common choice (e.g. Braile 1973, Kanasewich and Chiu 1985, Chiu and Stewart 1987, Phadke and Kanasewich 1990). In solving the linearized data equations (1), Zelt and Smith (1991) considered the overdetermined case by choosing to solve for fewer model parameters than travel time observations, and they sought a model update vector to satisfy the system of equations in the least-squares sense via a damped least-squares solution. The relevant objective function needs to be of the form: 0 (Am) = II AAm  —  Atir I 6 2 IIAmir - -  where 6 2 is a scalar trade-off parameter that is specified in order to set the relevant importance of the terms, and 11.II means the 12 norm. Extremizing 0 with respect to Am yields the linear algebraic equations Am= (ATA Fe 2 I) ATAt^ -  (2)  Following the method suggested by Lutter et al. (1990), if the trade-off parameter is defined as E 2 = Do- ? /0- 2 where = the standard deviation used to estimate the uncertainty of the i th travel time measurement, 0- • = the standard deviation used as an a priori estimate of the uncertainty of the i th model parameter,  42  D = an overall damping parameter usually equal to one,  then this becomes the Zelt and Smith (1991) damped least-squares solution to (1): Am= (A T C  1 A-1-DC;„;1 ) -1  ATC  l At^  (3)  where C t = the estimated data covariance matrix defined by diag (al ), and C m = the estimated model covariance matrix defined by diag fq  ).  The values of the standard deviation Cri were specified during travel time picking by assigning an error estimate to each observation, as mentioned in the previous section. The choice of the values used for the standard deviation a • will be included in the forthcoming discussion of the modelling procedure. It is important to observe that by removing some of the underdeterminancy of the problem by underparameterizing, the inverse solution becomes more stable, but also more dependent on the chosen model parameterization. It is therefore of utmost importance that a model parameterization appropriate to the subsurface geology is chosen (i.e. within the capacity of blocky trapezoids to represent the subsurface). The optimal parameterization uses the maximum number of nodes that the data are able to constrain, for if these nodes are prudently positioned they promise enhanced knowledge about the subsurface. Zelt and Smith (1991) suggested that the optimal node placement may be determined by running a series of tests based on the following criteria: 1. the travel time residual must be sufficiently small; 2. parameter resolution (as determined by the resolution matrix R) must be adequately high; and 3. rays must be traced to all observations.  43  Using these three criteria, a preliminary series of inversions was performed to attempt to optimize the number and placement of the velocity and boundary nodes used in inversion. In practice, this was possible for the boundary nodes; for the velocity nodes, however, while the number of nodes chosen could be optimized, in certain cases their placement could not. In these instances, which are discussed later, the choice of velocity node positioning as governed by the above guidelines was somewhat arbitrary. The resolution matrix R defined by Zelt and Smith (1991), R. (A  T Ci-1 A+DC;,I) -1  A T CT 1 A  (derived in Lutter et al. 1990), has as its diagonal elements numbers that range between values of 0 to 1. Non-zero diagonal elements indicate the degree of averaging or the linear dependence of the true model space as represented by the inverted model. Model parameters associated with diagonal elements of the resolution matrix that are ?_ 0.5 are considered meaningful and well-resolved (see Zelt and Smith 1991). The criteria determining optimal parameterization lead directly to the stopping criteria: 1. the final model chosen is that which gives the preferred trade-off between travel time residual and parameter resolution; and 2. the final model must permit rays to be traced to all observation points. The crux of (1) is that adding more nodes to a model typically allows the travel time residual to be decreased, but at the same time decreases the overall parameter resolution. It is essential to realize that the R matrix gives the resolution values for the parameters of the model update vector Am; this does not provide knowledge about the resolution of the parameters of the final model. In terms of tracing rays to all observation points, it is often not possible using a model obtained from an inversion. Zelt and Smith (1991) cited the cause of this as using more parameters than the data can adequately constrain. This may lead to shadow zones because  44  of the potential for unrealistic oscillations in layer velocities or boundary shapes, which cause rays to focus and scatter. A final model must therefore be rejected if it is not possible to trace rays to most observation points. Note that the ray take-off angles that connect sources and observed receiver locations are determined so that ray tracing is performed in the "two-point ray tracing" sense, whereby rays are traced only between sources and observation points. 3.2.3 The Modelling Techniques The approach used for modelling the data with the two algorithms was one of "layer stripping", by which all the observed travel times associated with a given layer are simultaneously inverted while the values of the parameters of the layers above are fixed. The model thus determined with the Zelt and Smith (1991) inversion algorithm was adjusted, if necessary, for the added constraints of amplitudes using the Zelt and Ellis (1988) forward modelling procedure. This process of alternating travel time inversion and amplitude forward modelling was continued until an acceptable fit of observed and calculated traveltimes and amplitudes was achieved for successively deeper layers. In this way, the final model was constructed. This method clearly simplifies and expedites the inversion process. Testing has shown that the method is valid, since, for example, PmP and Pn arrivals provide little constraint on upper and middle crustal structure (Zelt and Smith 1991). In practice, a trade-off between reducing the travel time residual and improving the computed amplitudes is inevitable. More emphasis was placed upon achieving a satisfactory travel time fit, because of the limitations involved in calculating amplitudes using ray theory. Due to the large amplitude variations within each phase (which is a common trait of data from refraction surveys on land), the comparison between computed and observed amplitudes was qualitative. In applying the modelling techniques in this manner, two assumptions were made: (1) the inversion algorithm is capable of resolving trade-offs between velocity and boundary depth; and (2) the effects of using a 2—D (straight-line) representation of the actual shot and receiver geometry are negligible. Zelt and Smith (1991) asserted that testing has revealed  45  assumption (1) to be valid. As well, assumption (2) seems valid in the case of Line 10, with the possible exception of the vicinity to the southeast of SP-48 (Figures 2 and 4).  3.2.4 Modelling the Refraction Data This section describes the application of the modelling procedure to the Line 10 data set in order to construct the final model (Figure 11a). Figures 5c-10c show the match between observed and calculated travel times, and the synthetic sections shown in Figures 5b-10b allow calculated and observed amplitudes to be compared. The location of the boundary and velocity nodes used in the inversion are shown in Figure 12, which also shows the numbers which were assigned to the boundaries and layers in the text. Figure 13 shows the ray coverage within (a) the upper crust and (b) the remainder of the crust and the upper mantle. All boundaries — including the surface, which represents the topography — in an inversion model are smoothed by uniform sampling at 250 points and the application of a three-point averaging filter (Zelt and Ellis 1988). The sub-horizontal boundaries shown by dashed lines in Figure 1 la represent velocity discontinuities (those which delineate changes in vertical velocity gradient have been omitted in the presentation). Since the velocity structure of the model suggested a natural division of the crust into thirds, the upper, middle, and lower crust were defined for the discussion. Table 6 provides a summary of the results of the inversion procedure.  3.2.4a The Upper Crust The uppermost layer in Figure 1 la represents a thin (P-2.5 km) near-surface stratum. Velocities within the layer were determined by inverting the Ps arrivals. The starting model had 6 velocity nodes on the surface at the shot point locations, which preliminary modelling had shown to be the optimal node locations; indeed these were the only zones of ray coverage. Velocity nodes representing values at the bottom of the surface  46  FIGURE 11: (a) Color contour plot of the velocities beneath Line 10 (vertical exaggeration approximately 4:1). Numbers represent average P-wave velocities in km/s between the dashed lines (not nodal values). The distance axis has SP-51 as the origin ("model distance"), and the depth axis has sea level as the origin. The shotpoint locations are shown by solid triangles, and positions of refraction Line 2 of SCORE '89 and LITHOPROBE reflection profiles 88-13, 88-14, and 88-17 are indicated. Thick solid lines represent areas which are constrained by reflected energy. The transitional Moho layer is indicated by "M". The shaded region around the model is unconstrained by the travel time data (The upper mantle reflector has been omitted here.) FIGURE 11: (b) A simplified schematic interpretation of (a) based partly on other studies (see text). Vertical exaggeration is approximately 2:1. Note that the interpreted extent of the CB batholiths does not mle out their presence below this depth.  NW  DISTANCE (km) 200 150 50^1001 1^ ^  0  ,,-----------:7 1---0-  0  ^-  0-  -  11  i  13---v-------a---5  °.  1^•  250  (  300  350  SE  •  5 6  6 0— C\2  7 8^el ^n --- v  0o--  0  N0  10  -  -  11  FIGURE 12: Location of depth (squares) and velocity (circles) nodes used in the inversion process. The numbers in the left column represent layer identification numbers used in the text, and those in the right column denote boundary numbers. There are no nodes shown on boundaries 3-5 because they were defined through forward modelling of amplitudes (see text). Boundary 1 represents the topography of Line 10, layer 1 is the near-surface layer, layers 2-5 constitute the upper crust, layer 6 is the middle crust, layers 7 and 8 make up the lower crust, layer 9 is the transitional Moho layer, and layers 10 and 11 are upper mantle layers separated by the upper mantle reflector shown as boundary 11. The shotpoint locations are shown by solid triangles. Vertical exaggeration is -3:1.  48 DISTANCE (km) 2 150^200  (a) 0  50  100  1  DISTANCE (km) 150^200  250  300  350  0  0  0 02 0  0 •TV 0  0  0 N  (b) FIGURE 13: Two-point ray tracing diagrams (all shots to all receivers) showing ray coverage in (a) the upper crust, and (b) the remainder of the crust and the upper mantle. The vertical exaggeration in (a) is approximately 9:1, and in (b) approximately 3:1. The shotpoint locations are shown by solid triangles.  49  Phase Inverted Ps Pg , layer 2 Pg, layer 3 Pg, layer 4 R1 R2 R3 PmP Pn R  RMS Travel Time Residual .01 .01 .12 .12 .06 .09 .11 .11 .09 .03  Normalized Chi-squared Misfit 7.7 4.5 1.9 2.6 0.7 1.0 1.3 1.4 1.1 0.2  Number of Velocity Nodes 11 12 0 0 0 6 6 6 1 0  Number of Depth Nodes 0 0 0 0 3 4 3 3 0 1  Table 6: Inversion results for the final velocity model. The RMS travel time residuals, chisquared misfits, and number of velocity and boundary nodes are displayed for each phase of the inversion procedure. Note that the RMS travel time residual is given in seconds. See Figure 12 for the placement of the boundary and velocity nodes. layer were initially placed 1 km beneath the shot point locations, a depth based upon the estimated thickness of the surface layer. The damping factor for the inversion was set at 1. A single value of a i (a,), 0.1 km/s was used for all velocities and a single value of  ai (a z ),  0.1 km was used for all boundary  nodes. These parameters were unchanged for each step of the layer-stripping procedure, except for the cases which are described below. The relative values of  ay  and a z determine  the trade-off between the size of velocity and boundary adjustments in the inversion. Since the Ps arrivals were limited in both number and spatial extent, extracting velocity gradient information was not feasible. The velocity gradient in this surface layer was therefore fixed at approximately 0.5 s -1 for the inversion procedure, based on a priori estimates of the velocities at the top of and at the base of the layer.  50  The resulting model had an acceptable travel time fit (i.e. the RMS travel time residual of 0.01 s was of the same order as the uncertainty of the travel time picks), and rays were traced to all observations. However, the resolution of the velocity nodes — 0.57, 0.84, 0.47, 0.42, 0.43, and 0.12 — was relatively low, a reflection of the fact that the shot-receiver geometry was inappropriate for obtaining extensive near-surface information. To determine the base of the Ps layer, and to constrain velocities in the upper crust, the 702 Pg arrivals were inverted. The modelling of the Pg arrivals was initiated by placing 6 boundary nodes beneath each shot point at the Ps-Pg interface (boundary 2 on Figure 12). Six velocity nodes were also placed at these locations, and six additional velocity nodes were added at the necessary intervals to invert for the velocity at the top of the Pg layer. This approach was adopted after initial modelling revealed that an acceptable travel time fit could not be obtained with fewer velocity nodes. Through forward modelling it was determined that the deepest turning point of the Pg first arrivals was 8 km, so lower velocity nodes were placed at this depth beneath the 12 upper nodes; this estimated depth of penetration was based solely upon the extent to which the Pg phase was visible as a first arrival on the field data. The boundary which appears at 8 km depth in Figure 12 (Boundary 5) therefore does not represent a structural feature; it indicates the depth to which velocity calculations in the crust were controlled by refracting rays. Preliminary inversion results indicated that inverting for the base of the Ps layer was not going to be successful — the boundary node positions oscillated with successive iterations instead of converging to a solution, and there were wild undulations in the boundary which led to shadow zones in ray coverage. As a result, the Ps-Pg boundary was determined through forward modelling of the Pg arrivals. This procedure required that 6 extra boundary nodes be placed along the Ps-Pg interface. Note that these boundary nodes (12 in total) do not appear in Figure 12, as they were determined through forward modelling alone. Velocities to 8 km depth were then determined by carrying out the inversion procedure.  51  During the course of the inversion procedure, some of the idiosyncrasies of the inversion program mentioned by Zelt and Smith (1991) became evident. The velocity gradient beneath SP-49 became sufficiently high that the Pg arrivals observed at large offsets for shot point 46 were not getting rays traced to them. This problem was solved when the velocity gradient below SP-49 at the bottom of the layer was held fixed at a more reasonable value for a portion of the inversion procedure. Also, the velocity below SP 51 was unreasonably large, so it was necessary to hold this value fixed during some iterations. For the final iterations, the damping factor was increased to as high as 25 in order to come up with the model that had the best travel time fit and had rays traced to the maximum number of observation points. While the inversion for the velocity in the Pg layer was taking place, it was necessary to continue the modelling of the Ps layer in order to account for "kinks" in the essentially linear Pg branch. It was apparent that these small undulations were caused by geological features, as their positions did not change on the six different record sections; as such, they were modelled as zones of higher or lower velocity material in the near surface layer. The resulting Ps layer was defined on its top and bottom by 11 velocity nodes; the bottom nodes do not appear in Figure 12 because the gradient in the overlying material was fixed, so the values attached to these nodes bottom depend upon values determined for the upper nodes through inversion. One of the "kinks" was relatively large — it appears at 50-60 km model distance on the Pg branches of SP's 49-51 (Figures 8a-10a). This feature corresponds with eight receiver stations in the northwest which traversed the Cadwallader terrane (Figure 4), a Coast Belt terrane of island-arc clastics and volcanics which may be linked with Stikinia (Wheeler et al. 1991); it was successfully modelled as a localized high velocity (5.0-5.6 km/s) zone. This location represents one of two times that the refraction line crossed a sliver of Cadwallader terrane; however, a similar kink in the data does not occur at the second crossing. Although there is no obvious explanation for this feature of the data, the fact that the second crossing  52  occurs where Line 10 deviated from a straight-line configuration may be significant. The velocity model determined down to 8 km in this manner was then used as the starting model for the amplitude forward modelling procedure. The initial model was quite problematic in terms of its amplitudes — there was a general trend on the synthetic seismograms for the amplitudes to increase with offset which was not seen on the seismic sections plotted from the field data. This indicated that the rays which penetrated more deeply should be passing through material with a lower gradient than the more shallow-diving rays. In order to accommodate this reduction in gradient with depth, two new boundaries were added to the Pg layer, at approximately 2 km and 4 km depth (see boundaries 3 and 4, Figure 12). These boundaries, across which there is no velocity discontinuity, delineate a change in velocity gradient. The gradients were then adjusted until the values gave a satisfactory fit to the amplitude data. Each time a velocity node value was changed in order to change the velocity gradient, forward travel time modelling was done to ascertain that there was no significant deterioration of the travel time fit. The amplitude fit was acceptable when the second, third, and fourth layers had vertical velocity gradients of about 0.027 s -1 , 0.009 s -1 , and 0.003 s -1 , respectively. The travel time data were then inverted again to solve for the velocities at the top of layer 2 with the velocity gradients fixed at these values. The resolution of the 12 velocity nodes in layer 2 is good (> 0.71) in all cases. Keeping the velocities in layers 2 — 4 fixed at their final values, the Ps and Pg arrivals were inverted again for all 11 of the upper velocity nodes in layer 1. The resulting RMS travel time residual is 0.01 s (Table 6). The resolution values are generally acceptable at the shot point locations (values were between 0.48 and 0.84), while resolution values at the intervening nodes range from 0.18 to 0.52. The low resolution values of certain nodes indicates that the model depends strongly on the parameterization. The modelling of the upper crust was completed by inverting for the R1 reflection. Tests  53  showed that the velocity in layer 5 (Figure 12) is not constrained by the data, so the velocity in this layer was modelled by fixing the vertical velocity gradient at 0.003 s  -1  as in layer  4, and having no velocity discontinuity between layers 4 and 5. Boundary 5 in Figure 12 is therefore a pseudo-boundary, as the velocities in the layer below are a continuation of those calculated for the layer above. Three boundary nodes were placed at a depth of 10 km, and were inverted to determine the depth to the reflecting horizon; the number and position of the boundary nodes were determined to be optimal through preliminary modelling. The inversion algorithm gave a model with the R1 reflector having an average depth of 10.3 km; the reflector is flat throughout most of the model, with an upward dip at the northwestern end. The three boundary nodes used in inversion have excellent resolution (0.94, 0.96, and 0.94), and the RMS travel time residual value is 0.06 s for the R1 arrivals (Table 6). In the case of R1, the boundary was extrapolated between the portions which are constrained by the reflections (see Figure 11a). This procedure was adopted for all boundaries determined through the inversion of reflected phases. The final upper crustal model has 23 velocity nodes and 3 depth nodes. The travel time fit was 0.096 when rays were traced for all the phases. Rays could not be traced to 10 of the 702 Pg data points for the reasons described earlier, and the corresponding travel times were therefore not used in the inversion.  3.2.4b The Middle Crust The mid-crustal structure of the model was determined by  inverting the R2 reflection. Preliminary modelling showed that the velocity gradient within layer 6 is not constrained by the data, so the velocity gradient was again fixed at 0.003 s -1 as in layers 4 and 5. This procedure was deemed reasonable because there was no evidence in the data for strong gradients in this layer, and it was consistent with similar studies in the Sierra Nevada (eg. Bolt and Gutdeutsch 1982, Pakiser and Brune 1980).  54  To invert this mid-crustal reflection, 4 boundary nodes were placed at a depth of 20 km. Velocity nodes were placed at the top of layer 6 (boundary 6 in Figure 12) beneath four of the shot points, to the exclusion of the two on the ends of the profile. Preliminary modelling had indicated that this was the optimal number of nodes, as inverting with less nodes led to a model which did not accurately fit the data, and using more nodes significantly decreased the resolution. The location of the nodes was more arbitrary, since the positioning criteria could be met equally by a number of placements, as long as the nodes were roughly evenly-spaced across the region constrained by the data. The inversion produced an initial model which had the R2 reflector at an average depth of 22.2 km, except an upwelling which reached  ,  19 km depth beneath SP 48, where the result  is constrained by data from both SP-46 and SP-47. For this initial model, the velocities were 6.30 km/s beneath SP's 50, 49 and 47, and 6.55 km/s beneath SP 48; the portions of layer 6 below SP's 51 and 46 are not constrained by the data, as indicated above. With the modelling of layer 6 (Figure 12) completed through the inversion of R2, the velocity contrast across boundary 6 was specified, so that the forward modelling of the R1 reflections was possible. This starting model was problematic — the R1 reflections from the northwestern part of the line were too low in amplitude, and those from the southeastern part had traces in which the amplitudes were too high. The velocities in layer 6 were adjusted until the R1 amplitudes on the synthetic seismograms matched those on the original data; holding these velocities fixed, the R2 arrivals were inverted to correctly position the boundary nodes at the base of layer 6. During the course of the R1 amplitude modelling, it was necessary to add a velocity node to the northwest of the one beneath SP-50. The resulting R2 reflector has an average depth of 22.8 km, except for the upwelling beneath SP-48, which reaches a depth of 19.7 km (Figure 11a). The velocities determined for the layer are --6.35 km/s beneath SP-50 and SP-49, 6.45 km/s beneath SP-48, and 6.40 ,  km/s beneath SP-47. The velocities beneath SP-51 and SP-46 are not constrained due to  55  inadequate data. The RMS travel time residual is 0.09 s for the R2 arrivals (Table 6), and the resolution of the boundary nodes is very good: 0.97,0.97,0.94, and 0.70. The resolution of the 4 velocity nodes used in the inversion is also good (0.71 or greater in all cases). 3.2.4c The Lower Crust and Upper Mantle The modelling of the lower crustal velocity  structure was initiated by inverting the R3 arrivals. The vertical velocity gradient was fixed as for the R2 phase. Three boundary nodes were placed at a depth of 30 km, and velocity nodes were placed beneath four of the shot point locations at the top of layer 7 (boundary 7 in Figure 12). Preliminary modelling was used to ascertain that this was the optimal number of nodes in terms of the criteria described above; again, the placement of the nodes was somewhat arbitrary, as described for the boundary 6 case. The initial model which resulted from the inversion procedure had velocities in layer 7 of 6.35, 6.40, 6.50, and 6.45 km/s beneath SP's 50, 49, 48 and 47, respectively, and had a flat R3 reflecting boundary with an average depth of 29.9 km. The portions of the model beneath SP's 51 and 46 are not constrained by the data. When the model based upon travel time inversion was used in the forward amplitude modelling procedure, the amplitudes on the synthetic sections were seen to be generally too small for the R2 reflections in comparison with those on the true data. The exceptions were the R2 reflections from SP's 46 and 47, which compared well in terms of amplitudes. The velocities in layer 7 were adjusted until the amplitudes of the R2 reflections on the synthetic seismic sections matched those on the real data sections; with these velocities held fixed, the R3 arrivals were inverted to model the base of the layer. The inversion model which produced an acceptable travel time fit (see Table 6) has velocities of 6.45 km/s beneath SP-50, 6.50 km/s beneath SP-49, and SP-48, and 6.45 ,  km/s beneath SP-47 (see Figure 11a). The velocities at the nodes beneath SP's 51 and 46 were arbitrarily fixed at the values of the adjacent velocity nodes. The R3 reflecting boundary is essentially flat at a depth of 30.0 km. The resolution values of the boundary nodes are  56  excellent (0.98,0.98,0.97), as are those of 3 of the 4 velocity nodes used in inversion (> 0.83). For the velocity node beneath SP-47, the resolution is poor (0.32) because the ray coverage was not high in the vicinity. However, added confidence may be ascribed to this value because the R2 reflection on the southwestern portion of SP-48 boundary 7 is constrained at the location of its node, and the amplitude of the reflection therefore provides additional constraint. Since both the PmP and Pn arrivals constrain the structure of the Moho, these arrivals were inverted simultaneously to model boundary 9 on Figure 12. According to the criteria described above for the R2 and R3 modelling procedures, the velocity gradient was fixed, 3 boundary nodes were located at the base of layer 8, and velocity nodes were placed beneath four of the shot point locations at the top of layer 8 (Figure 12). Since the Pn data were so limited (the upper mantle was sampled only in the upper kilometer, and the lateral extent of this was confined to within 150 and 300 km model distance, as shown in Figure 13b), it was necessary to fix the vertical velocity gradient in the upper mantle at a reasonable estimated value before performing the inversion. Amplitude modelling tests led to the selection of 0.003 s -1 — a gradient of 0.004 s -1 produced amplitudes of the Pn arrivals on synthetic seismic sections which were much larger than those on the observed sections, and the choice of 0.002 s -1 caused shadow zones in ray coverage. The velocity at the top of the upper mantle is represented by a single velocity node (i.e. this was modelled as a laterally homogeneous layer, since the information available was so limited). Preliminary modelling showed that inversion for the crust-mantle boundary resulted in so much relief that rays could not be traced to all sites for the Pn arrivals; the damping factor was increased to 10 to alleviate this problem. Inverting the PmP and Pn arrivals in this manner led to a preliminary model in which the velocity values above the Moho (layer 8 in Figure 12) ranged from 6.90 km/s beneath SP-50 to 6.45 km/s beneath SP-47. The crust-mantle boundary had an average depth of  57 34.5  km, except at the northwestern end of the line, where it sloped upward.  Forward modelling of R3 amplitudes led to the following changes for layer 8: the velocity node below SP-50 was moved roughly 25 km closer to SP-49, and the velocity at this node was decreased to 6.70 km/s, and the velocity at the nodes beneath SP's 49 and 47 were increased to 6.65 km/s and 6.60 km/s, respectively. The velocity of 6.70 km/s at the node beneath SP-48 was not changed. The velocities of the nodes centered under SP's 51 and 46 were arbitrarily fixed at the values of the adjacent velocity nodes. This overall decrease in velocity contrast from NW to SE means that the observed trend of decreasing R3 amplitude from NW to SE is replicated in the calculated data (compare Figures 5b-10b with Figures 5a-10a). With these velocities held fixed, the PmP and Pn arrivals were again simultaneously inverted for the boundary nodes at the base of layer 8 and the upper mantle velocity node. The initial model of the crust-mantle boundary which resulted from the inversion process was essentially flat with an average depth of 34.6, except for an upward-sloping section at the northwestern portion of the line, which reached a depth of 32.4 km at the northwestern extent of the portion of the boundary constrained by the data. The velocity at the top of the upper mantle was 8.05 km/s. The RMS travel time residual for the PmP and Pn phases was 0.110 with this model, the resolution values of the velocity nodes in layer 8 were acceptable (0.49-0.87), and those of the crust-mantle boundary nodes were excellent (> 0.97). When this model was used in the amplitude forward modelling program, the amplitudes of the PmP arrivals on the synthetic seismic sections seemed quite large in comparison with those on the true data sections. This situation is most readily explained by one of the following: 1. the velocity of the upper mantle is too high, 2. the Moho cannot be modelled simply as a single-order velocity discontinuity. In order to see if the first explanation was feasible, tests were performed in which the velocity of the upper mantle was fixed at every reasonable value lower than 8.05 km/s, i.e.  58  7.6 km/s, 7.65 km/s, etc... These tests showed that this lowering of upper mantle velocity did not sufficiently decrease the amplitude of PmP. The second explanation is consistent with amplitude modelling in previous studies (e.g. Benz et al. 1990, Zelt and Smith 1991, and Zelt et al. 1992a), which showed that the similarity between observed and calculated data increased when the Moho was modelled as a transition zone rather than a velocity discontinuity. The method used in these studies was to incorporate a transitional Moho layer with velocities ranging from 7.4-7.8 km/s; PmP and Pn were then modelled respectively as reflections from the top of, and refractions below, this transitional Moho, which ranged in thickness from 0.5-3.5 km (Zelt et al. 1992a) to 2.0-5.0 km (Benz et al. 1990). There were two other aspects of the data which suggested that the modelling should be done in this manner. First, the calculated critical points were  ,  10-20 km closer to the  source than the observed critical points. As was the case in amplitude modelling, decreasing the upper mantle velocities did not alleviate the problem. Also, the PmP phase of SP-48 NW had a significantly lower amplitude than the others. This amplitude anomaly could be modelled with lateral velocity variation in either upper mantle if the Moho was modelled as a single-order velocity discontinuity, or within the transitional Moho itself; it was more difficult to justify the former using this one feature of the data as a basis. The velocities in the transitional Moho (layer 9 in Figure 12) were determined through the forward modelling of PmP amplitudes. Velocity nodes where placed at the top of the Moho layer beneath each of the shot point locations for the forward modelling procedure, and the velocity node below SP-50 was assigned a value of 7.8 km/s, below SP's 49 and 48 a value of 7.4 km/s, and below SP 47 a value of 7.5 km/s. The nodes beneath SP's 51 and 46 where the data do not constrain the amplitudes were given the velocity values of the adjacent nodes. Since the Pn data were so limited (33 observed arrivals were picked), it was not possible  59  to model the base of the transition zone. It was therefore arbitrarily decided that the base of the layer (boundary 10 in Figure 12) would parallel the top boundary of the layer in shape, for a transitional Moho of uniform thickness (chosen to be 2 km in Figures 11 and 12). With layer 9 defined in this manner, the Pn arrivals were inverted to determine the velocity in the upper mantle; the resulting velocity value was 8.15 km/s. Obviously, this value depends upon the arbitrarily chosen thickness of the transition zone. The extent of this dependence was shown by a simple test in which the thickness of layer 9 was decreased to 1 km: inverting the Pn arrivals led to an upper mantle velocity of 8.10 km/s. This modelling of the crust—mantle boundary as a transitional layer improves the calculated amplitudes of all PmP arrivals, puts the critical point of all the arrivals into the appropriate position, and provides a reasonable means of incorporating the PmP amplitude anomaly of SP-48 NW into the model. The transitional Moho is therefore included in the final model, although its thickness is not constrained. This may provide a more realistic picture of the true nature of the crust-mantle boundary, which has been interpreted by some (see Mooney and Brocher 1987) as a zone of inter-layered materials of higher (ultramafic) and lower velocities (mafic). The upper mantle velocity structure was completed by inverting the arrivals from the sub-Moho reflection R visible on the seismic section from SP-46. One boundary node was used in the inversion procedure, which resulted in a sub-Moho reflector at 70.3 km depth in the final model. The resolution of the one boundary node used is very good (0.88). The depth of this reflection depends upon the velocity chosen for the upper mantle. If, for example, a 1 km thick transition zone is chosen, the upper mantle velocity is 8.10 km/s, and the depth to the R reflector is 68.3 km. When amplitude modelling was performed for the R reflection, it was shown that a small velocity increase (-Ai km/s) resulted in amplitudes on the synthetic seismograms which are comparable to those on the true data sections. Similarly, a negative velocity discontinuity  60  ( 0.2 km/s) provides a good fit to the observed data (see Figure 5e cf. Figure 5f). Since ,  the relative phase (polarity) of the arrival could not be distinguished, the data do not enable a preferred velocity change to be inferred.  61  IV. INTERPRETATION 4.1 The Final Model 4.1.1 Primary Features of the Velocity Model The principal features of the velocity model are shown in Figure 1 la. The uppermost stratum is a thin (< 2.5 km) near-surface layer with velocities averaging 3.9 km/s at the top and 4.9 km/s at the bottom. Within this average velocity structure, there are zones of higher and lower velocity responsible for fluctuations in the essentially linear Pg branch. The most significant of these is a high velocity (5.0-5.6 km/s) zone which is located below SP 50, and has been linked with a feature of the surface geology (the Cadwallader Terrane). The velocities at the top of the upper crust average 6.15 km/s, and those at the bottom 6.25 km/s, with higher velocities occurring in the southeast. The depth to the lower boundary of the upper crust (^40 km) is constrained for the northwestern third and for an  ,  30 km  segment at the southeastern part of the model by the R1 phase. In the northwest, the R1 reflection is from a relatively large velocity contrast of —6.15-6.20 to 6.35 km/s. The base of the upper crust in the center of the model where the velocity contrast is intermediate (-0.10 km/s) is not constrained by observations of reflected energy, while the R1 reflection in the southeast occurs where the velocity contrast between the upper and the middle crust has decreased significantly to 40.05 km/s (Figure 11a). This suggests that another source of reflectivity (e.g. fine-scaled layering which cannot be resolved in this experiment) other than simply a relatively small single-order velocity discontinuity may have a role in generating R1 in the southeast. Alternatively, the situation may be an artifact of the modelling procedure (i.e. the blocky model parameterization). The velocities in the middle crust average 6.35-6.45 km/s, with slightly higher values below SP 48, for an overall velocity structure of little lateral variation. The depth to the R2  62  reflecting horizon defines the base of the middle crust at 22 km; this boundary is well,  constrained for the center portion of the model (i.e. 100-250 km model distance). As in the upper crust, the two northwestern reflected phases are generated by a velocity discontinuity of at least 0.15 km/s, and the southeastern reflection by one of 0.05 km/s. Since the base ,  of the middle crust is not similarly constrained by observed R2 reflections in areas with an intermediate velocity contrast ( 0.10 km/s), it is again suggested that perhaps other sources ,  of reflectivity may be in effect. The lower crust, which has the Moho as its base, is divided into an upper and a lower layer by the R3 reflector. The lower crust is also homogeneous in its velocity structure within the resolution of the survey, with values averaging 6.45-6.50 km/s in the upper portion and 6.60-6.70 km/s in the lower layer. The depth to the top of the transitional Moho averages 34.5 kin, dips to 35.5 km at the southeastern end of the model, and shows minimal topography. The velocities within the transitional layer vary from 7.80 km/s in the northwest to 7.40 and 7.50 km/s in the center ,  and southeastern portions of the line, respectively. The top of the crust-mantle boundary is well-constrained by PmP reflections from 50-300 km (model distance), but its thickness is not constrained by the data. Modelling the Moho as a transitional layer increases the similarity between the observed and the calculated data; as well, if the Moho transition is explained as a complex series of interlayered high and low velocity materials, this presentation is consistent with recent studies of the character of the crust-mantle boundary (e.g. Hale and Thompson 1982, Braile and Chiang 1986). While the sparsity of the Pn arrivals and the unknown thickness of the Moho rule out constraining the upper mantle velocity, the available data suggest that velocities are relatively high (e.g. 8.15 km/s for a 2 km thickness, r-8.10 km/s for a 1 km thickness of the Moho ,  transitional layer) in this region. An upper mantle reflector which is not shown in Figure 1 la presentation was imaged at a depth of 70 km, a value which is also dependent upon ,  63  the Moho thickness.  4.1.2 Non-uniqueness of the Model Since two different models may both produce an acceptable fit to the data, be wellresolved, and appear physically acceptable, non-uniqueness is inherent to the final model. The scope of the non-uniqueness issue in this problem is outlined below. The non-unique nature of the solution is evident from the subjective choice of node placement. Furthermore, for a certain choice of node placement, different a priori estimates of parameter uncertainties, different damping parameters, and different starting models are associated with a different final model. Even within a certain parameterization, there is a sense of non-uniqueness associated with any nodal value output from the RAYINVR inversion — each calculation has an associated standard error, indicating a range of feasible values. Thus, it must be stressed that the final model presented above is a model which fits the data, not a unique solution of the inverse problem.  4.2 Spatial Resolution and Absolute Parameter Uncertainty Zelt and Smith (1991) described tests which provide estimates of (i) the spatial resolution of a model about a specific velocity or boundary node, and (ii) the absolute parameter uncertainty of a particular node (note that the term "spatial resolution" is used here to distinguish this quantity from the numerical resolution discussed in a previous chapter). The first step in estimating spatial resolution is to choose a model parameter and perturb it by a measure on the order of its estimated uncertainty o — the perturbation should be large enough to change travel time, but small enough that ray path pattern is not disturbed significantly. The next step is to trace rays through the perturbed model to calculate the travel times corresponding to the receiver locations and assign the observed pick uncertainty at each location. Finally, the parameter is reset to its non-perturbed value, and the calculated  64  data are inverted to determine all the model parameters that were computed at the same time as the selected parameter during the inversion for the final model. The spatial resolution about the selected parameter is indicated by the amount that the values of adjacent parameters differ from the corresponding values in the non-perturbed model. If the model is well-resolved about the selected parameter, then all other parameter values will resemble the non-perturbed values; poor resolution about selected parameter will result in the perturbation being "smeared" into the adjacent parameters. To estimate the absolute uncertainty of a model parameter, Zelt and Smith (1991) suggested another test, the first step of which again involves choosing a model parameter to perturb. The second step is to invert all other parameters that were determined at the same time as the selected parameter during the inversion for the final model, while holding the selected parameter fixed at its perturbed position. Finally, the size of the perturbation is continually increased until the data cannot be fit as well as they were in the non-perturbed model according to 1) the ability to trace rays to all observations, and 2) the result of an F—test comparing the chi-squared values of the two final models. The statistical F—test involves a ratio calculated such that values either >>1 or <<1 indicate very significant differences between two samples; in this application, it reveals in a statistical sense whether or not the fit of a perturbed model has deteriorated significantly from that of the final model. The size of the largest perturbation which leads to an equivalent fit to the data gives an indication of the absolute uncertainty of the selected parameter. This test for absolute parameter uncertainty, however, does not account for the added constraints of forward amplitude modelling. As acknowledged by Zelt and Smith (1991), this type of analysis on every node would take longer than the inversion procedure; they suggest that testing a number of representative velocity and boundary nodes may be sufficient. The estimated spatial resolution and absolute uncertainty values thus obtained for the final model are shown in Table 7. Generally, the resolutions decrease and the uncertainties  65  increase with depth in the model. As well, this testing procedure revealed a trade-off between resolution and uncertainty based upon node spacing: in comparison with widelyspaced nodes, those which were more closely spaced were shown to have both a decrease in resolution and an improvement in absolute uncertainty.  ^ Nodes Velocities in near-surface layer (layer 1) Velocities at top of upper crust (layer 2) Depth nodes at base of upper crust (boundary 6) Velocities at top of middle crust (layer 6) Depth nodes at base of middle crust (boundary 7) Velocities at top of lower crust (layer 7) Depth nodes of R3 in lower crust (boundary 8) Velocities within the lower crust (top of layer 8) Depth nodes at crust-mantle boundary (boundary 9) Velocities in the upper mantle  Estimated Estimated ^ lateral absolute ^ resolution uncertainty (km) 20-70 10-15 35 40 35 50 50 30 50 *  0.5 km/s 0.1 km/s 1.5 km 0.1-0.2 km/s 1.5 km 0.2-0.3 km/s 1.5 km 0.2-0.3 km/s 1.5 km 0.1 km/s**  * There were not sufficient Pn data to determine this quantity. **based upon an assumed value of Moho thickness. Table 7. Estimated lateral resolution of the final velocity model about the given nodes, and absolute uncertainty of velocity and depth nodes.  66  V. DISCUSSION AND CONCLUSIONS 5.1 Comparison with LITHOPROBE Refraction Interpretations The refraction lines of SCoRE '89 — Lines 1, 2 and 3 in Figure 2 — were all within the proximity of Line 10. As such, the interpretations of these profiles which exist at present may be used for comparison with the results obtained from the Line 10 study. The data from Line 1 and Line 3 have been interpreted by Zelt et al. 1992a and 1992b, respectively. A preliminary interpretation of Line 2 has been made by McLean (in preparation). The Line 2 data modelled by McLean (in preparation) were provided by a profile which ran from the eastern Insular Belt across the Coast Belt and into the western Intermontane Belt (see Figure 2). The proposed model features upper crustal velocities of - 6.25 km/s ,  (to -10 km depth), mid-crustal velocities of 6.35 —6.50 km/s to 17 km depth overlying materials of 6.55 km/s average velocity which extend to 24 kin depth, and lower crustal ,  velocities of 6.8 km/s. The Moho has an average depth of approximately 34 km, and the underlying upper mantle has velocity of 7.9 km/s. These preliminary results of the line 2 interpretation are important to this study since Line 2 intersected Line 10 (Figure 2; reflection profile 88-13 also intersected Line 10 at this location). It is important to ascertain that the interpretations of the two refraction data sets agree at their intersection point, or to provide a plausible explanation for any discrepancy. The Line 2 data interpreted by McLean (in preparation) exhibit the same distinct features —i.e. Pg phase confined to upper crustal region with underlying crust giving rise to three prominent reflection phases — as the Line 10 data, and the parameterization chosen in the modelling procedure was very similar. This suggests a comparison as shown in Table 8, where the approximate values of the main crustal features of the models at their intersection point (shown in Figure 11a) are displayed. Clearly, the two models agree closely, with the  67  ^ Model Feature average velocity above R1 depth to R1 reflector average velocity between R1 & R2 depth to R2 reflector average velocity between R2 & R3 depth to R3 reflector average velocity between R3 & Moho depth to Moho  Line 2^Line 10 6.30 km/s 6.25 km/s 10.0 km 10.0 km 6.40 km/s 6.40 km/s 22.0 km 20.0 Ian 6.50 km/s 6.55 km/s 29.0 km 24.5 km 6.70 km/s 6.80 km/s 34.0 km 35.0 km  Table 8. Velocity values in the crust and depth to reflectors for Line 2 of SCoRE '89 and Line 10 of SCoRE '90 at the intersection point of Line 10 and Line 2 (see Figure 2). exception of the depth to the R3 reflecting horizon. Explanations for this discrepancy include the 3—D nature of Line 10 at this location being ignored in the 2—D approach, the combination of the absolute errors resulting from the modelling procedures, the possibility that a different cycle of the phase was picked (testing has shown that a difference of one cycle in time would be equivalent to an "s1 km shift in depth), and the possibility of a structural variation in the reflector (although the near-coincidence of the region of R3 constraint on each of the two models would necessitate a rather abrupt inflection). Alternatively, the most reasonable explanation for this 4.5 km difference may be that the source of the McLean (in preparation) R3 reflection is an entirely different horizon than the one defined by the Line 10 data. Line 3 of SCoRE '89 extended from the Insular Belt through the southernmost Coast Belt to the western border of the Intermontane Belt (Figure 2), passing about 20 km south of the southeastern end of Line 10. An interpretation of the in-line data has been presented by Zelt et al. (1992b). The proposed model shows average velocities in the upper crust (to a depth of  ,  11 km) of 6.3 km/s in the west and 6.2 km/s in the east; in the middle crust (to  a depth of 22 km) of 6.7 km/s in the west and 6.35 km/s in the east; in the lower crust ,  of 7.0 km/s in the west and 6.8km/s in the east; and in the upper mantle of 7.65 km/s. The  68  Moho has an average depth of 37 km throughout most of the model, except in the west, where a depth of 32 km is reached at the Insular Belt-CB border. An important feature ,  of the model is an abrupt east-west variation interpreted in the Harrison Fault region; Zelt et al. (1992b) have identified this as the location of the collision zone between the Insular superterrane (which exists as Wrangellia in this region) and the Intermontane superterrane. The southward projection of Line 10 intersects Line 3 in the region of Harrison Fault, i.e. at the Insular superterrane-Intermontane superterrane boundary. Table 9 shows the velocity values in the upper crust (to 10 km depth), middle crust (10-22 km depth), and lower crust (22 km — Moho depth) for Line 10 and Line 3 at the intersection location. The Moho depth was interpreted for Line 3 to be -34.5, as was the case for Line 10. "Line 3W" and "Line ,  3E" refer to the portions of the Line 3 model directly to the west and east of the Harrison Fault, respectively. Table 9 shows that the velocities determined for the southeastern portion of Line 10 are much closer to Zelt et al.'s (1992b) values for the eastern side of the collision zone, i.e. those of the Intermontane superterrane.  Model Feature velocity of upper crust (km/s) velocity of middle crust (km/s) velocity of lower crust (km/s)  Line 3W 6.30 6.70 7.00  Line 10 6.20 6.40 6.55  Line 3E 6.20 6.35 6.8  Table 9. Velocity values in the crust for Line 10 of SCoRE '90 and Line 3 of SCoRE '89 at the intersection point of Line 10's projection onto Line 3 (see Figure 11a). The upper crust 0-10 km depth, the middle crust 10-22 km depth, the lower crust ti 22-34.5 km depth. The "W" and "E" refer to the portions of Line 3 to the west and east, respectively, of the intersection point. See Figure 14a for a pictorial representation. The Intermontane Belt to the east of the CB was the site of Line 1 of the 1989 Southern Cordillera Refraction Experiment (SCoRE '89, Figure 2). Zelt et al. (1992a) have presented a 2-D model for the along-strike line which shows (1) upper crustal and mid-crustal velocities varying laterally from 6.0 to 6.4 km/s; (2) lower crustal velocities averaging 6.5km/s and 6.7  69  FIGURE 14: (a) and (b) show comparisons of 1—D average velocity (in km/s) versus depth profiles for the refraction interpretations indicated (see Tables 9 and 10). The shaded region in the Line I (VISP '80) presentation represents an area that is poorly constrained. The "W" and "E" represent the portions of Line 3 to the west and east, respectively, of the intersection point. The Moho depths are indicated with a thick horizontal line.  70  km/s at the top and the base; (3) a transitional Moho of 1-3 km thickness at an average depth of 33 km overlying an upper mantle with velocities of 7.9 km/s in the south and 7.7 km/s in the north; and (4) an upper mantle reflector at --16 km below the Moho. thought to represent ,  the base of the lithosphere, which is considered to be thin throughout the Intermontane Belt. For the calculation of upper mantle velocities, Zelt et al. (1992a) had many observations of Pn to use in the inversion algorithm, while the upper mantle velocity determined using the Line 10 data is not well-constrained, as will be discussed later. The velocities beneath Line 1 of the Intermontane Belt as interpreted by Zelt et al. (1992a) bear a strong resemblance to the Line 10 velocities, as shown in Table 10. Furthermore, the interpretation of Spence et al. (1985) of VISP 80 Line I has shown Wrangellian velocities in the upper 18 km which are much higher (at 6.3-6.7 km/s, see Table 10) than the velocities determined through the Line 10 analysis.  Crustal Velocity (km/s) Line 10 ^ 6.20 (to 10 km) upper crust ^ 6.40 (to 22 km) middle crust ^ 6.55 (to 34.5 km) lower crust  Line 1 (SCoRE '89) 6.20 (to 12.5 km) 6.20-6.40 (to 23 km) 6.60 (to 33 km)  Line I (VISP '80) 6.30 (to 9 km) 6.70 (to 18 km)  * The Line I model of Spence et al. (1985) was not well-constrained below 18 km depth. Table 10. Velocity values in the crust (to depths approximated in parentheses) for Line 10 of SCoRE '90, Line 1 of SCoRE '89, and the easternmost part of VISP '80 Line I. See Figure 14b for a pictorial representation. Clearly, this overview of both concurrent and previous regional refraction studies shows that the velocities beneath Line 10 are similar to those derived in studies of the Intermontane superterrane, and significantly slower than those interpreted for elements of the Insular superterrane (i.e. Wrangellia); this information proved to be critical to the tectonic interpretation of results from this study.  71  The Zelt et al. (1992a) interpretation of the Line 1 upper mantle reflector being the base of the lithosphere is also significant. The review of geophysical, geological, and petrological studies from the Cordillera by Gough (1986) indicated that the lithosphere is thin throughout the Intermontane Belt, and tapers to deeper levels beneath the Coast Belt. This suggests that the upper mantle reflector below Line 10 could, at a depth of " 35 km beneath the Moho, ,  represent a western counterpart to Zelt et al.'s reflection from the top of the asthenosphere, an idea supported by the fact that a negative velocity contrast across the R reflecting horizon was shown to provide as good a fit to the observed data as a positive velocity contrast. Another horizon that was considered as a possible source of the upper mantle reflection R was the top of the subducting oceanic crust. However, the depth of the subducting plate has been inferred to be  ,  100 km at a point - 50 km west of Line 10 based on hypocenters ,  (Rogers et al. 1990) and the presence of the Garibaldi volcanic chain (Dickinson 1970). This precluded the top of the subducted plate as a contingent source of the upper mantle reflection R, which was imaged at a depth of 70 km. Further, the multidisciplinary review ,  of Gough (1986) has suggested that the subducted plate loses its identity before it reaches the region beneath Line 10.  5.2 Comparison with Other Geophysical Studies — Reflection, Heat Flow, Gravity, and Geomagnetic 5.2.1 Comparison with LITHOPROBE Reflection Data and Interpretations In 1988, nearly 350 km of deep seismic reflection data were acquired in the southern CB through a LITHOPROBE seismic reflection transect (profiles 11 to 18; Figure 1 of Varsek et al. 1992). Lines 88-13, 88-14, 88-17, and, to a lesser degree, 88-18 were in the vicinity of the present study area (Figure 2). The seismic sections generally show east-dipping upper crustal (i.e. above P-3.5 s) reflectors, which are truncated by the subhorizontal to west-dipping structures that characterize the middle and lower crust (Varsek et al. 1992). Monger and  72  Journeay (1992) have provided a preliminary interpretation of the data from these profiles, and Varsek et al. (1992) have made a more detailed interpretation. Figure 15 shows the reflection data for (a) profile 88-13, and (b) profiles 88-14 and 88-17; the location of Line 10 is indicated. As shown in Figure 2, the westernmost limit of 88-13 abutted against Line 10 at the southeast side of SP-48, and 88-14 and 88-17 provided a crossing of Line 10, with the eastern end of 88-14 coinciding with the refraction line and the western limit of 88-17 falling 11 km to the east of it; the locations of the reflection ,  lines with respect to Line 10 are shown on Figure 1 la. The proximity of these profiles to Line 10 makes it important to consider the nature of the reflection data in both areas. The section of Figure 15a shows that the western part of profile 13 is characterized by east-dipping structures (indicated by "A") in the upper crust. The middle crust shows zones which dip slightly to the west ("B") and truncate the east-dipping upper crustal reflectors. Nearly-linear features are apparent in the lower crust, as are the east-dipping reflections ("C") which are visible at -6.0 and 9.5 s underneath the Line 10 location. The pronounced ,  package of reflections from the Moho (denoted with an "M") begins below this at -10.9 s and has a travel time thickness of 0.3 s; these values vary laterally along the profile, with ,  the travel time thickness reaching up to "0.8 s. The overall crustal fabric for profile 88-17 (Figure 15b) is again one of upper crustal east-dipping reflectors ("A") truncated by west-dipping mid-crustal reflections ("B") which overlie flat reflectors in the lower crust. Prominent crustal reflections occur at the Line 10 site on the 88-14 reflection section in zones of enhanced reflectivity between approximately 1.5 and 3.5 s and 6.0 and 7.0 s, and at roughly 10.0 s (each of these is indicated with an "*"); all seem continuous with reflections in the 88-17 data. Reflections from the Moho are not visible for 88-14, but for 88-17 they again appear as a package of prominent reflections, with a travel time of 10.8 s at the top beneath Line 10 (denoted with an "M"). This package ,  of Moho reflections begins at -10.8 to 11.0 s across the profile, and is of variable thickness.  Line 10  FIGURE 15:(a) LITHOPROBE reflection profile 88-13, migrated and coherency filtered (Varsek et al. 1992), with the Line 10 location indicated at the surface. "A" indicates upper crustal east-dipping reflectors, "B" mid-crustal west-dipping reflectors, "C" east-dipping reflections embedded in the generally flat reflections of the lower crust, and "M" the Moho reflections.  Line 10  BKFS ---•-•%--,-  ''^  : -''.*:-...2",': .  _.. . _  •-•■•4•'...  , '.. ,^■ -7--7,--'-..--':  -,  _  . .r-  -.  -----i-.--i-..  5...:. 4'.,‘, • n -,—,__ .-..7 .. :. ,:k.  „-  -."--, -.-.-..--7,"'7: 2;F:::-.  .  .  --..%--^-_',' •  •\  -  .^.  -2,----,--7 :  :--- '. '--, ...-.,---  `....,,..4.,--F-2.t.  ''. ,,---,"---,•*- -.-'--  - 7^.1:.-^-, -'  2  .^.  . .o..  ,  ..m.  6  ..,--.' •,.._,---. , _^- ...,  "--  4  --_,-.  8 ''''..,----  -^.-r-_ -k'- •,,^- ,-  M  ---._--....--,... --,--•:-"-  ..-,---.,--  ,  -..,..--"..7,..,--7" .. -------- •: • • - -^_ .,-----,... ,_^.., .  ,  .  ..--,..' -.%-•  -^:•••-•-• "._, ..'".:--"----";:,-ir-F:.;_"- -.-.  10 12 14  •  50 km  (b)  FIGURE 15:(b) LITHOPROBE reflection profiles 88-14 and 88-17, migrated and coherency filtered (Varsek et al. 1992), with the Line 10 location indicated at the surface. BKFS indicates the Bralome-Kwoiek Creek fault system, and "*" zones of enhanced reflectivity discussed in the text; "A", "B", and "M" are as defined in (a).  75  Many of the upper crustal structures visible in the reflection sections undoubtedly arose from the Mesozoic crustal shortening adjustments that accommodated terrane accretion. The similarity of the data interpretations (i.e. Monger and Journeay 1992, Varsek et al. 1992) for the upper crust have shown that observed seismic reflectors may indeed be successfully correlated with the surface expressions of terrane boundaries and associated terrane properties. However, the extensions of the upper crustal reflection interpretations to lower crustal depths have not been certain; this is a point to which this discussion returns. Figure 16a shows (i) the 1—D velocity profile of the Line 10 interpretation at the location of profile 88-13 plotted versus two-way travel time, and (ii) a simplified schematic of the 88-13 reflection interpretation of Monger and Journeay (1992) at the location of Line 10. Figure 16b is the counterpart of Figure 16a, the difference being that part (ii) represents the interpretation of Varsek et al. (1992). Figure 17a shows (i) a simplified schematic of the 88-14 reflection interpretation of Monger and Journeay (1992) at the location of Line 10, (ii) the 1—D velocity profile plotted versus two-way travel time of the Line 10 interpretation at the location of the 88-14/88-17 crossing of Line 10, and (iii) a simplified schematic of the 88-17 reflection interpretation of Monger and Journeay (1992) at the location of Line 10. Figure 17b is the counterpart of Figure 17a, with parts (i) and (ii) representing the work of Varsek et al. (1992). Monger and Journeay (1992) did not define a terrane affiliation for all units; these have been added to Figures 16a and 17a where it was deemed appropriate. The conversion of the final model from depth to two-way travel time for Line 10 velocities was accomplished by the integration of reciprocal velocities from the surface downward. Since the refraction model changes very slightly between its intersection with 88-14 and its intersection with the projection of 88-17 (see Figure 11a), and structural continuity between the southwest end of profile 17 and profile 14 is clear from the reflection data, Figures 17a and 17b show a 2—D representation of the crossing. The first important implication of the Figure 16 and Figure 17 presentations is that  FIGURE 16: Figure 16a shows (i) the 1—D velocity profile plotted versus two-way travel time of the Line 10 interpretation at the location of profile 88-13, and (ii) a simplified schematic of the 88-13 reflection interpretation of Monger and Joumeay (1992) at the location of Line 10. Figure 16b is the same as Figure 16a, except that part (ii) represents the interpretation of Varsek et al. (1992). "Wr" indicates that the Insular superterrane exists as Wrangellia in this region.  77 Velocity (km/s) 3 4 5 6 7 8 9  Line 10  1^I^I^IJ 1-Gambier (Overlap Assemblage  Cadwallader Terrane  upper  Gambier Overlap Assemblage Insular Superterrane  idc le  u16  E  -  F  -----  --  R2  lntermontane Superterrane  Intermontane Superterrane  F13  lowe r a ist  10  12  Upper Mantle  14  (iii)  )  (a)  Line 10  88-14  Line 10  88-17 from Monger and Joumeay (1992)  Velocity (knits) 3 4 5 6789  Gambier Overlap Assemblage  Cadwallader Terrane Gambier Overlap Assemblage  Shuksan Terrane  6  E  Insular Superterrane —————— Shuksan Terrane  8  10  12 Upper Mantle  14  (b) 88-14  (iii) Line 10  88-17 from Varsek et al. (1992)  78  the conversions of the final model from depth to two-way travel time allow the refractionimaged reflections to be compared with the reflection data shown in Figure 15. Although the reflection lines intersected Line 10 in the cases of 88-13 and 88-14, a precise correspondence is not expected due to the following points: (i) the possibility of crustal anisotropy; (ii) the possibility of lateral averaging in the refraction velocity model; (iii) a 3-D effect is introduced to the 88-13 case due to its intersection point being located on the only portion of Line 10 which veers significantly from the straight line between its endpoints; and (iv) the quality of the reflection data suffers at the extreme ends due to the loss of fold. At the intersection of Line 10 and 88-13, R1 is not constrained (see Figure 11a); this is consistent with the fact that there is no definite individual reflection imaged at this level ( 4.0 s) in the profile 13 data (Figure 15a). The R2 and R3 reflections were modelled at ,  -6.7 and 9.7 s, respectively (see Figure 16a (i) or Figure 16b (i)), and correspond well with structures visible in the reflection data (Figure 15a) and described above. The transitional refraction Moho was modelled at -10.8 s which roughly coincides with the reflection Moho; further discussion of this comparison appears later in this section. At the location of the 88-14/88-17 crossing of Line 10, R1 is again not constrained (see Figure 1 la); however, an event in the reflection data occurs at N3.5 s (Figure 15b), the time of Ri. Varsek et al. (1992) have interpreted this reflecting horizon as the subsurface continuation of the northern extension of the Thomas Lake fault (Figure 4); the rocks juxtaposed at this boundary are of theological properties that could not be resolved with the refraction experiment. The R2 and R3 reflectors were imaged at r-7.2 and 9.8 s, respectively, by the Line 10 data (see Figure 17a (ii) or Figure 17b (ii)). The reflection data (Figure 15b) show that the base of a high-reflectivity zone appears at '7.0 s on both the 88-14 and 88-17 sides of the crossing, and that a strong event occurs similarly at 9.8 s. The second significance of Figures 16 and 17 is that the differences of the Monger and Journeay (1992) and the Varsek et al. (1992) models show the tenuous nature of the  79  interpretations of deeper (i.e. below .-3.5 s) reflections, so that while the upper crustal interpretations are alike, below this the solutions diverge. For example, the source of the reflection corresponding R2 at the westernmost end of 88-13 is identified in Varsek et al. (1992) (Figure 16b (ii)) as the top of the "Insular ramp", a proposed unit containing a suture separating the Insular and the Intermontane superterranes; a very different wedge-like structure appears in the Monger and Journeay (1992) version (Figure 16a (ii)). Also, while the Monger and Journeay (1992) interpretation does not put forward a structural origin for the R3 reflector, the Varsek et al. (1992) model shows R3 to be generated at the bottom of the Insular ramp. Further examples appear in the comparison of Monger and Journeay's (1992) interpretation of the Line 10 crossing of 88-14 and 88-17 (Figure 17a (i) and (iii)) with that of Varsek et al. (1992) (Figure 17b (i) and (iii)). This comparison shows that the former has the R2 reflecting horizon embedded in an Intermontane unit, while the latter has it within the Insular superterrane. And while Varsek et al. (1992) did not interpret a structural source for the R3 reflection, Monger and Journeay (1992) have shown this as the top of a thrust ramp within the Intermontane superterrane. The crustal structures shown in Figures 16 and 17 do show an important similarity: the juxtaposition of Insular and InterMontane materials. As well, the presence of rocks classified with the Shuksan terrane, which represents the vestiges of an inter-terrane ocean basin, appears in the Varsek et al. (1992) interpretation. This shows that both of the reflection studies have concluded that the Line 10 region was the site of the collision zone between the Insular and Intermontane superterranes, as was the case with the Zelt et al. (1992b) refraction study of SCoRE '89 Line 3. It is not surprising that the interpretation of seismic reflection data across such a structurally complex features has lead to ambiguities in linking deep reflections with known structural features; in fact, the anticipation of such ambiguities provided part of the motivation  80  for a study designed in the manner of SCoRE '90. In the region of this study, such ambiguities have led to fundamentally different views of the nature of the superterrane collision which produced the CB. Varsek et al. (1992) have constructed a collisional model based upon crustal imbrication, as they identified the source of reflectivity as the lithological layering resulting from the interfingering of terranes. The Monger and Journeay (1992) model is one of crustal delamination, whereby the sub-horizontal reflections characterizing the middle and lower crust were interpreted as deep counterparts to east-vergent faults identified in the shallower crust to the east (i.e. the Bralorne-Kwoiek Creek fault system, Figure 15b), with Wrangellia being dispaced along this fault system over the top of the terrane material below. Figure 18a shows the interpretation of Monger and Journeay (1992) and Figure 18b that of Varsek et al. (1992) over the entire width of the reflection profiles. A fundamental difference in these two viewpoints is the extent to which Wrangellia, whose easternmost surface expressions are shown as x 's in Figure 4, has penetrated the Cordillera. The Monger and Journeay (1992) proposal has Wrangellia being pushed upward over the Intermontane superterrane, so that all that remains of it below Line 10 is a wedge overlying the rest of the middle and upper crust (Figures 16a and 17a; Figure 18a). The imbrication model of Varsek et al. (1992) results in Wrangellia existing much farther to the east in the middle and lower crust, so that thick fingers of it are proposed beneath the Line 10 location (Figures 16b and 17b; Figure 18b). The comparison of refraction results in the previous section has shown that if Wrangellia exists under Line 10, it must be limited to zones whose extent could not be resolved; thus, the crustal delamination model of Monger and Journeay (1992) is the one which is consistent with results derived in the present study. Third and finally, the summary provided in Figures 16 and 17 allows the correspondence of the refraction-imaged transitional Moho (layer 9 in Figure 12) with the crust-mantle boundaries in the reflection data to be verified. The Moho has been defined by Steinhart (1967) as "that level in the Earth where compressional wave velocity increased rapidly or  •^ -^ •^  81  Line 10 '^-Mc ...... .-  :.  ./. '.-• '.  -"41..C .1"4  ..  * .0.1.,y_  •4 4-"•-•  ....'  '----' :'-114.- .^-.•^•^•  •  •  :45.....,.:  •  • ..: '^-.... •  -6  ii, •  ,•••' v. ; - c:C - 8  fIN7=- Insular k -_ -, -7.7,:-C-=- . . -  -  -- - -  -  - -  7- •  •.:""7-V7--i-•••^—  •  —  ^-  10 12  (a) -^  e  14  from Monger and Journeay (1992) —^•  2  —^ •  . •^  —'^ •-••^— s‘s•-••••-• •  E  •  I•cmt--.^  -  4cC  8  12 14 50 km  from Varsek et al. (1992)  FIGURE 18: Interpretations of the reflection data from the 88-14/88-17 crossing of Line 10 by (a) Monger and Joumeay (1992) and (b) Varsek et al. (1992) overlie the reflection data (see Figure 15). Black lines indicate reflections which have been interpreted in the studies, the dotted area represents Insular superterrane, and the line-shaded area represents Intennontane superterrane (see Figures 16 and 17 for more detail in the region of Line 10). Where appropriate, some identifications missing from the original studies have been added.  82  discontinuously to a value between 7.6 and 8.6 km/s"; thus, the Moho must by definition be determined through refraction studies, which yield velocity information. For reflection data, the Moho has been identified by Klemperer et al. (1986) as the "deepest, high-amplitude, laterally extensive reflection or group of reflections". The question of how the reflection Moho relates to the conventional definition remains one of current debate (Jarchow and Thompson, 1989). The duration of the zones of Moho reflectivity on the 88-13 (Figure 15a) and 88-17 (Figure 15b) sections are consistent with the interpretation of the crust-mantle boundary as an alternating series. of high and low velocity material. As shown in the seismic section of 88-13 (Figure 15a), the pronounced zone of Moho reflectivity at the site of Line 10 is between ,  10.9 and 11.2s. The top of this zone coincides well with the top of the transitional refraction  Moho, which was modelled at -40.8 s. The travel time thickness of s observed for these well-defined reflections is equivalent to -4.9 km, which may be an indication of the thickness of the transitional Moho. The package of Moho reflections visible on the profile from 88-17 (Figure 15b) also begins at a travel time (-41.0 s) which corresponds to the top of the refraction-imaged Moho layer (-41.0 s); the thickness of the reflectivity package in this case is more difficult to discern. Monger and Journeay (1992) and Varsek et al. (1992) placed the Moho at the position of the deepest prominent reflection; this interpretation coincides with the base of the refraction-imaged transitional Moho, if its thickness is assumed to be about 2.0 km as inferred from the well-defined 88-13 Moho reflections (Figure 15a). The correlation of the top of the refraction-imaged transitional Moho with the top of the zone of prominent Moho reflectivity observed in the 88-13 and 88-17 sections is suggested by the data; similarly, the base of this transitional layer may correspond to the deepest extent of this Moho reflectivity, although it is important to note that these observations are subject to the refraction modelling error (Table 7).  83  5.2.2 Heat Flow Data and Interpretations Heat flow in the broad region of subduction zones characteristically follows a pattern: a band of low heat flow reaches from the trench to the volcanic arc about 200 km inland, and a much higher than normal heat flow exists for a great distance inland (Keen and Hyndman 1979). The heat flow measured across southwestern British Columbia in an area encompassing the LITHOPROBE Southern Cordillera Transect site exhibits this characteristic pattern (Lewis et al. 1988, 1992). The data show high heat flow values for most of the CB (-75-95 mW/m 2 ), and values almost as high for the Intermontane Belt; in comparison, heat flow values are low (<50 mW/m 2 ) throughout the Insular Belt. In Phanerozoic regions, including the southern Canadian Cordillera, the lower crust typically exhibits bands of subhorizontal reflectors. Klemperer (1987), Weyer et al. (1987) and others have indicated that a correlation exists between the depths to the tops of these characteristic bands and regional heat flow. Heat flow studies such as these have commonly emphasized the depths to the 450 and 730 °C isotherms, which are identified respectively as the approximate temperatures of the brittle-ductile transition zone in the crust, and the transition from a mineralogy in equilibrium with coexisting free water to a dry mineralogy. In a number of Cordilleran reflection sections, the 450 °C isotherm has corresponded to the top of the characteristic, sub-horizontal reflective bands; the 730 °C isotherm has corresponded to the bottom of these characteristic bands (Lewis et al. 1992). The origin of the reflections has been the subject of much discussion (e.g. Warner 1990) and support has been provided for the model of layered porosity, whereby the region between the transition zones described above is identified as the portion of the crust where thin pores can stay open, and fluid trapped within the pores generates high reflectivity coefficients if the porosity is layered. In the LITHOPROBE reflection data from the Line 10 region, however, the top and bottom of such reflective bands are not distinguishable, as a strong reflectivity is present throughout the entire crust (Figures 15a and 15b).  84  Lewis et al. (1992) have shown that the level of heat flow in the Line 10 region was -94 mW/m 2 ; corresponding crustal temperatures suggested a depth to the 450 °C isotherm of , 10.8 km and to the 730 °C of , 20.0 km; these depths coincide well with the depths of the R1 and R2 reflectors imaged in the refraction data (see Figure 11a). In the case of the R1 reflecting horizon, although its depth (-40.0 km) correlates well with that of the 450 °C isotherm, its sporadic nature (Figure 11a) does not closely resemble continuous reflection described by Lewis et al. (1992) as coming from the top of a zone of layered porosity in the crust. Also, the fact that the regions where R1 is constrained are equivalent to the portions of Line 10 which crossed highly-faulted regions (Wheeler and McFeely 1991) suggests faulting as a more likely source of the reflectivity. In terms of the R2 horizon, its depth (-22.0 km), continuous, nearly linear nature, and the relatively small velocity contrast across it (0.05-0.15 km/s) support the idea of its correspondence with the transition zone linked with the 730 °C isotherm. However, the reflective nature of the crust below this depth, which is visible on both the refraction data (i.e. R3 in Figures 5-10) and the reflection data (Figures 15a and 15b), means that the abrupt decrease in reflectivity associated with this transition zone cannot be uniquely identified, so it is not appropriate to attempt a definitive correlation with the R2 reflecting horizon. Nevertheless, the corresponding depths of the R2 horizon and the 730 °C isotherm do suggest that the interpreted structural boundaries „should be regarded with caution.  5.2.3 Gravity Data and Interpretations The Bouguer gravity anomaly map for the study area is shown in Figure 19 (from Riddihough and Seeman 1982); the map shows that this is a region of low gravity values. The gravity field across nearly all continental margins where subduction is taking place also shows a characteristic pattern (Keen and Hyndman 1979). This is typified by a band of very low gravity over the subduction trench, and a parallel band of high gravity extending inland for approximately 100 km. Much farther inland, there is usually a large region of low Bouguer gravity; this corresponds to the region of the present study. Riddihough (1979)  85  FIGURE 19: The Bouguer gravity anomaly map for the study area (from Riddihough and Seeman 1982), with contours of 10 mgal overlying the topography. The location of Line 10 is indicated, and the high gravity anomaly at the southeastern part of the profile is denoted with an "H".  86  pointed out that this low gravity could be attributed to low mantle density, possibly caused by high temperatures related to the high heat flow in the area. The gravity data coincident with Line 10 show that the profile is one of high negative anomalies, with lower negative values on the ends (-124 mgal in the northwest, and —87 mgal in the southeast) than in the middle, where values reach -200 mgal. Using the Barton (1986) velocity-density relationship, the final velocity model was converted to a density model for the purpose of gravity modelling. [Barton (1986) gave the mean, and the maximum and minimum bounds, of the measurements used for the conventional Nafe and Drake (1963) relationship; as well, Barton suggests slightly higher densities for seismic velocities greater than 6.5 km/s.] The algorithm available for gravity modelling was the interactive Fortran program SAKI (Webring 1985), which calculates and compares theoretical gravity responses of a structural model with profiles of observed data. The input model was constructed as a set of polygonal prisms of finite lateral extent, with the resulting format being 2.5—D. Gravity modelling in 2.5—D is suitable for the extrapolation of values to either side of a line; it is not suitable for modelling off-line changes. The application of the 2.5—D SAKI modelling algorithm to the data set from Line 10 was therefore limited, as the profile ran along the strike of a major geological feature, with quite different structures flanking it. Testing indicated that the velocity model could be converted into a density model which generally provided an accurate fit to the gravity data. The exception occurred in the modelling of the southeastern end of the profile, where out-of-plane effects impossible to account for with the 2.5—D modelling program rendered the algorithm inappropriate (Figure 19). Although the gravity modelling was not applicable, the increase in gravity at the southeastern portion of the line (Figure 19) must at least have a qualitative interpretation. Since the crustal thickness along Line 10 is relatively constant (Figure 11a), a possible explanation for the observed increase in values is that the densities of the rocks in this area  87  increase. The geological map of Wheeler and McFeely (1991) shows that the nature of the plutonic rocks which dominate the surroundings of Line 10 (Figure 4) changes in this southeastern region: the surface rocks are predominantly Mid-Cretaceous quartz diorites, as opposed to the Early Tertiary granodiorites/quartz diorites and Jura-Cretaceous quartz diorite packages which predominate at the northwest and central portions of the line, respectively. The gravity study of Dehler (1991) included the measured densities of rock samples from the CB. The wide range of density values observed per rock unit (e.g. 2600-2900 kg/m 3 for quartz diorite) lends plausibility to this possible interpretation.  5.2.4 Geomagnetic Data and Interpretations Natural geomagnetic induction studies (e.g. geomagnetic deep soundings, Caner 1970; magnetotellurics, Gough and Majorowicz 1991, Jones et al. 1992) provide information about the electrical conductivity of the crust and upper mantle. The early work of Caner (1970) showed that the entire region of the Cordillera in southern British Columbia to be one of high conductivity, or low resistivity. This was explained in the magnetotelluric (MT) investigation of Gough and Majorowicz (1991) as being caused by the wetness of Cordilleran rocks, which is linked in part with the ongoing subduction in the region. The Gough and Majorowicz (1991) work comprises one of two recent MT modelling studies in the general region of Line 10, the other being that of Jones et al. (1992). Both studies have interpreted subsurface structure in terms of electrical resistivity versus depth, and have shown an upper crust characterized by zones of high resistivity (i.e. > 300 S2•m) overlying a conducting middle and lower crust; the resistivity variations have been explained as resulting from variations in fracture densities. The resistive zones in the upper crust have been identified as granitoid plutons, with low fracture densities. The depth extent of the batholiths interpreted in the two MT studies ranges between 10 and 22 km. The vertical homogeneity of the Line 10 velocity model at this depth level suggests that the extent of the CB batholiths has not been determined using the refraction  88  data. This does not eliminate the possibility of the change in rock properties at the limit of a batholith generating reflectivity, since reflections identified as  the R1 and the R2 phase  come from depths ("40.0 and 22.0 km, respectively) that could correspond with the base of MT-imaged batholiths. A good example of this is provided by the intersection point of one of the Gough and Majorowicz (1991) MT profiles with Line 10, which occurs at approximately the site of SP-48 (the "E" site on their profile). Gough and Majorowicz image a resistive structure of which they say there is "little doubt that the ... mass ... is the main bulk of the Coast Plutonic granodiorites"; the structure reaches a depth of , 22.0 km, which corresponds with the depth of the R2 reflecting horizon at this location ( , 21.0 km). Such a scenario is less likely for the R1 reflections, however, since the northwestern region where most of the constraint on the R1 phase was located (Figure 11a) corresponds to the region of Line 10 where the surface rocks were terranes, as opposed to large granitic plutons (Figure 4).  5.3 Conclusions Through the application of successive iterations of inverse travel time modelling and forward amplitude modelling to the refraction data from Line 10 (Figures 2 and 4), a velocity structure was developed for the central Coast Belt (Figure 11a). The model features a thin ( , 2 km) near-surface layer with a velocity of 3.9 km/s at the top and a high (0.5 s -1 ) velocity gradient. Zones of high and low velocity overprint this average velocity structure, with the most pronounced high velocity zone possibly corresponding to a geological feature (the Cadwallader terrane). The surface layer overlies three major crustal strata: (i) the upper crust, which has an average depth extent of 10 km, velocity of 6.2 km/s, and velocity gradient of 0.02 s -1 ; (ii) the middle crust, which extends to 22 km depth and has a velocity of 6.35 km/s; and (iii) the lower crust, which consists of an 8 km thick upper unit with a velocity of 6.5 km/s overlying a unit of 4 km thickness and 6.7 km/s velocity. The crustal package dips slightly downward to the southeast, with the average depth to the Moho being 34.5 km. In order that synthetic seismograms provided a close match to the observed data, the Moho  89  was modelled as a thin transition layer, with velocities varying laterally from -7.40 to 7.80 km/s. The upper mantle velocity ranged from 8.05 to 8.15 km/s in accordance with a Moho thicknesses of between 0 to 2 km (0 km thickness represents a first-order discontinuity.) This upper mantle velocity estimate must be regarded with caution, however, as the observed data used for the inversion were sparse (only 33 observations), and the refraction data do not constrain the thickness of the transitional Moho. Furthermore, this relatively high velocity estimate is not consistent with the high heat flow, low gravity, and low electrical resistivity observed in the area. An upper mantle reflection was imaged at a depth ranging from -66 to 70 km, in keeping with a Moho thickness of 0 to 2 km. The reflection Moho observed on each of the 88-13 and 88-17 sections appears as a band of enhanced reflectivity with two-way travel time thickness varying from r-0.2 to 0.8s (or equivalently 0.7 to 3.0 km), with an upper extent corresponding with the top ,  of the refraction-imaged Moho. This observation is consistent with the interpretation of the refraction Moho as a band of alternating high and low velocity materials (Mooney and Brocher 1987), and suggests that its thickness, which is not constrained by the refraction data, could be constrained by the thickness of the bands. The Moho thickness of 2 km shown in the final model (Figure 11a) is consistent with the thickness of the Moho at the Line 10/88-13 intersection point, where the thickness of the pronounced band of reflectivity is clearly-defined. In general, the lateral velocity variations determined along the refraction line are minimal — this may be related to the approximately along-strike orientation of the profile (Figure 2, inset). Despite this lateral velocity homogeneity, the limited study of regional gravity has indicated that the density of the materials below Line 10 may increase in the southeast. This may be consistent with a change in surface geology in this region (Wheeler and McFeely 1991), in which the age of the plutonic rocks which dominate the region of Line 10 becomes Mid-Cretaceous, as opposed to the predominantly Jura-Cretaceous age of plutons to the north.  90  The strongest wide-angle reflections present in the SCoRE '90 Line 10 data are generated from the Moho; the origin of the three weaker crustal reflections, however, cannot be ascertained through the refraction data alone. The comparisons with other geophysical studies made in the preceding section have indicated that three possible types of reflecting horizons may exist at the depth of those imaged in the crust: (i) structural features (e.g. fault zones), (ii) the transition zones at the top and bottom of a region of layered porosity within the crust, and (iii) lithological contrasts (e.g. the boundary of a batholith and the surrounding terrane material). It is probable that the source of the R1 reflections was fault zones; the source of the R2 reflections could not be uniquely determined, and the only possible origin of the R3 reflections uncovered in this study are deep thrust faults which have been linked to the collision of the Insular and Intermontane superterranes (Monger and Journeay 1992). The upper mantle reflection R may represent the base of the lithosphere. Interpretations of seismic profiles which crossed, or nearly crossed, Line 10 have been provided by the SCoRE'89 Line 3 refraction study of Zelt et al. (1992b), and the Monger and Journeay (1992) and Varsek et al. (1992) studies of LITHOPROBE reflection profiles 88-13, 88-14, and 88-17 (Figure 2). These studies have determined that the location of Line 10 coincides with the collision zone between the Insular and the Intermontane terranes, a geological suture whose surface expression has been obscured by the granitic intrusions of the Coast Belt. Recent refraction interpretations for the southern Cordillera regions surrounding Line 10 include the Spence et al. (1985) study of VISP '80 Line I of the Insular Belt, and the Zelt et al. (1992a) modelling of SCoRE '89 data from the Intermontane Belt, as well as the study of Zelt et al. (1992b) (Figure 2). Spence et al. (1985) and Zelt et al. (1992b) have determined velocities for the Insular superterrane — which exists as Wrangellia in this region — that are significantly higher than those interpreted for the Line 10 profile (Tables 9 and 10, respectively); in contrast, the velocities interpreted by Zelt et al. (1992a) and Zelt et  91  al. (1992b) for the Intermontane superterrane are very similar to the Line 10 values (Tables 9 and 10). Although the presence of Wrangellia is to be expected in the proposed InsularIntermontane collision zone under Line 10, the refraction velocity model (Figure 11a) shows that it must be limited to zones whose extent cannot be resolved using the refraction data. Thus, the results of this study suggest that Line 10 represents approximately the easternmost extent reached by Wrangellia in its complex emplacement against the ancient margin of North America, an idea which is expressed in the interpretation shown in Figure 1 lb. The collisional model of Monger and Journeay (1992) involves the upward displacement of Wrangellia over the Intermontane rocks of the ancient North American margin, while that of Varsek et al. (1992) consists of the imbrication of Insular and Intermontane material (Figures 18a and 18b). Since the extension of Wrangellia to the east of Line 10 in the Varsek et al. (1992) model entails the local presence of thick Wrangellian layers, the Monger and Journeay (1992) interpretation which implies its presence as a wedge confined to the upper middle crust is favored by the results of this study.  92  BIBLIOGRAPHY  AM, K. and RICHARDS, P.G. 1980. Quantitative Seismology, vol. 2, W.H. Freeman, New York. 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, 87: 195-208. BENZ, H.M., SMITH, R.B., and MOONEY, W.D. 1990. Crustal structure of the north western Basin and Range Province from the 1986 Program for Array Seismic Studies of the Continental Lithospheric seismic experiment. Journal of Geophysical Research, 95: 21823-21842. BERRY, M.J. and FORSYTH, D.A. 1975. Structure of the Canadian Cordillera from seismic refraction and other data. Canadian Journal of Earth Sciences, 12: 182-208. BOLT, B.A., and GUTDEUTSCH, R. 1982. Reinterpretation by ray tracing of a transverse seismic profile through the California Sierra Nevada, Part I. Bulletin of the Seismological Society of America, 72: 889-900. BRAILE, L.W. 1973. Inversion of crustal seismic refraction and reflection data. Journal of Geophysical Research, 78: 7738-7744. BRAILE, L.W. and CHIANG, C.S. 1986. The continental Mohorovi6i6 discontinuity: Results from near-vertical and wide-angle seismic reflection studies. In Reflection Seismology: A Global Perspective, Geodynamics Series, edited by M. Barazangi and L. Brown, Washington, D.C., vol.13, pp. 257-272. BURIANYK, M.J.A., KANASEWICH, E.R., ELLIS, R.M., and CLOWES, R.M. 1992. SCoRE '90: Southern Cordillera Refraction Experiment Field Report and Data Set Description. LITHOPROBE Report No. 21, LITHOPROBE Secretariat, University of British Columbia, Vancouver, 30 pp. CANER, B., 1970. Electrical conductivity in western Canada and petrological implications. Journal of Geomagnetism and Geoelectricity, 22: 113-129 tERVENt V., MOLOTKOV, I., and PkNefK, I. 1977. Ray method in seismology. Charles University Press, Prague. CHIU, S.K.L., and STEWART, R.R. 1987. Tomographic determination of three-dimensional seismic velocity structure using well logs, vertical seismic profiles, and surface seismic data. Geophysics, 52, 1085-1098. CLOWES, R.M. (editor) 1989. LITHOPROBE Phase III Proposal — The Evolution of a Continent, LITHOPROBE Secretariat, The University of British Columbia, Vancouver. 213 pp. CLOWES, R.M., COOK, F.A., GREEN, A.G., KEEN, C.E., LUDDEN, J.N., PERCIVAL, J.A., QUINLAN, G.A., and WEST, G.F. 1992. LITHOPROBE — New perspectives on crustal evolution. Canadian Journal of Earth Sciences. (In press.)  93  COLES, R.L., and CURRIE, R.G. 1977. Magnetic anomalies and rock magnetizations in the southern Coast Mountains, British Columbia: possible relation to subduction. Canadian Journal of Earth Sciences, 14: 1753-1770. CONEY, P.J. JONES, D.L., and MONGER, J.W.H. 1980. Cordilleran suspect terranes. Nature, 288: 329-333. COOK, F.A., GREEN. A.G., SIMONY, P.S., PRICE, R.A., PARRISH, R.M., MILKEREIT, B., GORDY, P.L., BROWN, R.L., COFLIN, K.C., and PATENAUDE, C., 1988. LITHOPROBE seismic reflection structure of the southeastern Canadian Cordillera: initial results. Tectonics, 7:157-180. COOK, F.A., VARSEK, J.L., CLOWES, R.M., KANASEWICH, E.R., SPENCER, S.C., PARRISH, R.R., BROWN, R.L., CARR, S.D., JOHNSON, B.J. and PRICE, R.A. 1991. LITHOPROBE crustal reflection cross section of the Southern Canadian Cordillera: Foreland thrust and fold belt to Fraser River fault. Tectonics. (In press.) COOK, F.A., VARSEK, J.L., CLOWES, R.M., KANASEWICH, E.R., SPENCER, C.S., PARRISH, R.R., BROWN, R.L., CARR, S.D., JOHNSON, B.J., and PRICE, R.A. 1992. LITHOPROBE crustal reflection cross section of the southern Canadian Cordillera I: foreland thrust and fold belt to Fraser River fault. Tectonics, 11: 12-35. CUMMING, W.B., CLOWES, R.M., and ELLIS, R.M. 1979. Crustal structure from a seismic refraction profile across southern British Columbia. Canadian Journal of Earth Sciences, 16: 1024-1040. DICKINSON, W.R. 1970. Relations of andesites, granites, and derivative sandstones to arc-trench tectonics. Review of Geophysics and Space Physics, 8: 813-860. DAVIS, G.A., J.W.H. MONGER, and BURCHFIEL, B.C. 1978. Mesozoic construction of the Cordilleran "collage", central British Columbia to California. In Mesozoic paleography of the Western United states, Pacific Coast Paleography Symposium 2, edited by D.G. Howell and K.A. McDougall, Society of Economic Paleontologists and Mineralogists, Pacific Section, Los Angeles, California, pp.1-32. DEHLER, S.A. 1991. Intergrated geophysical modelling of the northern Cascadia subduction zone. Ph.D. thesis, University of British Columbia, 151 pp. DREW, J.J., and CLOWES, R.M. 1990. A re-interpretation of the seismic structure across the active subduction zone of western Canada. In Studies of laterally heterogeneous structures using seismic refraction and reflection data, edited by A.G. Green, Geological Survey of Canada, Paper 89-13, pp.115-132. FRIEDMAN, R.M. and ARMSTRONG, R.L. 1992. Jurassic and Cretaceous geochronology of the southern Coast Belt, B.C., 49°-51°N. Geological Society of America Memoir. (Submitted.) GABRIELSE, H., and YORATH, C.J. 1989. DNAG #4. The Cordilleran Orogen in Canada. Geoscience Canada, 16: 67-83.  94  GEHRELS, G.E., and SALEEBY, J.B. 1985. Constraints and 'speculations on the displacement and accretionary history of the Alexander-Wrangellia-Peninsular superterrane {abstract}. Geological Society of America, Abstracts with Programs, 17: 356. GOUGH, D.I. 1986. Mantle upflow mechanics in the Canadian Cordillera. Journal of Geophysical Research, 91: 1909-1919. GOUGH, D.I. and MAJOROWICZ 1992. Magnetotelluric soundings, structure, and fluids in the southern Canadian Cordillera. Canadian Journal of Earth Sciences, 29: 609-620. HALE, L.D. and THOMPSON, G.A. 1982. The seismic reflection character of the continental Mohorovi6id discontinuity. Journal of Geophysical Research, 87: 4625-4635. JARCHOW, C.M., and THOMPSON, G.A. 1989. The nature of the Mohorovitid discontinuity. Annual Review of Earth and Planetary Sciences, 17: 475-506. JOHNSON, S.H., COUCH, R.W., GEMPERLE, M., and BANKS, E.R. 1972. Seismic refraction measurements in southeast Alaska and western British Columbia. Canadian Journal of Earth Sciences, 9: 1756-1765. JONES, A.G., GOUGH, D.I., KURTZ, R.D., DeLAURIER, J.M., BOERNER, D.E., CRAVEN, J.A., ELLIS, R.G. and McNIECE, G.W. 1992. Electromagnetic images of regional structure in the southern Canadian Cordillera. Geophysical Research Letters. (In press.) JOURNEAY, J.M. 1990. Structural and tectonic framework of the southern Coast Belt, British Columbia. In Current research, part E. Geological Survey of Canada, Paper 90-1E, pp. 183-197. JOURNEAY, J.M., and FRIEDMAN, R.M. 1992. Structural setting and geochronometry of the Lillooet River Fault system, southwestern British Columbia. Tectonics. (In press.) HOUGH, S.E., and ANDERSON, J.G. 1988. High-frequency spectra observed at Anza, California: implications for Q structure. Bulletin of the Seismological Society of America, 78: 692-707. KANASEWICH, E.R. 1981. Time Sequence Analysis in Geophysics, third edition: University of Alberta Press, Edmonton. KANASEWICH, E.R., and CHIU, S.K.L. 1985. Least-squares inversion of spatial seismic refraction data. Bulletin of the Seismological Society of America, 75: 865-880. KEEN, C.E. and HYNDMAN, R.D. 1979. Geophysical review of the continental margins of eastern and western Canada. Canadian Journal of Earth Sciences, 16: 712-749. KLEMPERER, S.L., HAUGE, T.A., HAUSER, E.C., OLIVER, J.E., and POTTER, C.J. 1986. The Moho in the northern Basin and Range province, Nevada, along the COCORP 40°N seismic reflection transect. Geological Society of America Bulletin, 97: 603-618. KLEMPERER, S.L. 1987. A relation between continental heat flow and the seismic reflectivity of the lower crust. Journal of Geophysics, 61: 1 - 11. LEWIS, T.J., BENTKOWSKI, W.H., DAVIS, E.E., HYNDMAN, R.D., SOUTHER, J.G. and WRIGHT, J.A. 1988. Subduction of the Juan de Fuca Plate: thermal consequences. Journal of Geophysical Research, 93: 15207- 15225.  95  LEWIS, T.J., BENTKOWSKI, W.H., and HYNDMAN, R.D. 1992. Crustal temperatures near the Lithoprobe Southern Cordillera Transect. Canadian Journal of Earth Sciences, 29: 1197-1214. McLEAN, N.A. 1993. Velocity structure from seismic refraction across the Coast Plutonic and Intermontane Belts. M.Sc. Thesis, University of Victoria. (In preparation.) LUTTER, W.J., NOWACK R.L., and BRAILE, L.W. 1990. Seismic imaging of upper crustal structure using travel times from the PASSCAL Ouachita experiment. Journal of Geophysical Research, 95: 4621-4631. McMECHAN, G.A., and SPENCE, G.D. 1983. P-wave velocity structure of the Earth's crust beneath Vancouver Island. Canadian Journal of Earth Sciences, 20: 742-752. MONGER, J.W.H. 1992. Canadian Cordilleran Tectonics: from Geosynclines to Crustal Collage. Canadian Journal of Earth Sciences. (In press.) MONGER, J.W.H., and PRICE, R.A. 1979. Geodynamic evolution of the Canadian Cordillera — progress and problems. Canadian Journal of Earth Sciences, 16: 770-791. MONGER, J.W.H., PRICE, R.A., and TEMPELMAN-KLUIT, D.J. 1982. Tectonic accretion and the origin of the two major metamorphic and plutonic welts in the Canadian Cordillera. Geology, 10: 70-75. MONGER, J.W.H. and JOURNEAY, J.M. 1992. A field guide to accompany the Penrose Conference on the Tectonic evolution of the Coast Mountains orogen, May 1992. 97 pp. MOONEY, W.D. and BROCHER, T.M. 1987. Coincident seismic reflection/refraction studies of the continental lithosphere: a global review. Reviews of Geophysics, 25: 723-742. NAFE, J., and DRAKE, C. 1963. Physical properties of marine sediments. In The Sea. Edited by M.N. Hill, pp. 794-828. PAKISER, L.C., and BRUNE, J.N. 1980. Seismic models of the root of the Sierra Nevada. Science, 210: 1088-1094. PHADKE, S., and KANASEWICH, E.R., 1990. Seismic tomography to obtain velocity gradients and three-dimensional structure and its application to reflection data on Vancouver Island. Canadian Journal of Earth Sciences, 27: 104-116. RIDDIHOUGH, R.P. 1979. Gravity and structure of an active margin — British Columbia and Washington. Canadian Journal of Earth Sciences, 16: 350-363. RIDDIHOUGH, R.P. and SEEMAN, D.A. 1982. Juan de Fuca Plate Map: JFP-8 Gravity Anomaly, Pacific Geoscience Center, scale 1:2 000 000. RODDICK, J.A. 1983. Geophysical review and composition of the Coast Plutonic Complex, south of latitude 55°N. In Circum-Pacific Plutonic Terranes. Edited by J.A. Roddick. Geological Society of America, Memoir 159, pp. 195-211. ROGERS, G.C., SPINDLER, C. and HYNDMAN, R.D. 1990. Seismicity along the Vancouver Island Lithoprobe Corridor, in Proceedings of the Project Lithoprobe: Southern Canadian Cordillera Transect Workshop, Calgary, 3-4 March 1990,166-169.  96  SHERIFF, R.E. and GELDART, L.P. 1982. Exploration Seismology, v.1, Cambridge University Press, New York. SPENCE, G.D., CLOWES, R.M., and ELLIS, R.M. 1985. Seismic structure across the active subduction zone of western Canada. Journal of Geophysical Research, 90: 6754-6772. SPENCER, C., ASUDEH, I., and COTE, T., 1989. SEGY-LDS Version 2.0 Format Reference Document, Manual Version 1.00, Geological Survey of Canada, Ottawa. 17 pp. STEINHART, J.S. 1967. Mohorovitie discontinuity. In International Dictionary of Geophysics, edited by S.K. Runcorn, Oxford, vol.2, pp. 991-994. VAN DER HEYDEN, P. 1992. A middle Jurassic to early Tertiary Andean-Sierran arc model for the Coast Belt of British Columbia. Tectonics 11: 82-97. VARSEK, J.L., COOK, F.A., CLOWES, R.M., JOURNEAY, J.M., MONGER, J.W.H., PARRISH, R.R., KANASEWICH, E.R., and SPENCER, C.S. 1992. LITHOPROBE crustal reflection structure of the southern Canadian Cordillera II: Coast Mountains transect. Tectonics. (In press.) WARNER, M. 1990. Basalts, water, or shear zones in the lower continental crust? Tectonophysics, 173: 163-174. WEBRING, M. 1985. SAKI: A Fortran program for generalized linear inversion of gravity and magnetic profiles. US Geological Survey, Open File Report 85-122. WEVER, T., TRAPPE, H., and MEISSNER, R. 1987. Possible relations between crustal reflectivity, crustal age, heat flow, and viscosity of the continents. Annales Geophysicae, 5B: 255-266. WHEELER, J.O. and McFEELY, P. (comp.) 1991. Tectonic Assemblage Map of the Canadian Cordillera and adjacent parts of the United States of America, Geological Survey of Canada, Map 1712A, scale 1:2 000 000. WHEELER, J.O., BROOKFIELD, A.J., GABRIELSE, H., MONGER, J.W.H., TIPPER, H.W. and WOODWORTH, G.J. (comp.) 1991. Terrane Map of the Canadian Cordillera, Geological Survey of Canada, Map 1713A, scale 1:2 000 000. WICKENS, A.J. 1977. The upper mantle of southern British Columbia. Canadian Journal of Earth Sciences, 14: 1100-1115. ZELT, B.C., ELLIS, R.M., and CLOWES, R.M. 1990. SCoRE '89: The Southern Cordillera Refraction Experiment Description of Dataset. LITHOPROBE Report No. 20, LITHOPROBE Secretariat, University of British Columbia, Vancouver, 38 pp. ZELT, B.C., ELLIS, R.M., CLOWES, R.M., KANASEWICH, E.R., ASUDEH, I., LUETGERT, J.H., HAJNAL, Z., IKAMI, A., SPENCE, G.D., and HYNDMAN, R.D. 1992a. Crust and upper mantle velocity structure of the Intermontane belt, southern Canadian Cordillera. Canadian Journal of Earth Sciences, 29: 1530 1548. ZELT, B.C., ELLIS, R.M., and CLOWES, R.M. 1992b. Crust and upper mantle velocity structure in the eastern Insular and southernmost Coast belts, British Columbia, Canada. Canadian Journal of Earth Sciences. (Submitted.) -  97  LUX, C.A., and ELLIS, R.M., 1988. Practical and efficient ray tracing in 2—D media for rapid traveltime and amplitude forward modelling. Canadian Journal of Exploration Geophysics, 24: 16-31. ZELT, C.A., and SMITH, R.B., 1991. Seismic traveltime inversion for 2—D crustal velocity structure. Geophysical Journal International, 108: 16-34.  Appendix A Terrane Descriptions Terranes of the Insular Superterrane  The Alexander terrane is composed of Upper Proterozoic to Triassic volcanic and sedimentary rocks in various depositional settings (including ocean arc, back arc, platform, rift, trough, and offshelf) and comagmatic intrusions (Monger and Journeay 1992). Wrangellia is a submarine arc terrane consisting of Middle to Late Triassic tholeiitic basalts and associated mafic intrusions, calcareous sedimentary rocks of similar age, and Early to Middle Jurassic terrigenous clastics and volcanics (Coney et al. 1980). Terranes of the Intermontane Superterrane  The following terranes of the Intermontane Superterrane are of the island arc type. The Stikine terrane includes interbedded arc volcanics, with ages ranging from Lower Devonian  to Permian, overlain unconformably by Upper Triassic arc volcanics and granitic rocks, and unconformably overlying Lower Jurassic andesite arc volcanics and intercalated sedimentary rocks (Monger and Journeay 1992). The Bridge River terrane, where metamorphism is not predominant, is characterized by extremely disruptive ribbon chert, argillite and basalt, subordinate rocks including siltstones and greywackes, as well as mafic intrusives. The age of the terrane ranges from Mississippian to late Middle Jurassic, and metamorphic grade is mostly sub-greenschist to greenschist (Monger and Journeay 1992). The Cache Creek terrane is partly disrupted and partly coherent, and is composed of Permian to Lower and  Middle Jurassic chert and pelite, with melanges of basalt, ultramafic rocks, and carbonates of Middle Pennsylvanian and Permian age. The Quesnel terrane is comprised of Upper Triassic to Early Jurassic arc volcanics, granitic and alkaline intrusions, and clastic rocks; overlying these unconformably are the clastic rocks (ca. 190-160 Ma) of the Ashcroft Formation (Monger and Journeay 1992). 98  Terranes of the Coast Belt The Cadwallader terrane is an island arc terrane consisting of Permian-aged pioneer  greenstone, diorite, trondhjemite, gabbro, and alpine-type ultramafic rocks, overlain by Upper Triassic arc volcanics, carbonates, and distinctive clast conglomerates (of the Hurley Formation). These rocks are again overlain, by Jura-Cretaceous elastic rocks (Monger and Journeay 1992). The Shuksan terrane is an oceanic terrane which includes Upper Triassic and Lower Jurassic oceanic crust and sediments metamorphosed to greenschist and blueschist facies and Jurassic near-arc marginal basin crust and sediments. The Harrison terrane features Middle Triassic cherry argillites and mafic volcanics overlain by thick Lower Jurassic andesitic and dacitic volcanics. Strata of sedimentary and andesitic volcanic rocks of Middle and Late Jurassic age are also present within this island arc terrane (Monger and Journeay 1992).  99  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items