Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Structure of the Cascadia subduction zone off Vancouver Island : new evidence from seismic refraction… Wang, Carl X. 1997

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

Item Metadata

Download

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

Full Text

STRUCTURE OF T H E CASCADIA SUBDUCTION ZONE OFF VANCOUVER ISLAND : NEW EVIDENCE FROM SEISMIC REFRACTION DATA By Carl X. Wang B. Sc. (Mech. Engineering) Beijing Institute of Technology A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R OF SCIENCE in THE FACULTY OF GRADUATE STUDIES EARTH AND OCEAN SCIENCES - GEOPHYSICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 1997 © Carl X. Wang, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. EARTH AND OCEAN SCIENCES - GEOPHYSICS The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract The Cascadia subduction zone forms the present convergent margin between the North American plate and the oceanic Gorda, Juan de Fuca and Explorer plates from ap-proximately 40°N to 50°N. This thesis examines a subset of the offshore-onshore seis-mic refraction/wide-angle reflection data recorded during the RECAP'89 (Reflection-refraction Experiment: Cascadia Accretionary Prism). RECAP'89 was a "piggyback" experiment in which The University of British Columbia deployed four continuous record-ing seismographs at eight land sites on the west coast of Vancouver Island to record airgun shots fired during an offshore multichannel reflection survey acquired by the Geological Survey of Canada. 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. With the aid of interpretation of corresponding reflection data, velocity structures were developed for Line 89-02/03 and 89-09 in Cascadia subduction zone. The principal tectonic elements for which velocity structures are determined are (i) the Pacific Rim with an average velocity of 5.0-6.2 km/s, and Crescent terranes with an average velocity of 4.5-6.4 km/s, which He along the coast against and beneath Wrangellia, (ii) the Accreted Wedge with an average velocity of 4.1-5.0 km/s below Tofino Basin, (iii) the subducting Juan de Fuca plate with an average velocity of 6.1-7.0 km/s. The new velocity models are consistent with the gravity field constraints. The final velocity model of Line 89-02/03 indicates a possible position for the updip Hmit of the seismic rupture zone for a potential great thrust earthquake in the area. This new improved velocity structure would also help to determine more accurate earthquake locations for the region. n Table of Contents Abstract 1 1 List of Tables vi List of Figures vii Acknowledgement 1 X 1 INTRODUCTION 1 1.1 Tectonic History 1 1.2 Previous Geophysical and Geological Studies 3 1.3 Experiment and Objectives 4 2 DATA ACQUISITION AND PROCESSING 7 2.1 The Multichannel Reflection Experiment 7 2.2 The RECAP'89 Refraction Experiment 8 2.3 Processing of Refraction Data 10 2.3.1 Initial Data Reduction 10 2.3.2 Noise Reduction 11 2.3.3 Predictive Deconvolution 11 2.3.4 Amplitude Corrections 13 3 DATA ANALYSIS AND INTERPRETATION PROCEDURE 16 3.1 Interpretation of Reflection Data 16 m 3.1.1 Line 89-02/03 16 3.1.2 Line 89-09 17 3.2 Initial Analysis of Refraction Data 21 3.2.1 Data Characteristics and Quality 21 3.2.2 Travel Time Picking 26 3.3 Interpretation Methods of Refraction Data 27 3.3.1 The Modelling Algorithms 27 3.3.2 Traveltime Inversion Procedure 29 3.3.3 The Modelling Techniques 33 4 MODELING AND INTERPRETATION OF REFRACTION DATA 36 4.1 Modelling the Refraction Data 36 4.1.1 Initial Model 36 4.1.2 Line 89-02/03 37 4.1.3 Line 89-09 44 4.2 Primary Features of the Final Models 53 4.2.1 Line 89-02/03 53 4.2.2 Line 89-09 55 4.3 Non-uniqueness of the Model 56 4.4 Model Resolution and Parameter Uncertainty 57 5 DISCUSSION AND CONCLUSION 60 5.1 Testing of Consistency with Gravity Data and Comparison with Density Structures 60 5.1.1 Comparison between Line 89-02/03 and Line IAS 62 5.1.2 Comparison between Line 89-09 and Line 2A 67 5.2 Does the Velocity Model Relate to Potential Great Thrust Earthquakes? 68 iv 5.3 Conclusions References List of Tables 2.1 RECAP '89 Locations of Onshore Refraction Stations 9 4.1 Estimated lateral resolution and absolute uncertainty 59 v i List of Figures 1.1 Regional tectonic setting of the study area 2 1.2 Location Map 6 2.1 Amplitude spectrums before and after filtering 12 2.2 Data before and after deconvolution 14 3.1 Multichannel reflection Une 89-02/03 18 3.2 Multichannel reflection line 89-09 20 3.3 Observed record of line 89-02/03 at 2A 22 3.4 Observed record of line 89-02/03 at 2B 23 3.5 Observed record of line 89-09 at 4A 24 3.6 Observed record of line 89-09 AT 4B. . 25 4.1 Color contour plot of velocities beneath line 89-02/03 38 4.2 Observed record and synthetic of line 89-02/03 at 2A 39 4.3 Ray tracing and traveltimes of line 89-02/03 at 2A 40 4.4 Observed record and synthetic of Une 89-02/03 at 2B 41 4.5 Ray tracing and traveltimes of Une 89-02/03 at 2B 42 4.6 Color contour plot of velocities beneath Une 89-09 45 4.7 Observed record and synthetic of Une 89-09 at 4A 46 4.8 Ray tracing and traveltimes of Une 89-09 at 4A 47 4.9 Observed record and synthetic of Une 89-09 at 4B 48 4.10 Ray tracing and traveltimes of Une 89-09 at 4B 49 vn 4.11 Model nodes location for both lines 50 4.12 Ray coverage for both lines 51 5.1 Bouguer gravity anomaly map 61 5.2 Laboratory measurements of P-wave seismic velocity and density in rocks, 63 5.3 Gravity modeling of Line IAS and 89-02/03 65 5.4 Gravity modeling of Line 2A and 89-09 66 5.5 Location map and model profile of thermal analysis 69 5.6 Seismicity cross section 74 vm Acknowledgement I wish to thank my supervisor Dr. Ron Clowes and Department Head Dr. Robert M. Ellis for giving me the opportunity to come to this lovely country, and to Dr. Ron Clowes for providing the funding to allow this study to proceed and guiding me all the way through this study. I also would like to thank Brad Isbell and Sonya Dehler for the hard work of collecting the field data. I especially appreciate the persistent effort of Brad Isbell in data processing to bring the analog data to digital files. My thanks also go to Dr. George Spence at the University of Victoria for providing the shot origin times. I would also like to thank my teammates: manager John Amor, Phil Hammer, Dierdre O'Leary, Mike Perz and Dave Baird for their help. The other departmental members have made this a great place to study - I am very grateful for all of them for their help and friendship. Financial support for my student stipend and analysis costs derived from NSERC Research Grant OGP0007707 to Dr. Ron Clowes. i x Chapter 1 I N T R O D U C T I O N The Cascadia subduction zone forms the present convergent margin between the North American plate and the oceanic Gorda, Juan de Fuca and Explorer plates from approxi-mately 40°N to 50°N. The potential for a devastating megathrust earthquake along the northern part of the subduction zone (Rogers 1988) has focused attention on the need to fully understand the structure and tectonic history of this convergent margin. 1.1 Tectonic History The modern plate-tectonic regime of the continental margin off southwestern Canada and northwestern United States is dominated by the motions of three main lithospheric plates: the large Pacific and American plates and the intervening Juan de Fuca plate system which converges with North America at a present rate of 40-50 mm-a - 1 (Figure 1.1). The plate regime for the past 10 Ma, over which time convergence in the study area has not changed in a significant way, is well known. However, the earlier regime is much less certain. In the early Tertiary there may have been a triple junction south of Vancouver Island with transform motion to the north and subduction to the south. The North America plate has been advancing westward toward the Pacific at a rate of about 22 mm-a - 1 for at least the past 10 Ma and probably at a similar rate through the Tertiary. Thus, the downgoing Juan de Fuca plate is being overridden by the North American continent. The formation of the Tertiary margin appears to have begun with the initiation of Figure 1.1: Regional tectonic setting of the study area (outlined by rectangle) within the Cascadia subduction zone off Vancouver Island. Solid triangles locate volcanoes. Chapter 1. INTRODUCTION 3 subduction in the Early Eocene. This was followed by a shift in the locus of under-thrusting to the west, probably as a result of a major reorganization of the plate-tectonic regime. Two small terranes, the Mesozoic metasedimentary Pacific Rim Terrane and the Eocene ophiolitic Crescent Terrane were emplaced against and beneath the Insular Superterrane, comprising the previously accreted volcanic Wrangellia terrane, at about that time. The modern accretionary complex has formed beneath and against the Crescent Ter-rane as a result of scraping off incoming sediments on the subducting Juan de Fuca plate. Subduction since the Eocene has been relatively simple, the continuous record of accretion showing no evidence for the arrival of allochthonous terranes. The T o r i n o Basin covers most of the continental shelf off Vancouver Island to near the shelf edge. This forearc basin was formed by deposition of Eocene to recent marine clastic sediments over the Pacific Rim and Crescent terranes and the inner portion of the modern accreted sedimentary wedge. A narrow zone of outcrop occurs along the south coast of Vancouver Island as the Carmanah Group. Torino basin was the site of petroleum exploration in the 1960s and a number of offshore wells were drilled and logged (Shouldice 1971). However, there has been little interest in the prospectiveness of Tofino basin since then. Hyndman et al. (1990) and Hyndman (1995) provide recent reviews of the structural and tectonic history of the northern Cascadia subduction zone. 1.2 Previous Geophysical and Geological Studies A n extensive geoscience data base has been compiled over the past three decades of geophysical and geological investigation of this region. Regional coverage is provided by gravity and magnetic data and geological mapping, while local measurements include Chapter 1. INTRODUCTION 4 heat flow, borehole, magnetotellurics, and seismicity. Both seismic reflection and seismic refraction lines provide profile coverage across and along the margin. Of relevance to this study are the ocean bottom seismograph refraction results of Waldron et al. (1989) who provided a velocity model from the deep ocean to the shelf edge, and the offshore-onshore refraction results of Spence et al. (1985) and Drew and Clowes (1990) who interpreted velocity structural models from the ocean basin across Vancouver Island to the mainland. Of particular note is the fact that the logistics of the experiment provided few velocity constraints from the shelf edge to western Vancouver Island and from the sediments of Torino basin to the top of the subducting oceanic plate. This part of the margin is the focus of the present study. Hyndman et al. (1990) and Hyndman (1995) provide a summary of much of the data available for studying of the northern Cascadia margin, interpretation of these data and their significance for understanding tectonic evolution. Dehler and Clowes (1992) present an overview of the major geological and geophysical studies used in the detailed modelling and analysis of the gravity and magnetic data off southwestern Vancouver Island. Most recently, Clowes et al. (1996) present a series of two-dimensional density cross-sections along the margin and a block model of density structure developed through the integration of numerous geophysical and geological data and interpretations. The model illustrates variations in extent, structure and properties of the major tectonic components along most of the Vancouver Island margin and provides a framework within which to examine previously uninterpreted seismic data. 1.3 Experiment and Objectives This thesis examines a subset of the offshore-onshore seismic refraction/wide-angle reflec-tion data recorded during the RECAP'89 (Reflection-refraction Experiment: Cascadia Chapter 1. INTRODUCTION 5 Accretionary Prism). RECAP '89 was a "piggyback" experiment in which The University of British Columbia deployed four continuous recording seismographs at eight land sites on the west coast of Vancouver Island to record airgun shots fired during an offshore multichannel reflection survey acquired by the Geological Survey of Canada (Figure 1.2). The primary objectives of RECAP '89 are to provide velocity-depth structural sections for a part of the margin, particularly below the continental shelf, for which there is a lack of such information, thereby providing information useful for better locations of earth-quakes in the region; and to relate the structural sections to the tectonic development of the margin. The principal tectonic elements for which velocity structures are determined are (i) upper and western Wrangellia, (ii) the Pacific Rim and Crescent terranes which He along the coast against and beneath Wrangellia, (iii) the accreted wedge below Tofino Basin, and (iv) the subducting Juan de Fuca plate (Figure 1.2). 1 2 8 - 1 2 6 . 1 2 4 ° Figure 1.2: Location Map. Center straight lines show 1989 offshore reflection profiles recorded at seismograph sites on land (stars). Bold lines and corresponding bigger stars indicate the locations of present study (89-02/03 at 2A and 2B; 89-09 at 4A and 4B). Terranes and other tectonic features are indicated. Drillholes in Tofino Basin are shown by the "dryhole" symbols. Chapter 2 D A T A ACQUISITION A N D P R O C E S S I N G 2.1 The Multichannel Reflection Experiment From September 1 to 15, 1989, a marine multichannel seismic survey was shot off the west coast of Vancouver Island (Figure 1.2). A total of 13 lines (722 km) of reflection data was obtained over the Tofino Basin, the continental slope, and the deformation front of the accretionary wedge. The data were collected by Digicon (Canada) Inc., using the vessel M / V Geo Tide equipped with Digicon's DSS-240 cable system. The data were processed by a commercial contractor, Geophoto Inc. Both data acquisition and processing were under contract to the Geological Survey of Canada as part of site surveys for the scientific international Ocean Drilling Program (Spence et al., 1991) The objective of the reflection program was to determine the structural architecture and tectonic history of the Vancouver Island continental margin. Of particular interest were (1) the regional structure and stratigraphy of the Tofino Basin and of the deeper structures within the accretionary prism which provide basin controls; (2) the three-dimensional variation in fold and thrust structures at the deformation front; (3) physical properties related to the shallow bottom-simulating reflector beneath the lower slope, previously interpreted as the base of a methane hydrate layer; and (4) the nature of the basal thrust or detachment that could produce a large earthquake. The 144 channel system recorded data for 14 s at a sampling rate of 4 ms. The 3600 m cable, towed at a depth of 12 m, was configured at a group interval of 25 m, with 14 7 Chapter 2. DATA ACQUISITION AND PROCESSING 8 hydrophones in each group. The seismic source was a tuned airgun array fired at 50 m intervals to provide 36-fold data. The array was made of 4 strings of 10 guns each (3 active strings and one spare), with a total airgun volume of 7820 in3 operating at 1900 p.s.i. and a depth of 10 m. The offset from the center of the gun array to the center of the first group was 183 m. Navigation for the survey was by Starfix augmented by satellite navigation. The data were processed in Calgary by Geophoto Inc. Processing procedures in-cluded: (1) true amplitude recovery, with the application of exponential gain and spher-ical divergence corrections; (2) f-k (frequency-wavenumber) velocity filter to suppress low velocity backscatter; (3) zero-phase source designature plus f-k demultiple; (4) derivation of stacking velocities from semblance analyses for sub-bottom times less than about 5 s, combined with a general refraction velocity model for the area for greater depths; (5) gapped deconvolution after stacking; (6) f-k migration, using hand-smoothed N M O func-tions with 10 percent velocity reduction to prevent over-migration. Spence et al. (1991) provide a summary of the experiment and preliminary interpretation of the data. 2.2 The RECAP'89 Refraction Experiment The airgun source from the offshore reflection survey was used in a piggyback experiment to record refraction and wide-angle reflection data on sites on Vancouver Island both inline and oblique to the reflection profiles which are perpendicular to the strike of the margin (Figure 1.2). The onshore refraction stations are numbered 1A, IB, 2A, 2B, 3A, 3B, 4A, and 4B. Isbell (1989) provides a detailed summary of the U B C field experiment. A l l of the onshore stations (Table 2.1) were deployed on bedrock. Most of the stations were very secluded but several sources of seismic noise were apparent. Most of these sources were localized and discontinuous, e.g. wind, traffic, rock-fall, creek, waves etc. Chapter 2. DATA ACQUISITION AND PROCESSING Onshore • C o o r d i n a t e s E l e v Map D i s t a n c e S t a t i o n N o r t h i n g E a s t i n g (m) (92) From L i n e L a t Long Zone 10 1A 5390050 392050 170 C/9 30.4 KM 48.656° -124.466° (1.2 S) L i n e i 48°39'22" -124°27'57" IB 5404300 413500 620 C/16 56.1 KM 48.788° -124.178° (1.4 N) 4-8'°47 ' 16" -124°10'39" 2A 5406160 352130 330 C/14 19.3 KM v 48.793° -125.013 (0.2 S) L i n e n 48°47'34" -125° 0*48" 2B A 5418450 372050 630 C/15 42.6 KM/ 48.908° -124.746° (0.3 S) 48°54'28" -124°44'46" 3A 5432000 313600 85 F/4 18 . 4 KM 49.015° -125.549° (0.2 N) L i n e c 49° 0'53" -125°32'57" 3B D 5450650 328300 235 F/3 41.3 KM 49.187° -125.356° ! 49°11'12" -125°21 123" (7.4 N) Zone 9 4A 5471800 698725 0 E/8 12.6 KM 49.368° -126.263° (0.4 N) L i n e Q 49°22' 6" -126°15'46" 4B 5481200 714325 5 E/8 30.7 KM 49 . 448° -126.043° (1.8 S ) 49°26'51" -126° 2'35" Table 2.1: R E C A P '89 Locations of Onshore Refraction Stations. From Isbell (1989) Chapter 2. DATA ACQUISITION AND PROCESSING 10 Four FM analog slow speed (15/160 ips) 7-channel seismic recording systems with 1 Hz vertical and horizontal geophones, the outputs being recorded at two different gain settings to extend the dynamic range, were used. Filters on the amplifiers were set for a 1-12 Hz pass band, consistent with the slow speed recording. For timing, the radio time code WWVB was recorded on one channel, this time base also being related to the times of the airgun shots. The recorded analog data were reduced, digitized and compiled in SEGY format. The focus of this study is the analyses of data from lines 89-02/03 at stations 2A/2B, and 89-09 at stations 4A/4B. 2.3 Processing of Refraction Data 2.3.1 Initial Data Reduction The initial data reduction was carried out by Brad Isbell, a research assistant working with Dr. Ron Clowes, after the acquisition experiment. More details about this process can be found in Isbell (1990). Data reduction for the refraction data was divided into four phases: (1) Looking at the FM analog tapes using the Sanborn tape system to determine data quality, the precise locations of the data on each tape, unusual problems with the data, and the appropriate seismic channel to be digitized for each recorded line. (2) acquiring the navigation and shot time data, correcting errors in these data sets, and combining these data sets into one single file for each line which includes both navigation data and the associated WWVB shot times in milliseconds. (3) checking the digitizing system and making any necessary adjustments and modific-ations to the existing system, including modifying the flutter compensation procedure and existing computer programs. (4) digitizing the tapes and storing the data in SEGY format. Backup tapes are also Chapter 2. DATA ACQUISITION AND PROCESSING 11 made at the same time. 2.3.2 Noise Reduction A commercial seismic processing software package, INSIGHT from former Inverse Theory and Application (ITA) Ltd. (now Landmark/ITA), installed on a Sun Microsystem SPARCstation 10, provided the base for most of the application processing in this study. 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 6-12 Hz range (Figure 2.1a), whereas that of the noise was mainly outside of this range. A 5-15 Hz band-pass filter was applied to eliminate low-frequency micro-seismic noise and high-frequency instrumentation noise (Figure 2.1b). These filtered sections were used concurrently with the respective unfiltered sections to make travel time picks. Data Summing Additional noise reduction was achieved by summing or binning the traces within 150 m bin widths using a velocity of 6.5 km/s, which is the reducing velocity for plotting traces. This velocity is chosen by testing with an aim to boost secondary arrivals which have an average apparent velocity close to 6.5 km/s. 2.3.3 Predictive Deconvolution Most of the wide-angle land recordings show a prominent ringing of the signals (Figure 2.2a), which often makes the identification of later arrivals difficult. It is therefore nec-essary to deconvolve the traces to remove the ringing and sharpen the arrivals. Most of the ringing is probably due to seabed multiples which also show up on the corresponding Chapter 2. DATA ACQUISITION AND PROCESSING 12 Figure 2.1: (a) Amplitude spectrum of data before filtering, (b) Amplitude spectrum of data with a bandpass filter of 5-15 Hz applied. Chapter 2. DATA ACQUISITION AND PROCESSING 13 reflection data (Spence et al., 1990). If so, predictive deconvolution using a short op-erator should work because seabed multiples are generally short period. However, tests using short operators and short prediction lengths (100 ms) yielded poor results. There are two factors limiting this kind of deconvolution. First, multiples are aperiodic for non-zero offsets and this indicates most of the multiples are not strictly predictable. Second, these data are extremely bandlimited which will substantially degrade the performance of predictive deconvolution with short prediction lengths. After extensive testing of prediction length, operator length, and percentage of noise to be added, I found that that using a prediction length of 150 ms, an operator length of 800 ms with 2 percent of added noise gave good results (Figure 2.2b). After deconvo-lution, the ringing is substantially decreased and a higher resolution is achieved. Thus the identification of phases, especially secondary phases, becomes more reliable. The predictive deconvolution program in the Landmark/ITA INSIGHT processing package was used. 2.3.4 Amplitude Corrections Amplitude corrections were applied to the data to account approximately for decrease in amplitudes as the waves propagate to larger distances. After some testing to make the amplitude of the whole data section balanced, each trace was scaled by a factor proportional to d15, where d is the shot-receiver distance. This process enhanced the amplitudes of the far-offset traces with respect to those near the source. A point to be noted here is that the trace normalized data can only be used to look at relative amplitudes between phases on the same trace, whereas amplitude corrected data enable relative amplitude with offset to be compared with those from a theoretical calculation. This is important because amplitudes are sensitive to velocity gradients. The final step in preparing the data for interpretation was editing out traces which Chapter 2. DATA ACQUISITION AND PROCESSING 14 Figure 2.2: Data before and after deconvolution. (a) Example of refraction data before predictive deconvolution. (b) Same data after predictive deconvolution with a predictive length of 150 ms, an operator length of 800 ms and 0.02 of noise added. The reduction velocity used is 6.5 km/s. Chapter 2. DATA ACQUISITION AND PROCESSING 15 were so noisy that they obscured the signal visible on adjacent traces when the data were displayed with the amplitude corrections described above. V Chapter 3 D A T A ANALYSIS A N D I N T E R P R E T A T I O N P R O C E D U R E 3.1 Interpretation of Reflection Data Although the focus of this study is the modeling and interpretation of the refraction data, a very brief description of the interpretation of the offshore reflection data is presented here since this interpretation will be used as a guide for the model construction based on refraction data. More detailed analysis and interpretation of the reflection data can be found in Spence et al. (1991). 3.1.1 Line 89-02/03 Multichannel seismic line 89-02/03 was recorded in the southern part of the study area (Figure 1.2) from the deep ocean, across the deformation front, up the slope and across Tofino basin underlain by the accreted wedge, Crescent and Pacific Rim terranes. The seismic reflection profile of line 89-02/03 is shown in Figure 3.1 with variable-density display. The clearly layered and strongly reflective sediments at the seaward end of the profile (around 140 km mark) represent the undeformed deep ocean basin; east of that is the Deformation Front (DF). In this area (120-140 km), many folds and thrust belts, which can not be seen clearly in this section, were found on other adjacent reflection lines (Spence et al., 1991). The strong basement reflection underneath the deep ocean basin represents the top of the subducting oceanic Juan de Fuca plate. 16 Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCED URE 17 The sedimentary layer becomes severely deformed and therefore only weakly reflective east of the ocean basin (95-120 km). This is the deformed accreted wedge sediments below which the signature of the subducting oceanic plate can not be continuously seen due to the highly deformed upper crust, although pieces of it appear at shallower times than beneath the ocean basin because of the velocity pull-up effect in this time section. The gap between Line 89-02 and Line 89-03 also made it difficult to follow these events. Further to the east, partly coherent, moderate to high-amplitude reflections within the Tofino basin (40-80 km in distance) are in marked contrast to the weakly reflective to non-reflective underlying continental crust. The very high-amplitude reflections in the lower part of the basin fill at some locations are probably generated by interbedded elastics and volcanics (Spence et al, 1991). One deep dipping event around 5.0-7.0 s at 35-45 km distance may represent the Crescent thrust or lower boundary of the igneous Crescent Terrane. Another shallow dipping event around 3.5 s at 40 km distance is probably the top of the Crescent Terrane which is right beneath the base reflection of the Tofino basin. No clear reflection corresponds to the Pacific Rim Terrane which exists underneath the sedimentary layer at the east end of the profile. A set of deeper eastward-dipping reflections beneath the continental margin (from 25 to 80 km in distance) is interpreted as the down-going oceanic plate. The interpretation of the lower boundary of the upper sedimentary layer (Sed in Figure 3.1), which makes it easier to set up the velocity model, and the top of the subducting oceanic plate are used in the construction of the initial velocity model for the modeling of the refraction data of Line 89-02/03. 3.1.2 Line 89-09 Multichannel seismic line 89-09 was recorded in the northern part of the study area (Figure 1.2) from the deep ocean, across the deformation front, up the slope and across Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 18 Distance (km) Figure 3.1: Multichannel reflection line 89-02/03. (1) and (2) are enlargements of portions of the sedimentary layer. Sed, upper sedimentary boundary used in the modeling of refraction data; JdF, location of the top of the Juan de Fuca plate used in the modeling of refraction data; Mult, multiples. All distance scales in this study are shot-receiver distances from the first onshore receiver site at each line (2A or 4A) unless otherwise specified; DF, Deformation Front; CT, Crescent Terrane; PRT, Pacific Rim Terrane. Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 19 Tofino basin underlain by the accreted wedge, Crescent and Pacific Rim terranes. The seismic reflection profile of fine 89-09 is shown in Figure 3.2 with variable-density display. Similar to Line 89-02/03, the clearly layered and strongly reflective sediments at the seaward end of the profile (around 120 km mark) are the undeformed deep ocean basin, and east of that is the Deformation Front (DF). In this area (115-125 km), many folds and thrust belts, which can not be seen clearly in this section, were found on other adjacent reflection lines (Spence et al., 1991). The strong basement reflection underneath the deep ocean basin represents the top of the subducting oceanic Juan de Fuca plate. Also similar to Line 89-02/03, the sedimentary layer becomes deformed and therefore only weakly reflective east of the ocean basin (70-105 km). This is the deformed accreted wedge sediments below which the signature of the subducting oceanic plate can be seen more continuously than on Line 89-02/03 because the upper crust is less deformed here. These events appear at shallower times than beneath the ocean basin because of the velocity pull-up effect in the time section. Further to the east, partly coherent, moderate to high-amplitude reflections within the Tofino basin (35-70 km in distance) are in marked contrast to the weakly reflective to non-reflective underlying continental crust. The very high-amplitude reflections in the lower part of the basin fill at some locations are probably generated by interbedded elastics and volcanics (Spence et al, 1991). No clear reflections corresponding to Crescent Terrane and Pacific Rim Terrane can be clearly identified. A set of deeper eastward-dipping reflections beneath the continental margin (from 15 to 70 km in distance) is interpreted as the down-going oceanic plate. Due to the fact that the subducting oceanic plate is shallower on this line, some reflection events below the oceanic plate may represent the oceanic moho (marked as O M in Figure 3.2). Similar to Line 89-02/03, the interpretation of the lower boundary of the upper sed-imentary layer (Sed in Figure 3.2), which makes it easier to set up the velocity model, Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 20 Distance (km) Figure 3.2: Multichannel reflection line 89-09. (1) and (2) are enlargements of portions of the sedimentary layer. Sed, 1 boundary used in the modeling of refraction data; Jdf, location of the top of the Juan de Fuca plate used in the modeling of refraction data; Mult, multiples; DF, Deformation Front; OM, Oceanic Mono of JdF crust. Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 21 and the top of the subducting oceanic plate are used in the construction of the initial velocity model for the modeling of the refraction data of Line 89-09. 3.2 Initial Analysis of Refraction Data 3.2.1 Data Characteristics and Quality In this section, a brief description of the data characteristics and quality for each identified phase is given. Acquisition of the field data resulted in a data set of 7044 traces (1761 traces after summing). The four record sections are displayed in Figures 3.3-3.6. The general quality of the data is good, especially for first arrival times, but some secondary arrivals are interfered with by high-level random noise and "ringing" (multiple cycles from the source and/or internal reflections). The branches of the phases, which were identified by close observation and preliminary forward modelling, are labelled in Figures 3.3-3.6. In total, four phases were identified on each section: event A , which represents the turning rays through the upper layers of the continental crust; event C, which includes reflections from the top of the oceanic plate; event M l , which denotes energy reflected from the oceanic Moho; and event M2, which represents energy refracted through the oceanic mantle. Because of low signal-to-noise ratio, instrument error, and survey limitation (offshore reflection survey), not all the phases are clearly observed on all record sections. The duration of each phase also varies from one section to another. Event A is clearly visible on all sections. But in section 4A and 4B, the second half of the phase is difficult to follow because of (1) instrument error (e.g. battery breakdown), and (2) severely deformed upper crust (i.e. with little coherent reflectivities), in which rays are heavily defracted and attenuated. The variations of sediment thickness are primarily responsible for the apparent velocity changes of these first arrivals. The Distance (km) 140 135 130 125 120 115 UO 105 100 95 90 85 80 75 70 65 60 55 50 45 40 35 30 25 20 Figure 3.3: Observed record of line 89-02/03 at station 2A. The reduction velocity used was 6.5 km/s, and a bandpass filter (5-15 Hz) has been applied. (1) and (2) are enlargements of certain phases of the data. Specified phases, identified by alphanumeric labels, are described in the text. Distance (km) Figure 3.4: Observed record of line 89-02/03 at station 2B. The reduction velocity used was 6.5 km/s, and a bandpass filter (5-15 Hz) has been applied. (1) and (2) are enlargements of certain phases of the data. Specified phases, identified by alphanumeric labels, are described in the text. K> Distance (km) 120 100 80 BO 40 20 Figure 3.5: Observed record of line 89-09 at station 4A. The reduction velocity used was 6.5 km/s, and a bandpass filter (5-15 Hz) has been applied. The framed section is an enlargement of certain phases of the data. Specified phases, identified by alphanumeric labels, are described in the text. to Distance (km) Figure 3.6: Observed record of line 89-09 at station 4B. The reduction velocity used was 6.5 km/s, and a bandpass filter (5-15 Hz) has been applied. The framed section is an enlargement of certain phases of the data. Specified phases, identified by alphanumeric labels, are described in the text. Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 26 amplitude of A is normally high, which indicates a large velocity gradient in the upper layer of the continental crust. Event C appears as a secondary arrival on all sections but with a low signal-to-noise ratio. It is quite difficult to follow on sections 4A and 4B (Figure 3.5 and 3.6). The amplitude of C is relatively high, but it is not prominent due to high noise level presumably associated with scattered energy. Event C (both travel time and amplitude) is the principal constraint for the velocity contrast across the top of the subducting oceanic plate. Event M l appears as a later secondary arrival on sections 2A and 2B differently from sections 4A and 4B. It is partly due to the gap between line 89-02 and 89-03. Another factor is the highly deformed upper crust sedimentary layer. It is much easier to follow on sections 4A and 4B than sections 2A and 2B, on which M l is hidden in the ringing "tail" of phase C . Some "dead" traces make it more difficult to follow on section 4B. Event M2 is visible on all sections with the best quality on section 4B. On section 4A, it has a low signal-to-noise ratio probably caused by higher local noise level. On sections 2A and 2B, it is somewhat difficult to follow because of the deformed sedimentary layer right under the shots. M2 has an apparent velocity of about 7.8-8.2 km/s which is reasonably consistent with the expected velocity of the upper oceanic mantle. Although the amplitude of M2 is low, it still indicates that there is some velocity gradient in the upper oceanic mantle. 3.2.2 Travel Time Picking Identifying and picking travel times of the phases described in the previous subsection in-volved looking at plots (of different scales) as well as computer screen presentations using ITA interactive software. The procedure followed for travel time picking was dependent upon the nature of the phase, as outlined below. Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 27 First, the A and M2 first arrivals were picked from large scale presentations of un-filtered traces using the interactive software package. For most parts of the sections, this was a relatively straightforward task because of the high quality of the data and the uniform nature of the waveforms. Each travel time pick was assigned an uncertainty dependent upon the difficulty of pinpointing the first motions. For phase M2 on section 4A, it was necessary to use a filtered section due to high noise level. The picking of phases C and M l involved more detailed inspection because of the low signal-to-noise ratio, variations in the waveforms, and the complicated interference among phases, their multiples and scattered energy. To enhance the phase coherency, different methods have been tried including various filtering, predictive deconvolution, plotting with different reducing velocities, and plotting in common maximum amplitude format. Through some preliminary modelling, it was found that further adjustment of some picks, even re-picking in some cases, was necessary because it might happen that instead of the proper arrival, some later phase was picked, especially for secondary phases C and M l . This fact required that the whole interpretation process should be at least partly interactive in order to enable the possibility of correcting the errors in picking even between iterations. 3.3 Interpretation Methods of Refraction Data 3.3.1 The Modelling Algorithms The modelling procedure applied to the refraction/wide-angle reflection data combined the following approaches: (1) the 2-D rapid ray-trace travel time and amplitude forward modelling algorithm of Zelt and Ellis (1988); Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 28 (2) the 2-D ray-trace travel time inversion algorithm of Zelt and Smith (1991). Both the inversion and the forward modeling algorithm are based on the numerical so-lution of the ray-tracing equations of zero-order asymptotic ray theory (Cerveny 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 uses the Runge-Kutta method with error control as suggested by Cerveny 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. The 2-D travel time forward modelling method (Zelt and Ellis 1988) was used to gen-erate a starting model because it is more interactive and easy to incorporate the structural and velocity information from the analysis of corresponding reflection data. The trial-and-error forward modelling also helped to identify the less obvious, later phases before including them in the inversion. Then the Zelt and Smith (1991) 2-D travel time inver-sion was chosen primarily because it represents a more objective and efficient procedure to fit the real data. In addition, it provides useful estimates of resolution, uncertainty and non-uniqueness of model parameters (boundary depths and layer velocities). Fur-thermore, the high density shot coverage and presence of both refracted and reflected crustal and upper mantle arrivals on all record sections make this dataset well suited for inversion. Again 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. 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. Since the structure and reliability of the final model are strongly dependent upon Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 29 the nature of the inversion algorithm, the aspects of the inversion procedure which are important in this study are discussed below. 3.3.2 Traveltime Inversion Procedure The travel time t between a source and receiver along a ray path L is given in integral form for a continuous velocity field v(x,z) as = / -F^—.dl (3.1) J v(x,z) t L The discrete form t = t~ (3-2) is used in practical applications, where U and V{ are the path length and velocity of the ith ray segment, respectively. Therefore, travel time is a linear combination of slowness (reciprocal of velocity), but travel time inversion is a non-linear problem since the ray path is velocity dependent. 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 (3.3) where A t is the travel time residual vector, A m is the model parameter adjustment vector, and A is a matrix of partial derivatives. A contains the elements dti/dm-i , where U is the ith observed travel time and rrij is the jth model parameter selected for inversion (i.e. either a velocity value or the z-coordinate of a boundary node). To solve equation (3.3), rays are traced through the current model during a particular iteration and A t and A are computed, with partial derivatives calculated analytically using the method described in Zelt and Smith (1991). This determines the model updated Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 30 vector A m , which is added to the current model for further iterations or for stopping if the stopping criteria are met. Crustal refraction/wide-angle reflection experiments often give rise to a system (3.3) which is overdetermined (actually mix-determined), unstrained, 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, Phadke and Kanasewich 1990). In solving the linearized system of equations (3.3), 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: (f>(Am) = \\AAm - At\\2 + e2 | | A m | | 2 (3.4) where e2 is a scalar trade-off parameter that is specified in order to set the relevant importance of the terms, and means the l2 norm. Extremizing <f> with respect to A m yields the linear algebraic equations: A m = [ATA + e2l) _ 1 AT At (3.5) Following the model suggested by Lutter et al. (1990), if the trade-off parameter is defined as e2 = Dof/cr 2 , where (Ti — the standard deviation used to estimate the uncertainty of the i t h travel time measurement, <Tj = the standard deviation used as an a priori estimate of the uncertainty of the j t h model parameter, D = an overall damping parameter usually equal to one, Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 31 then this becomes the Zelt and Smith (1991) damped least-squares solution to (3.3): A m = (ATCr1A + DC'1)'1 ATC;1^t (3.6) where C t = the estimated model covariance matrix defined by diag[of], and C m = the estimated model covariance matrix defined by diag[«r|]. The values of the standard deviation <7i 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 o-j 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 upon 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 number and placement may be determined by running a series of tests based on the following criteria: (1) the best trade-off between travel time residual and parameter resolution is obtained; (2) rays must be traced to all observations. Using these two 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; however, while the number of nodes Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 32 chosen could be optimized, in certain cases their placement could not. In these instances, the choice of velocity node positioning as governed by the above guidelines was somewhat arbitrary and sometimes according to the availability of a priori information (e.g., results of velocity analysis of reflection data). The resolution matrix R defined by Zelt and Smith (1991) is: (A^A + D C - ^ A ^ A (3.7) The diagonal elements of this resolution matrix range between 0 to 1. Non-zero 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. The criteria determining optimal parameterization lead directly to the stopping criteria: (1) the final model chosen is that which gives the desired trade-off between travel time residual and parameter resolution; and (2) the final model must permit rays to be traced to all observation points. In general, travel time residuals can be reduced by adding more model parameters; however, overall parameter resolution is concurrently reduced. Thus there is a need to achieve an appropriate trade-off between these two measures. It is also essential to real-ize that the R matrix gives the resolution values for the parameters of the model update vector A m ; this does not provide knowledge about the resolution of the parameters of the final model. Often it is not possible to trace rays to all observation points using a model obtained from an inversion. A common cause of this problem is that the model consists of too many parameters, so the data are unable to adequately constrain all model pa-rameters, causing unrealistic oscillations in layer velocities or boundary positions, which cause rays to focus or scatter. A final model must therefore be rejected if it is not possible to trace rays to most observation points. Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 33 3.3.3 The Modelling Techniques (a) Pick Errors and Uncertainties: When picking the arrivals on seismic sections, we can commit at least two principal errors which may distort the solution. First, it may happen that instead of the proper arrival, some later phase is picked (e.g. if the true arrival is hidden in the noise or in the multiples of other arrivals). Second, it may happen that especially for wide-angle events, the type of wave may be wrongly assigned. These facts entail the requirement that the modelling process should be at least partly interactive in order to have the possibility to remove these types of errors even between iterations. This is another reason to use the interactive travel time forward modelling method as a compensation in the process. Pick uncertainties are required to allow for the appropriate data weighting during travel time inversion. For most of the data used here, pick uncertainties ranged from 30-90 ms, much larger than the actual timing uncertainties associated with the data. Data with much higher than normal pick uncertainties are eliminated to ensure that a consistent set of data, i.e., one that can be fit within its uncertainties without introducing artificial lateral heterogeneities, is used for modelling. (b) Smooth Layer Boundary Simulation: The 2-D isotropic velocity structure in both routines is parameterized into layers composed of variably sized trapezoids. Each layer boundary is specified by an arbitrary number of boundary nodes, with arbitrary spacing, connected by straight Une segments. A smooth layer boundary simulation (Zelt and Smith 1991) is used to reduce the effects of ray scattering and focusing which can arise with this type of parameterization, particularly for rays incidents on a boundary near a node separating two boundary segments with different slopes. The smooth layer boundary simulation affects only the incident and emergent ray angles where the slope of the smoothed boundary is used. Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 34 (c) Node Number and Locations: As mentioned before, the number and locations of nodes were determined by several inversion tests so that an optimum trade-off between travel time fit and model parameter resolution was achieved while still allowing rays to be traced to all observations. In general, it is desirable to look for the simplest (i.e. least structure) model that replicates the observations; thus the one with the fewest parameters was chosen. The number of independent model parameters can be reduced by fixing boundary and velocity nodes for some iterations, specifying a zero velocity discontinuity across a boundary and fixing the vertical velocity gradient within a layer. These methods are generally used to increase the stability of the inversion when the travel time data do not provide strong constraints on a particular feature of the velocity structure but one which may be inferred from either the corresponding reflection data (e.g., geometry of an upper crustal boundary) or the amplitude modelling which is sensitive to velocity gradient. (d) The "Layer Stripping" Approach: A "layer stripping" approach was followed here to model the observed data. Typically, travel time data were inverted to solve for boundary depths and velocities simultaneously within each layer while the values of the parameters of the layers above were fixed. The relative amplitudes of phases used in the inversion and previously modelled phases were compared qualitatively with observed amplitudes, and velocity gradients and/or velocity contrasts across reflecting boundaries were adjusted by forward modelling to improve the match. In some cases, velocities were fixed as required by the amplitude data. In other cases, small adjustments to the parameters determined by travel time inversion were required to remove artificial shadow zones or to correct physically unrealistic situations. It was necessary to make sure that all these adjustments were within the estimated absolute parameter uncertainties. This procedure of alternating travel time inversion and amplitude forward modelling was repeated until a satisfactory fit of observed and calculated travel times and amplitudes Chapter 3. DATA ANALYSIS AND INTERPRETATION PROCEDURE 35 was achieved for successively deeper boundaries. This method clearly simplifies and expedites the inversion process. But it is also essential to note: (1) the errors in parameters of upper layers cannot be further resolved by inversion of deeper layers; (2) this method is better suited for a set of thin layers which is not always the case in the models with which we are dealing. So in some cases, it is necessary to go back and check the sensitivities of model parameters in the upper layers. (e) Damping and Uncertainty of Model Parameters: Overall damping and a priori estimates of the uncertainty of model parameters are important in a damped least squares inversion. The damping parameter determines the overall trade-off between the resolution and uncertainty of model parameters (Zelt and Smith 1991). Unfortunately, it can only be decided by trial-and-error. Various levels of damping were tried, but for the final model a damping factor of 1.0 (i.e., equal weight placed on model resolution and uncertainty) was used. The estimated uncertainties of model parameters were 0.1 km/s and 1.0 km for velocities and boundary depths, respectively. (f) Other Considerations: In practice, a trade-off between reducing the travel time residual and improving the computed amplitudes is inevitable. More emphasis has been placed on achieving a satisfactory travel time fit, because of the limitations involved in calculating amplitudes using ray theory. In applying the modelling techniques following the procedure described above, one basic assumption is made: the inversion algorithm is capable of resolving trade-offs be-tween velocity and boundary depth. Zelt and Smith (1991) asserted that testing has revealed this assumption to be valid, though in practice, velocities are generally better determined than boundary positions. Chapter 4 M O D E L I N G A N D I N T E R P R E T A T I O N OF R E F R A C T I O N D A T A 4.1 Modelling the Refraction Data 4.1.1 Initial Model A brief description of the basis for initial model construction is given here. As mentioned in the interpretation of reflection data, upper sedimentary boundaries and sediment veloc-ities obtained from reflection data were used in the initial refraction model, and actually fixed throughout the modeling procedure. Because deep velocity structure is the focus of this study, the effect of local sedimentary velocity changes on the refraction arrivals are not emphasized. The boundaries of Crescent terrane and Pacific Rim terrane were introduced on the basis of previous studies; i.e., the interpretation of 1985 deep seismic reflection data (Clowes et al. 1987b), and then integrated with extensive modeling of potential field data (Dehler and Clowes 1992, Clowes et al. 1996) The emphasis of this study has been put on the velocities rather than the boundary positions. Because of this and the fact that the refraction data provide little constraints on the boundaries of Crescent and Pacific Rim terranes, they were fixed during the modeling procedure. The positions of the subducting oceanic plate, interpreted from the reflection data, were used in the initial model. But they were not fixed during the modeling, since the transform from two-way traveltime to depth is dependent on the velocity structure. Also the boundary positions are constrained by the wide-angle reflection data. 36 Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 37 4.1.2 Line 89-02/03 This section describes the application of the modelling procedure to the Line 89-02/03 data set in order to construct the final velocity structure model (Figure 4.1). It points out the degree and source of constraints supporting each model feature, and the fit of traveltime and amplitude of each phase. Figures 4.2a and 4.4a show the observed data section recorded at site 2A and 2B, respectively, while Figures 4.2b and 4.4b show synthetic seismograms determined from the final velocity model for comparison with the observed data. Figures 4.3a and 4.5a show exemplary rays traced from receiver sites 2A and 2B, respectively, with Figures 4.3b and 4.5b showing the match between observed and calculated travel times. The locations of depth and velocity nodes are shown in Figure 4.11a. Figure 4.12a shows the total ray coverage of line 89-02/03. Note that possible sources of error in the 2-D interpretation are out-of-plane effects due to 3-D variation in velocity and reflector surfaces, particularly since the recording of 89-02 is not strictly inline with 89-03. Although these effects are expected to be minimal, this assumption cannot be assured in any meaningful way. The final velocity model is divided into three regions: the continental crust, the oceanic crust and the oceanic upper mantle. (a) The Continental Crust Only part of the continental crust is constrained by phases A and C. The features of interest to which the data apply include: (1) Pacific Rim Terrane, the upper portion of which is constrained by A; (2) Crescent Terrane, the upper part of which is constrained by A and the lower part constrained by C; (3) the landward end of Accreted Wedge, the upper part of which is constrained by A and the lower part constrained by C; and (4) the weakly constrained upper and western corner of Wrangellia. 38 2 A 2 B L I N E 8 9 - 0 3 | | L I N E 8 9 - 0 2 | # * -140 -120 -100 -80 -60 -40 -20 0 20 i i i i i i I I i 1.48 3.00 4.00 5.14 5.71 6.29 7.00 8.00 8.40 Velocity (km/s) Figure 4.1: (a) Color contour plot of velocities beneath line 89-02/03 (vertical exaggeration approximately 4:1). Numbers within model represent P-wave velocities in km/s at various locations. Heavy solid lines represent locations of reflectors constrained by the traveltime data. * represents onshore receiver site. The shaded region in the model is unconstrained by the traveltime data, (b) 1:1 version of model in (a). SED, sedimentary layer; AW, accreted wedge; CT, Crescent terrane; PRT, Pacific Rim terrane; WR, Wrangellia; JdF, Juan de Fuca plate. D I S T A N C E ( k m ) - 1 4 0 - 1 2 0 - 1 0 0 - 8 0 - 6 0 - 4 0 - 2 0 0 2 0 D I S T A N C E ( k m ) Figure 4.3: Two-point ray tracing diagram showing ray coverage (a) and comparison between observed (short vertical bars representing uncertainty of picks) and calculated (thin lines) traveltimes (b) for line 89-02/03 at station 2A. Distance (km) Figure 4.4: (a) Observed record of line 89-02/03 at Station 2B. (b) Synthetic section for line 89-02/03 at Station 2B. DISTANCE ( k m ) —i 1 1 1 1 1 i i r - 1 4 0 - 1 2 0 - 1 0 0 - 8 0 - 6 0 - 4 0 - 2 0 0 2 0 DISTANCE ( k m ) Figure 4.5: Two-point ray tracing diagram showing ray coverage (a) and comparison between observed (short vertical bars representing uncertainty of picks) and calculated (thin lines) traveltimes (b) for line 89-02/03 at station 2B. Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 43 The boundary of the upper sedimentary layer in the uppermost crust was fixed be-cause it is mainly constrained by the corresponding reflection data. The velocities of the sedimentary layer are also fixed using the values derived from the velocity analysis of reflection data. The boundary locations of Pacific Rim Terrane, Crescent Terrane and Accreted Wedge were also fixed, because they cannot be resolved effectively by the data. The resulting model has an acceptable traveltime fit for phases A and C. The amp-litude match is overall good except some local mismatch. This usually indicates the requirements of introducing more layers, but one of our modeling objective is to make the model simple based on known geometry, which also helps to maintain the integrity of the inversion algorithm. The high amplitude of A indicates a strong vertical velocity gradient in the upper crust (i.e. upper portions of Pacific Rim Terrane, Crescent Terrane and seaward end of Accreted Wedge). This is consistent with the results of velocity ana-lysis of the reflection data (Spence et al., 1991). There is no constraint on the velocity gradient in the lower crust, which is therefore fixed at approximately 0.1 s - 1 for the inversion procedure, based on a priori estimates of the velocities at the top and at the base of the lower crust. (b) The Oceanic Crust Only part of the oceanic crust was constrained by the data set. A small portion of the top of the oceanic plate was constrained by reflection event C. The position of the top of the oceanic plate was resolved by modelling of refraction data constrained by the interpretation of corresponding reflection data. The velocity of the part of the oceanic crust which is below the seaward Accreted Wedge is constrained by M l , which also provides a constraint for the position of oceanic Moho at the base of the area. The M l data do not constrain a vertical velocity gradient, and since there was no evidence for strong gradients, the gradient in the oceanic crust Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 44 was set to zero. But a two-layer structure was imposed to be consistent with the general understanding of the oceanic crust (Clowes and White, 1990). The traveltime fit of phase M l is acceptable on both sites 2A and 2B. But the ampli-tude fit on site 2A becomes difficult to monitor because of the heavy ringing in the data section. On the other hand, the amplitude fit of site 2B is reasonably good. (c) The Oceanic Upper Mantle Both the velocity and velocity gradient of the oceanic upper mantle were constrained by refraction event M2. But the coverage of the data was limited due to the gap between lines 89-02 and 89-03, and the low signal-to-noise ratio in the central part of the phase. The traveltime fit of phase M2 is overall acceptable for both sites 2A and 2B although some local misfits are probably caused by abrupt velocity changes in the highly deformed upper continental crust. The amplitude fit is reasonable for phase M2 on both sites 2A and 2B. 4.1.3 Line 89-09 This section describes the application of the modelling procedure to the Line 89-09 data set in order to construct the final model (Figure 4.6). It point out the degree and source of constraints supporting each model feature, and the fit of traveltime and amplitude of each phase. Figures 4.7a and 4.9a show the observed data section recorded at site 4A and 4B, respectively, while Figures 4.7b and 4.9b show synthetic seismograms determined from the final velocity model for comparison with the observed data. Figures 4.8a and 4.10a show exemplary rays traced from receiver sites 4A and 4B, respectively, with Figures 4.8b and 4.10b showing the match between observed and calculated travel times. The locations of depth and velocity nodes are shown in Figure 4.11b. Figure 4.12b shows the total ray coverage of line 89-09. 45 4A 4B LINE 89-09 | % % Distance (km) _._ : — ] ~ : . i I t I I I i I I 1.48 3.00 4.00 5.14 5.71 6.29 7.00 8.00 8.40 Velocity (km/s) Figure 4.6: (a) Color contour plot of velocities beneath line 89-09 (vertical exaggeration approximately 4:1). Numbers within model represent P-wave velocities in km/s at various locations. Heavy solid lines represent locations of reflectors constrained by the traveltime data. * represents onshore receiver site. The shaded region in the model is unconstrained by the traveltime data, (b) 1:1 version of model in (a). SED, sedimentuiy layer; AW, accreted wedge; CT, Crescent terrane; PRT, Pacific Rim terrane; WR, Wrangellia; JdF, Juan de Fuca plate. 80 \1R D I S T A N C E ( k m ) - 1 2 0 - 1 0 0 - 8 0 - 6 0 - 4 0 - 2 0 0 20 D I S T A N C E ( k m ) Figure 4.8: Two-point ray tracing diagram showing ray coverage (a) and comparison between observed (short vertical bars representing uncertainty of picks) and calculated (thin lines) traveltimes (b) for line 89-09 at station 4A. Figure 4.10: Two-point ray tracing diagram showing ray coverage (a) and comparison between observed (short vertical bars representing uncertainty of picks) and calculated (thin lines) traveltimes (b) for line 89-09 at station 4B. Figure 4.11: (a) Locations of depth (solid squares) and velocity (solid circles) nodes used to define the line 89-02/03 model. Empty circles are fixed velocity nodes, (b) Node locations of line 89-09. Vertical exaggeration is approximately 4:1. C T , Crescent Terrane; PRT, Pacific Rim Terrane. 51 Figure 4.12: (a) Two-point ray tracing diagrams showing ray coverage of line 89-02/03. (b) Ray coverage of line 89-09. Vertical exaggeration is approximately 2:1 in both (a) and (b). Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 52 Similar to Line 89-02/03, the final velocity model of Line 89-09 is divided into three regions: the continental crust, the oceanic crust and the oceanic upper mantle. (a) The Continental Crust Only part of the continental crust is constrained by phases A and C of Line 89-09. The features of interest to which the data apply include: (1) Pacific Rim Terrane, the upper portion of which is constrained by A; (2) Crescent Terrane, the upper part of which is only weakly constrained by A due to a significantly decreased thickness; (3) the landward end of Accreted Wedge, which is constrained by A and C; and (4) the weakly constrained upper and western corner of Wrangellia. Both the boundary and velocity of the sedimentary layer were fixed using the results derived from reflection data. The boundary locations of Pacific Rim Terrane, Crescent Terrane and Accreted Wedge were also fixed due to lack of constraints from the data. The resulting model has an overall acceptable traveltime fit for phases A and C. But some mismatch exists on phase A due to the apparent "fade out" of data on phase A from 25-35 km distance on both sites 4A and 4B. This might be caused by the strong attenuation effect of the highly deformed upper crust sediments since a similar "white out" can be seen in the same data section on the corresponding reflection data of line 89-09 (Figure 3.2). This may have also affected the amplitude of both events A and C; therefore, the obvious mismatch on part of both phases is difficult to correct with the theoretical approach used. The amplitude information indicates that the gradient in the upper crust is not as high as Line 89-02/03. This is not confirmed due to the "fade out" of data in phase A and lack of detailed velocity analysis of the reflection data. Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 53 (b) The Oceanic Crust Compared with line 89-02/03, a smaller portion of the top of the oceanic plate was con-strained by reflection event C. Similarly, the position of the top of the oceanic plate was resolved by modeling of refraction data constrained by the interpretation of corresponding reflection data. The traveltime fit of phase Ml is acceptable on both 4A and 4B. But the amplitude fit becomes very difficult to accomplish due to the heavy attenuation effect in the data section as mentioned above. (c) The Oceanic Upper Mantle Similar to Line 89-02/03, both the velocity and its gradient of the oceanic upper mantle were constrained by refraction event M2. The traveltime fit of phase M2 is overall acceptable for both sites 4A and 4B although some local misfits are probably caused by abrupt velocity changes in the highly deformed upper continental crust. The amplitude fit is reasonable for phase M2 on both sites 4A and 4B. 4.2 Primary Features of the Final Models A brief description of the primary features of the final models for Line 89-02/03 and Line 89-09 is given below. Further discussion and comparison with other data sets will be presented in Chapter 5. 4.2.1 Line 89-02/03 The principal features of the model are shown in Figure 4.1. Below the water layer, the uppermost sedimentary stratum is a thin (<4 km) layer with a trend of velocity increase from west to east. Beneath the continental slope, the velocities average 3.6 km/s at the Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 54 top and 4.0 km/s at the bottom, with a contrast to the uppermost layer of continental crust where the velocities average 5.3 km/s at the top and 6.2 km/s at the bottom. The structure of this sedimentary layer is mainly constrained by the reflection data. The velocities of the crust below the western part of the continental shelf are assumed homogeneous with 4.0 km/s in the upper portion and 4.5 km/s in the lower portion, since this part of the model is not well constrained by the data. The model structure beneath the eastern continental shelf is complex compared to that of the west. The crustal thickness increases landward from 15 to 17-20 km. The detailed boundaries of the Accreted Wedge, Crescent and Pacific Rim Terranes are in-serted according to previously known information. At the top of the crust, the velocity increases eastward from 4.3 km/s in the Accreted Wedge Terrane to 4.9 km/s in the Pacific Rim Terrane. At the bottom of the crust, there is a high velocity zone right at the bottom of the Crescent Terrane with a velocity of 6.4 km/s. The model shows a velocity of 5.5 km/s for the lower portion of the landward end of the Accreted Wedge. Due to poor constraints on the seaward portion of the Accreted Wedge provided by the data, the rapid velocity change from 4.5 km/s to 5.5 km/s may not be real. But the velocity of 5.5 km/s is reasonable for the landward portion of the Accreted Wedge which contains very consolidated sediments lying below and against the Crescent terrane. With a velocity of about 4.5 km/s at its top, the Crescent terrane has a very high velocity of 6.4 km/s at its bottom part which is well constrained by the data. This may suggest that this part of Crescent terrane was highly compressed during and after the accretion since it has been serving as a backstop of the Accreted Wedge sediments. The alternative suggestion is that the Crescent terrane may have been truncated at certain time of its history, which was previously suggested by Clowes et al. (1987a). In other words, the actual Crescent terrane might be smaller than it appears in this final model. The Pacific Rim terrane has a higher velocity of 4.9 km/s at the top than the Crescent Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 55 terrane. The bottom part appears to have a velocity around 6.2 km/s although this part of the model is not well constrained by the data. The velocities of the continental crust are assumed homogeneous with 6.2 km/s at the top and 6.4 km/s at the bottom, since this part of the model is only poorly constrained by the data set. The position of the top of the subducting oceanic plate is determined by examination of both reflection and refraction data. The upper crust of the oceanic plate has velocities which increase from 5.5 km/s in the west to 6.5 km/s in the east. The lower crust also has a trend of velocity increase which is 6.5 km/s in the west and 7.5 km/s in the east. No velocity gradient is assigned since this part of the model is only partly constrained by the data. The lateral velocity change in the oceanic crust is interesting since this may be related to the earthquake source zone indicated by Hyndman and Wang (1995). Further discussion on this topic will bb presented in Chapter 5. The oceanic Moho was determined from the M l and M2 phases. We found a normal oceanic upper mantle velocity around 8.0 km/s which has an unusual lateral variation and a vertical gradient around 0.01 s - 1 . But this part of the model has very low resolution due to the unstable inversion results. So these variations may not be real. 4.2.2 Line 89-09 The principal features of the model are shown in Figure 4.6. The uppermost sedimentary layer is thinner than that of line 89-02/03, but the velocity also increases from 3.0 km/s in the west to 4.0 km/s in the east below the continental slope. The structure is determined by the reflection data. Similar to line 89-02/03, the velocities of the crust beneath the western part of the continental slope are assumed to be 4.0 km/s at the top and 4.5 km/s at the bottom. The structure of the crust beneath the eastern continental shelf is more complex. The Crescent Terrane becomes much smaller than that of line 89-02/03 if it exists at all, which cannot Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 56 be confirmed in the corresponding reflection data. Crustal velocities increase landward, from 4.2 km/s to 5.0 km/s at the top, and from 5.5 km/s to 6.2 km/s at the bottom. A high velocity zone is found at the east lower corner of the Accreted Wedge with a velocity of 5.5 km/s, which is reasonable for this very consolidated part of Accreted Wedge sediments. The Pacific Rim terrane has a velocity of 5.0 km/s at the top, and 6.0 km/s at the bottom which is not well constrained by the data. The velocities of the continental crust are assumed homogeneous with 6.0 km/s at the upper portion and 6.2 km/s at the lower portion, since this part of the model is only poorly constrained by the data set. The top of the oceanic plate is about 2 km shallower than that of line 89-02/03 below the oceanic basin and 3-4 km shallower under the continental slope. The velocity of the upper crust increases from 6.0 km/s in the west to about 6.5 km/s in the east. It is interesting that there is no big lateral velocity variation in the lower crust as for Line 89-02/03. This may indicate a difference in the nature of earthquake source zone between Line 89-09 and Line 89-02/03 which will be discussed further in Chapter 5. The position of oceanic Moho was determined from the M l and M2 phases. We found a normal oceanic upper mantle velocity around 8.0 km/s which has an unusual lateral variation and vertical gradient around 0.01 s - 1 . Similar to Line 89-02/03, the resolution is low in this part of the model and these variations may not be real. 4 . 3 Non-uniqueness of the Model Non-uniqueness is inherent to the final model because an acceptable fit to the data can be produced by two different models which are both well-resolved and physically acceptable. The nature of this non-uniqueness problem is outlined below. The non-uniqueness of the solution is evident from the subjective choice of the node Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 57 placement in the modelling procedure. Furthermore, for a certain choice of node place-ment, different a priori estimates of parameter uncertainties, different damping para-meters and various starting models are associated with a different final model. Even within a certain parameterization, there is a sense of non-uniqueness with any nodal value output from the R A Y I N V R inversion since each calculation has an associated standard error, indicating a range of feasible values. Therefore, it must be stressed here that the final model presented above is a model which reasonably fits the data, not a unique solution of the inverse problem. 4.4 Model Resolution and Parameter Uncertainty Zelt and Smith (1991) described tests which provide estimates of (1) the spatial resolution of a model about a specific velocity or boundary node, and (2) the absolute parameter uncertainty of a particular node. 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 diag[o"i] - the perturbation should be large enough to change travel time, but small enough that the ray path pattern is not disturbed significantly. The next step is to trace rays through the perturbed model to calculate the travel time 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 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 of 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 the other parameter values will resemble the non-perturbed values; poor resolution of the selected Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA 58 parameter will result in the perturbation being "smeared" into the adjacent parameters. Another test was suggested by Zelt and Smith (1991) to estimate the absolute un-certainty of a model parameter. The first step is again to choose a model parameter to perturb, then to invert all other parameters which 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 con-tinually 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-square values of the two final models. The statistical F-test involves a ratio calculated such that values either >>1 or <<1 indicate very signific-ant 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 obtained for the final models are shown in Table 4.1. Normally, the resolution decreases and the uncertainty increases with depth in the model. This testing procedure also revealed a trade-off between resolution and uncertainty based on node spacing; compared with widely spaced nodes, the more closely spaced nodes were shown to have a decrease in resolution but an improvement in absolute uncertainty. Chapter 4. MODELING AND INTERPRETATION OF REFRACTION DATA Estimated Estimated lateral absolute resolution uncertainty Nodes (km) Velocities at top of AW, CT and PRT 20-25 0.2-0.3 km/s Velocities at bottom of AW, CT and PRT 10-20 0.2-0.3 km/s Depth nodes at top of oceanic plate 25 1.5 km Velocities of upper oceanic crust 25-35 0.5 km/s Velocities of lower oceanic crust 25-35 0.5 km/s Depth nodes at oceanic Moho 50 2.0 km Velocities of oceanic mantle >50 0.5 km Table 4.1: Estimated lateral resolution and absolute uncertainty of velocity and depth nodes in the final models. Chapter 5 DISCUSSION A N D C O N C L U S I O N Any model of the Earth should be consistent with all available data, and a model derived from seismic data must be consistent with the observed gravity field. Based on this we convert our final velocity models to density and generate corresponding gravity profiles to compare with the observed gravity data. We also compare our density model with the integrated block models of density structure developed by Clowes et al. (1996). Then we discuss our results with respect to the theory of great thrust earthquakes in the Cascadia subduction zone as presented by Hyndman and Wang (1993, 1995) from thermal and deformation data. A summary of what we have achieved in this study will be given at the end of this chapter. 5.1 Testing of Consistency with Gravity Data and Comparison with Density Structures Clowes et al. (1996) present a series of interpreted block models of density structure developed through the integration of numerous geophysical and geological data. The models illustrate variations in extent, structure and properties of the major tectonic components along the margin, including the Crescent terrane, Pacific Rim terrane, the top of the Juan de Fuca plate, etc., and provide a good framework within which to examine recent seismic and other data. The Bouguer gravity anomaly map for the study area and locations of gravity lines IAS and 2A of Dehler and Clowes (1992) are shown in Figure 5.1. 60 -210-185-160-135-110-85 -60 -35 -10 15 40 65 90 Bouguer/Free-Air Gravity (mgal) Figure 5.1: The Bouguer gravity anomaly map for the study area. The locations of gravity lines IAS and 2A, .seismic lines 89-02/03, 89-09 and corresponding onshore stations are indicated. Scale is approximately same as Figure 1.2. Chapter 5. DISCUSSION AND CONCLUSION 62 As the first step in testing the consistency of our velocity models with density models derived from gravity data, we have to convert the velocity structural models to density models with a constant density value for each polygonal block. This requires the conver-sion of the velocity model with velocity gradients in most blocks to a model with constant velocity in each block. To keep the conversion simple and reasonably accurate, we calcu-lated the velocity at the midpoint (midpoint from the top to bottom of the layer) of each velocity node using the known velocity gradient, then used the average of all "midpoint" velocities within one block as the constant velocity for each block. Then we converted the resulting velocity model to a density model using the traditional velocity-density relation-ship (Ludwig, Nafe and Drake 1970) as shown in Figure 5.2. A forward modeling routine (Webring, 1985) was used iteratively to calculate the associated gravity anomalies and compare these with observed values. The minimum-maximum limits of Barton (1986) were used to impose bounds on density variations during the modeling procedure. We did not attempt to match every detail of the observed gravity data because the purpose of our calculation was to test consistency with the observed gravity field. So we kept the number of iterations to a minimum (less than 3). The resulting fit with the observed data is shown in Figure 5.3a and 5.4a for line 89-02/03 and 89-09, respectively; the resulting density model is shown in Figure 5.3c and 5.4c. The reasonable fit shows clearly that our final velocity models are generally consistent with the observed gravity anomaly. But the ease with which a model fit could be achieved indicates the possible low resolution and non-uniqueness of this converted model. This aspect has been discussed by Barton (1986). 5.1.1 Comparison between Line 89-02/03 and Line IAS The continuous and extensive coverage of gravity and magnetic data (Dehler and Clowes, 1992 and Clowes et al, 1996) provides an opportunity to extend interpretations based Chapter 5. DISCUSSION AND CONCLUSION 63 T o o o o J ! ! ! ! I 2 3 4 5 Density , gem *.3 Figure 5.2: Laboratory measurements of P-wave seismic velocity and density in rocks, after Ludwig, Nafe and Drake (1970). The thin line in the center of the scatter is the mean velocity-density relationship commonly used for gravity calculations based on seismic models. The heavy lines show the maximum and minimum density values for each seismic velocity. Open circles represent measurements of metamorphic and igneous rocks; closed circles represent sediments. Chapter 5. DISCUSSION AND CONCLUSION 64 on single-point and profile data across and along the margin. The gravity anomaly field (Figure 5.1) is reasonably two-dimensional along the margin, allowing for variability in the geometry of the sedimentary wedge, thereby justifying the 2.5-D modeling approach. The Crescent terrane in density model matches reasonably well with the velocity model, although the detailed velocity model (Figure 4.1) shows that higher velocity/density rocks maybe located in the lower part of the terrane. With strong positive magnetic anomalies clearly associated with the Eocene volcanics of the Crescent terrane, an adjacent negative anomaly appears to characterize the meta-morphosed and deformed sedimentary rocks of the Pacific Rim terrane (Clowes et al. 1996), whose existence is indicated in both density and velocity models. But the model values from velocity modeling and gravity modeling are not perfectly matched. We be-lieve it is caused by the low resolution in modeling procedure, especially in the gravity modeling. The location of the top of the subducting Juan de Fuca plate in each model matches reasonably well down to about 18 km. The density values of oceanic crust and upper mantle in both models also agree with each other, although the lateral variations can not be well resolved from either data set due to the data quality and fundamental ambiguities in modelling procedure. Overall the values of density converted from velocity model are lower than densities for the model derived from gravity data. It can also be noted in Figure 5.3a that the gravity anomalies calculated from the densities of the converted model are consistently less than the observed values. These are probably results of the limited reliability of the relationship between velocity and density. But in general, the achieved velocity model reasonably fits the observed gravity field. Chapter 5. DISCUSSION AND CONCLUSION 65 Figure 5.3: Gravity modeling of Line IAS and 89-02/03. (a) Calculated gravity anomalies for converted density model of Line 89-02/03 (shaded line) and density model of Line IAS (solid line), compared with observed data ("+" dotted line), (b) Density model of line IAS (Dehler and Clowes 1992). (c) Converted density model of line 89-02/03 using the mean Nafe-Drake velocity-density conversion. In (b) and (c): CT, Crescent terrane; Jdf, Juan de Fuca Plate; PRT, Pacific Rim ter-rane. Heavy bars in (c) were boundaries added to the converted density model. Vertical exaggeration is 2:1. Chapter 5. DISCUSSION AND CONCLUSION 66 5.1.2 Comparison between Line 89-09 and Line 2A The density model of Line 2A (Figure 5.4b) was constructed by using the step-wise modeling method (Dehler and Clowes 1992) to provide structural estimates for a margin segment poorly-constrained by the data. The accumulation of accretionary material lying beneath the continental slope be-comes much thinner, which is similar to what the velocity model has indicated. Also the gravity profile is characterised by a much broader and lower amplitude negative anomaly than Line IAS. Within the wedge, the Crescent terrane still appears to be present in the density model of Line 2A. The velocity model of Line^ 89-09 indicates a much less prominent Crescent terrane than Line 89-02/03, although it does appears to exist according to the corresponding reflection image. The Crescent terrane is absent on the lines further north (Spence et al. 1991). On the density model, the adjacent negative anomaly supports the continued presence of the Pacific Rim terrane, but again at a slightly reduced width, which can also be seen in the velocity model. The top of the subducting Juan de Fuca plate becomes shallower in Line 2A, which agrees very well with the velocity model of Line 89-09. Also, the density value for oceanic crust and upper mantle in both models agree with each other very well. Similar to line 89-02/03, the gravity anomalies calculated from the converted model are overall less than the observed values caused by the limited reliability of the rela-tionship between density and velocity. But in general, the achieved velocity model is consistent with the observed gravity constraints. Chapter 5. DISCUSSION AND CONCLUSION 67 Figure 5.4: Gravity modeling of Line 2A and 89-09. (a) Calculated gravity anomalies for converted density model of Line 89-09 (shaded line) and density model of Line 2A (solid line), compared with observed data ("+" dotted Une). (b) Density model of Une 2A (Dehler and Clowes 1992). (c) Converted density model of Une 89-09 using the mean Nafe-Drake velocity-density conversion. In (b) and (c): CT, Crescent terrane; Jdf, Juan de Fuca Plate; PRT, Pacific Rim terrane. Heavy bars in (c) were boundaries added to the converted density model. Vertical exaggeration is 2:1. Chapter 5. DISCUSSION AND CONCLUSION 68 5.2 Does the Velocity Model Relate to Potential Great Thrust Earthquakes? The Cascadia subduction zone off Vancouver Island is unusual in having experienced no historical great thrust earthquakes. But most similar margins around the world have had such very damaging events. The comprehensive geophysical data base and detailed struc-tural and tectonic information for the southern Vancouver Island Lithoprobe corridor has proven invaluable in studies of the subduction zone great earthquake hazard (Hyndman, 1995). The region now has the most complete data documenting the nature and extent of great earthquakes (M > 8) on the Cascadia margin, including paleoseismicity, geodetic data and thermal data (e.g., Hyndman, 1995). Extensive paleoseismicity data indicate great earthquakes have occurred in the past at irregular intervals averaging about 600 years, the last having occurred about 300 years ago (e.g., Clague and Bobrowsky, 1994, 1995). Strong shaking of the continental shelf and slope region is inferred from turbidite sediment layers, found in Cascadia Basin deep sea channels, produced by mobilization of sediment on the continental slope and shelf (Adams 1990). Also, there is a record of abrupt vertical downdrop of coastal regions in buried intertidal marsh surfaces for Vancouver Island (Clague and Bobrowsky, 1994, 1995). The average interval is roughly 500 years, in general agreement with the average interval between the deep sea turbidite layers. The most recent event has been consistently radiocarbon dated at approximately 300 years (e.g., Atwater et al., 1991). Geodetic data have shown that the region is short-ening in the direction of plate convergence, as expected for a locked subduction thrust fault (e.g., Hyndman and Wang, 1995). Other detections of present flexural uplift and shortening across the margin (e.g., Savage et al., 1991; Dragert et al., 1994) also provide strong evidence that the Cascadia subduction thrust fault is indeed locked and probably will fail in future great earthquakes. It becomes very important to estimate the downdip and updip limit of the locked zone Chapters. DISCUSSION AND CONCLUSION 69 0 50 100 150 Distance from trench axis (km) Figure 5.5: Location map and model profile of thermal analysis, (a) Location map with the widths of the locked and transition zones from thermal ana-lysis (thick short lines) compared to those from current deformation (stippled area). The bands on the thermal results (350 and 450 degree) indicate the estimated uncertainties (after Hyndman and Wang. 1995). (b) Model temperatures across the northern Cascadia margin (Southern Vancouver Island profile in (a)) showing the positions on the subduc-tion thrust fault of the 350 and 450 degree downdip limits of the locked and transition zones. Chapter 5. DISCUSSION AND CONCLUSION 70 along the subduction thrust fault since such information would lead to a much better understanding of the subduction thrust system and further indication of the magnitude and resulting damage of the future great earthquake. Hyndman and Wang (1995) have estimated the downdip landward limit of the seismic rupture zone (Figure 5.5) on the subduction thrust fault for the whole Cascadia margin from (1) the locked zone from dislocation modeling of current deformation data, and (2) the thermal regime, taking the downdip limit of seismic behavior on the fault to be controlled by temperature. The resulting downdip width of the seismic rupture zone is about 60 km locked and 60 km transition. However they did not infer a convincing updip seaward limit of the rupture zone since it is poorly constrained by the data, although they mentioned that the updip limit could be at the deformation front. This updip limit could become very important if it indicates the updip aseismic zone is sufficiently wide such as to significantly reduce the rupture zone width. Their conclusion suggesting that the landward limit to the seismic rupture zone extended little, if at all, beneath the coast of Vancouver Island, limits the ground motion from great subduction earthquakes at the larger cities, such as Victoria, Vancouver and Seattle, that lie 100-200 km in land. The narrow width also limits the earthquake size but events of magnitude well over 8 are still possible; the maximum depends on the along-margin length. The location of our final model of Line 89-02/03 (Figure 4.1) is very close to what is shown in Figure 5.5. Can the final velocity model provide any further information or constraints with regard to this seismic rupture zone? The nature of the crust of subducting oceanic plate - the extent to which it is consolidated and metamorphosed -will have a great deal to do with the frictional behavior of the subduction thrust fault interface (Scholtz 1990). And the velocity information of the subducting oceanic crust is a good indication of its degree of consolidation. Scholtz (1990) also suggested that the Chapter 5. DISCUSSION AND CONCLUSION 71 critical transition of the seismic rupture zone is actually that from seismogenic velocity-weakening to stable-sliding velocity-strengthening behavior. This indicates that the part of subducting oceanic crust with unconsolidated, weak material (i.e. with lower velocity) tends to be in a stable transition and velocity-strengthening situation with less stress-buildup. This means that this part of the oceanic crust is not likely being locked. On the other hand, the part of oceanic crust with very consolidated material (i.e. with higher velocity) is in an unstable and velocity-weakening situation with high stress-buildup. It is more likely that this part of the oceanic crust could be locked and would eventually fail in the great thrust earthquake. But it has to be emphasized that higher velocity itself does not point to the conclusion that the zone is locked, which also relies on other factors such as thermal and deformation information. More detailed discussion can be found in Scholtz (1990). The above discussion suggest that we may be able to estimate the boundary of the seismic rupture zone with the help of velocity information. In the case of the velocity model of Line 89-02/03 (Figure 4.1), there is an apparent west-to-east lateral velocity change from 5.5-6.5 km/s to 6.5-7.5 km/s in the subducting oceanic crust, as mentioned in the last chapter. This appears to be a transition from a velocity-strengthening zone with unconsolidated weak material (lower velocity) to a velocity-weakening zone with consolidated material (higher velocity). Thus the lateral velocity change may provide an indication concerning the updip Hmit of the seismic rupture zone in this area. Assuming the preceding approach is valid, the final model in Figure 4.1 indicates a possible updip (seaward) boundary for the seismic rupture zone at about -60 km distance. Comparing this position with the seismic rupture zone in Figure 5.5, this possible updip Hmit Hes within the "transition zone". Considering resolution and uncertainty of our velocity model, the least indication we can make is that the width of the total seismic rupture zone (locked zone -f- transition zone) would be significantly reduced. This is Chapter 5. DISCUSSION AND CONCLUSION 72 very important since it could further indicate that the maximum earthquake size (the actual size also depends on the along-margin length) would be substantially limited if our hypothesis is valid (Hyndman 1995). It is interesting that in the final velocity model of Line 89-09 (Figure 4.6), we see a lateral velocity change only in the upper oceanic crust but not in the lower oceanic crust. This may indicate that the width of the seismic rupture zone changes along the margin. But it is difficult to make such an inference without the opportunity to compare results with a corresponding profile generated from thermal and geodetic data. 5.3 Conclusions Through the application of successive iterations of inverse travel time modelling and forward amplitude modelling to the refraction data, with the aid of interpretation of corresponding reflection data, velocity structures were developed for Lines 89-02/03 and 89-09 in the Cascadia subduction zone (Figure 4.1 and 4.6). From these two models, a number of important results concerning the geotectonic structure of the area have been derived: (1) The velocity models for both lines enable clear identification of a narrow strip of allochthonous Pacific Rim terrane. This indicates its probable existence along the southern portion of the west coast of Vancouver Island, which is consistent with the recent density model derived from gravity anomaly data (Dehler and Clowes, 1992; Clowes et al., 1996). Such a geometry is also consistent with emplacement of the terrane by strike-slip faulting or transpression. (2) The velocity structure along Une 89-02/03 has allowed the volcanic Crescent ter-rane, probably Eocene oceanic crust formed in a marginal basin (Massey, 1986), to be clearly identified adjacent to Pacific Rim terrane. The velocity model also suggests higher Chapter 5. DISCUSSION AND CONCLUSION 73 velocity for the lower half of Crescent terrane, with a high velocity of 6.4 km/s at the bottom of the terrane which is almost the same velocity as the underlying upper oceanic crust. This may suggest that this part of the Crescent terrane was highly compressed during and after the accretion since it has served as a backstop of the Accreted Wedge sediments. This is also consistent with its possible origin as a piece of oceanic crust un-derthrust by Juan de Fuca plate (e.g. Massey, 1986; Hyndman et al., 1990). A possible alternative suggestion is that the Crescent terrane may have been truncated at certain time of its history which was previously discussed by Clowes et al (1987a). (3) The velocity model along line 89-09 suggests that the Crescent terrane might fade away toward north part of the margin, which is also indicated by the magnetic and gravity anomaly data (Dehler and Clowes, 1992). The corresponding reflection data also suggests a much weaker possible signature for the Crescent terrane. (4) The velocity structure of the Accreted Wedge is very similar in both lines, with a high velocity zone of 5.5 km/s right beneath the adjacent Crescent terrane but lower velocities further west. The velocity of 5.5 km/s is reasonable for the landward portion of the accreted wedge which contains very consolidated sediments lying below and against the Crescent terrane. This is also consistent with the density model from gravity anomaly data (Dehler and Clowes, 1992). This part of the accreted wedge may represent the older component of this accretionary wedge (Tabor and Cady, 1978; Hyndman et al., 1990). (5) From the velocity model, the subducting oceanic plate is gently dipping from west end of the model (at a depth of about 6 km) to the east end beneath the Vancouver Island (at a depth of about 26 km). Also the angle of dipping of the oceanic plate increases from 3-4 degrees at the west end to 10-15 degrees starting at the edge of the shelf to beneath the west coast of Vancouver Island. These are consistent with the corresponding reflection data and previous studies (Hyndman, 1995). The velocities of the oceanic crust are only partly resolved due to limited data coverage, but the velocities of the oceanic Chapter 5. DISCUSSION AND CONCLUSION 74 Deformation 100' 1 1 1 ' ' ' 1 1 1 1 1 ' ' > ' 1 1 • • 1 1 i : 1 1 » i • 0 50 100 150 200 250 300 Distance (km) Figure 5.6: Seismicity cross section from a 100 km wide corridor across the southern Vancouver Island margin from the deformation front to just seaward of the volcanic arc (from Hyndman, 1995). upper mantle are better estimated at velocities of 8.0-8.4 km/s, which are quite consistent with the density model from gravity data (Dehler and Clowes, 1992). (6) The final velocity model of Line 89-02/03 indicates a possible position for the updip limit of the seismic rupture zone for a potential great thrust earthquake in the area. This further indicates the width of total seismic rupture zone would be reduced significantly; therefore the maximum size of the potential great thrust earthquake could be substantially limited. In general, we have developed new velocity-depth structural sections for a part of the margin, particularly below the continental shelf and slope. This would also help to calculate better and more accurate earthquake locations in this region. A seismicity Chapter 5. DISCUSSION AND CONCLUSION 75 cross section from a 100 km wide corridor across the southern Vancouver Island margin is shown in Figure 5.6 (Hyndman, 1995). Both the positions (epicentres) and the depths are well determined for events beneath Vancouver Island and adjacent continental area, but the accuracy decreases rapidly seaward to the west, especially the depth accuracy (Hyn-dman, 1995). This is because these earthquake locations are determined from land-based stations, and depend largely on the velocity structure of the area. Thus an improved velocity structure model would help to improve the accuracy of earthquake locations and to assist in the interpretation of regional tectonics. In this case (Figure 5.6), the new velocity structure from this study would help to increase the accuracy of earthquake locations seaward to the west, especially beneath the eastern continental shelf. References [1] Adams, J. 1990. Paleoseismicity of the Cascadia subduction zone: evidence from turbidites off the Oregon-Washington margin, tectonics, 9: 569-583. [2] Atwater, B .F . , Stuiver, M . , and Yamaguchi, D .K . 1991. A radiocarbon test of earth-quake magnitude at the Cascadia subduction zone. Nature (London), 353: 156-158. [3] 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. [4] Braile, L .W. 1973. Inversion of crustal seismic refraction and reflection data. Journal of Geophysical Research, 78: 7738-7744. [5] Clague, J .J . , and Bobrowsky, P.T. 1994. Evidence for a large earthquake and tsunami 100-400 years ago on western Vancouver Island, British Columbia. Quaternary Re-search, 41: 176-184. [6] Clague, J .J . , and Bobrowsky, P.T. 1995. Tsunami deposits beneath tidal marshes on Vancouver Island, British Columbia. Geological Society of America Bulletin, 106: 1293-1303. [7] Cerveny, V . , Molotkov, I., and Psenak, I. 1977. Ray method in seismology. Charles University Press, Prague. [8] Clowes, R . M . , Brandon, M.T . , Green, A . G . , Yorath, C .J . , Sutherland Brown, A . , Kanasewich, E.R., and Spencer, C. 1987a. L I T H O P R O B E - southern Vancouver 76 .References 77 Island: Cenozoic subduction complex imaged by deep seismic reflections. Canadian Journal of Earth Sciences, 24: 31-51. [9] Clowes, R . M . , Yorath, C.J . , and Hyndman, R.D. 1987b. Reflection Mapping ac-ross the convergent margin of western Canada. Geophysical Journal of the Royal Astronomical Society, 89: 79-84. [10] Clowes, R . M . , Baird, D.J . , and Dehler S.A. 1996. Crustal structure of the Cascadia subduction zone, southwestern British Columbia, from potential field and seismic studies. Canadian Journal of Earth Sciences, 34: 317-335. [11] Dehler, S.A., and Clowes, R . M . 1992. Integrated geophysical modelling of terranes and other structural features along the western Canadian margin. Canadian Journal of Earth Sciences, 29: 1492-1508. [12] Dragert, H . , Hyndman, R.D. , Rogers, G.C., and Wang, K . 1994. Current deform-ation and the width of the seismogenic zone of the northern Cascadia subduction thrust. Journal of Geophysical research, 99: 653-668. [13] 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. [14] Hyndman, R.D. , Yorath, C.J . , Clowes, R . M . , and Davis, E .E . 1990. The north-ern Cascadia subduction zone at Vancouver Island: seismic structure and tectonic history. Canadian Journal of Earth Sciences, 27: 313-329. [15] Hyndman, R.D. , and Wang, K . 1993. Thermal constraints on the zone of major thrust earthquake failure: the Cascadia subduction zone. Journal of Geophysical References 78 Research, 98: 2039-2060. [16] Hyndman, R.D., and Wang K. 1995. The rupture zone of Cascadia great earthquakes from current deformation and the thermal regime, Journal of Geophysical research, 100: 22,133-22,154. [17] Hyndman, R.D. 1995. The Lithoprobe corridor across the Vancouver Island conti-nental margin: the structural and tectonic consequences of subduction. Canadian Journal of Earth Sciences, 32: 1777-1802. [18] Isbell, B., 1989. Summary Report of Land Based Refraction Experiment. Technical Report, University of British Columbia, Vancouver, B.C. [19] Isbell, B., 1990. The Layman's Approach to Digitizing FM Analog Tapes: A Prac-tical Guide. Technical Report, University of British Columbia, Vancouver, B.C. [20] 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. [21] Ludwig, J.W., Nafe, J.E., and Drake, C.L., 1970. Seismic Refraction, in The Sea, vol. 4, 53-84, ed. Maxwell, A.E., Wiley, New York. [22] 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. [23] Massey, N.W.D. 1986. The Metchosin Igneous Complex, southern Vancouver Island: ophiolite stratigraphy developed in an emergent island setting. Geology, 14: 602-605. [24] Phadke, L.C., and Kanasewich, E.R. 1990. Seismic tomography to obtain velocity gradients and three-dimensional structure and its application to reflection data on References 79 Vancouver Island. Canadian Journal of Earth Sciences, 27: 104-116. [25] Rogers, G.C. 1988. An assessment of the megathrust earthquake potential of the Cascadia subduction zone. Canadian Journal of Earth Sciences, 25: 844-852. [26] Savage, J . C , Lisowski, M . , and Prescott, W . H . 1991. Strain accumulation in western Washington. Journal of Geophysical Research, 96: 14,493-14,507. [27] Scholtz, C H . 1990. The mechanics of earthquakes and faulting, Cambridge Uni-versity Press, New York. [28] Shouldice, D.H. 1971. Geology of the western Canadian continental shelf. Canadian Society of Petroleum Geologist Bulletin, 19: 405-436. [29] 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. [30] Spence, G.D., Hyndman, R.D. , Lanton, S., Yorath, C.J . , and Davis, E .E . 1991. Multichannel seismic reflection profiles across the Vancouver Island continental shelf and slope. Geological Survey of Canada, Open File 2391. [31] Tabor, R.W., and Cady, W . M . 1978. The structure of the Olympic Mountains, Wash-ington - analysis of a subduction zone. United States Geological Survey, Professional Paper 1033. [32] Waldron, D.A. , Clowes, R . M . , and White, D.J . 1989. Seismic structure of a subduct-ing oceanic plate off western Canada. In studies of laterally heterogeneous structures using seismic refraction and reflection data. Editedby A . G . Green. Geological Survey of Canada. Paper 39-13, pp. 105-114. References 80 [33] Webring, M. 1985. SAKI: A Fortran program for generalized linear inversion of gravity and magnetic profiles. US Geological Survey, Open File Report 85-122. [34] Zelt, 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. [35] Zelt, C.A., and Smith, R.B. 1991. Seismic traveltime inversion for 2-D crustal veloc-ity structure. Geophysical Journal International, 108: 16-34. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items