UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Crustal structure of the northern Juan de Fuca plate McClymont, Alastair F. 2004

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

Item Metadata

Download

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

Full Text

C R U S T A L S T R U C T U R E O F T H E N O R T H E R N J U A N D E F U C A P L A T E By Alastair F . McClymont B. Sc. (Hons.), Victoria University of Wellington, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA March 2004 © Alastair F . McClymont, 2004 Library Au thor i za t ion In presenting this thesis in partial fulfillment 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. Name of Author (please print) Date (dd/mm/yyyy) Title of Thesis: C\f\L$\o\ ^T^McWt. JAJE\ Ai>AkeXv\ Degree: M- , Y e a r ; 0.00 M-Departmentof £ g / V U gv\(\ Qcc?a<\ Cc[ev\Ce%> The University of British Columbia Vancouver, BC Canada Abstract Anomalous crustal structure of the northern Juan de Fuca plate is revealed from wide-angle seis-mic and gravity modelling. A 2-D velocity model is produced for refraction line II of the 1980 Vancouver Island Seismic Project (VISP80). The refraction data were recorded on three ocean bottom seismometers (OBSs) deployed at the ends and middle of a 110 km line oriented parallel to the deformation front of the Cascadia subduction zone and immediately south of the Nootka fault zone. The velocity model is constructed via ray tracing and conforms to travel time ob-servations of direct, converted and reflected phases. Synthetic seismograms are produced from calculated first-arrival phases for comparison with the recorded seismogram sections. Travel time picks for turning ray phases with shot offsets of > 60 km indicate significantly different lower crustal and upper mantle velocity structure toward each end of the profile. The variation in travel times is modelled by a pronounced increase in igneous crustal thickness of the Juan de Fuca plate, from ~ 7 km at the southern part of the model to ~ 10 km at the northern part (a thickness increase of 43%). A complementary 2-D gravity model using the geometry of the velocity model and velocity-density relationships characteristic of oceanic crust is produced. The densities required to match the gravity field indicate the presence of serpentinized peridotites rather than excess gabbroic crust. Anomalous travel time delays and unusual reflection characteristics observed from proximal seismic refraction and reflection experiments suggest a broader zone of serpentinized peridotites that coincides with the trace of a pseudofault. This correlation implies that the inferred partial serpentinization of the upper mantle may be a consequence of slow-spreading at the tip of a propagating rift. ii T a b l e o f C o n t e n t s A b s t r a c t i i L i s t o f T a b l e s v i i L i s t o f F i g u r e s v i i i A c k n o w l e d g e m e n t s x i Dedication xii 1 I N T R O D U C T I O N 1 1.1 The 1980 Vancouver Island Seismic Project 2 1.2 Tectonic Setting 2 1.2.1 The Explorer triple junction 4 1.3 Formation of Oceanic Lithosphere . 8 1.4 Composition of the Oceanic Lithosphere .9 1.4.1 Marine sediments 9 1.4.2 Igneous basement 10 1.4.3. Upper mantle . . . 11 1.5 Propagating Rift Pseudofaults 12 1.6 Previous Work 12 1.7 Thesis Objectives . . . 15 2 A N A L Y S I S O F R E F R A C T I O N D A T A 17 2.1 Data Acquisition and Processing 17 iii 2.2 First Arrivals: Interpretation and Picking 19 2.2.1 Characteristics of Pc phase 21 2.2.2 Characteristics of Pn phase 22 2.3 Secondary Arrivals: Interpretation and Picking 23 2.3.1 Characteristics of PPS phase 24 2.3.2 Characteristics of PSS phase 24 2.3.3 Characteristics of PPP 1 phase 25 2.4 Picking and Uncertainties 26 2.5 Summary 29 3 MODELLING OF REFRACTION DATA 30 3.1 RAYINVR Model Setup 30 3.1.1 Velocity model parameterization 30 3.1.2 Ray tracing through the model 32 3.2 Developing the 1-D Starting Model 33 3.3 Constructing the 2-D Model ; 35 3.3.1 2-D forward modelling 37 3.3.2 Inversion subsequent to forward modelling 38 3.3.3 Incorporating secondary phase observations 41 3.3.4 Synthetic seismograms 44 3.4 Features of the Velocity Model 49 3.4.1 S-wave velocities and Poisson's ratios . 50 3.5 Model Uncertainties 52 3.5.1 Alternative models 52 3.5.2 Validating with amplitude data 54 3.5.3 Consequences of thick crust and low mantle velocities 56 iv 3.6 Summary 57 4 GRAVITY MODELLING 60 4.1 Data acquisition and processing 60 4.1.1 Characteristics of the gravity field . 61 4.1.2 Previous gravity studies 62 4.1.3 Resolution of the grid data 63 4.2 Development of the gravity model 64 4.2.1 Data procedures and modelling method . 66 4.2.2 Rock densities 66 4.3 Initial model 67 4.3.1 Source of the calculated anomaly 69 4.4 Serpentinization of peridotite mantle 70 4.4.1 Incorporating a serpentinized body 71 4.4.2 Reducing the misfit of short wavelength variations . . . 74 4.5 Summary 75 5 DISCUSSION AND CONCLUSIONS 77 5.1 Supporting Evidence 77 5.1.1 The 1977 Nootka fault zone wide-angle experiment results 77 5.1.2 The 1985 Juan de Fuca plate reflection experiment results 80 5.1.3 Comparison with free-air gravity map . 82 5.2 Serpentinization of Peridotite in Ocean Basins 84 5.2.1 Thermal regime of oceanic lithosphere 85 5.2.2 Serpentinization models 86 5.3 Timing of Serpentinization . 88 5.3.1 Crustal Age Data 88 v 5.3.2 Propagating rift history on the northern Juan de Fuca plate 89 5.3.3 Effects of ridge propagation on the crustal formation process 90 5.3.4 Implications of the initiation of the Nootka fault zone 91 5.4 Serpentinization model for the northern Juan de Fuca plate 93 5.4.1 Are propagating rift events elsewhere associated with serpentinization of the upper mantle? 96 5.5 Conclusions . . . 97 Bibliography 9 9 Appendices 107 A Calculation of gravity response from 2 - D modelling 107 B Conductive cooling model of the oceanic lithosphere 110 vi List of Tables 3.1 Layer P-wave velocities and thicknesses for the 1-D starting model 35 3.2 TRMS and xl residual travel times of first arrivals for 1 -D subsediment velocity model 36 3.3 TRMS and xl residual travel times of first arrivals for the 2-D velocity model before and after inversion routine applied 40 3.4 TRMS and xl residual travel times of primary and secondary phase observations for the 2-D velocity model from OBS-2 airgun data. 45 3.5 P-wave attenuation parameters assigned to the velocity model 46 3.6 Layer thickness, velocity and Poisson's ratio parameters for the final velocity model 53 3.7 TRMS and xl travel time residuals of first arrivals for the favoured 2-D velocity model and alternative models 55 4.1 Layer thickness, density and Poisson's ratio parameters for the initial gravity model. . . 68 4.2 TRMS and xl residual traveltimes of first arrivals for the 2-D velocity model before and after reconciliation with the density model 74 vii List of Figures 1.1 Location map for VISP 80 experiment 3 1.2 Juan de Fuca microplate system 4 1.3 Oblique perspective of the northern Juan de Fuca plate 5 1.4 The Explorer triple junction 7 1.5 Schematic of a spreading ridge and the structure of normal oceanic crust. . . . 9 1.6 Isochron and fracture zone data for the Juan de Fuca plate 14 2.1 VISP II refraction survey profile . 18 2.2 Vertical component receiver gathers for OBS-1 and OBS-4 20 2.3 Vertical component receiver gather for OBS-2 . 22 2.4 Pc phase on OBS-4 receiver gather 23 2.5 Pn phase observations on OBS-1 (left) and OBS-4 (right) 24 2.6 Multiple phases recorded on OBS-2 receiver gather 25 2.7 Terminology of converted and reflected raypath phases 26 2.8 OBS-4 hydrophone receiver gather . 27 2.9 PPP picks and uncertainties for OBS-4 vertical component record section . . . 28 3.1 RAYINVR velocity model parameterization 31 3.2 Initial 1-D subsediment velocity model and forward model results 34 3.3 Pn rays traced through the 2-D velocity model 38 3.4 Velocity model changes resulting from application of the inversion routine. . . 41 3.5 Observed primary and secondary arrivals from OBS-1 and OBS-2 record sec-tions and those calculated from the 2-D model 42 viii 3.6 Observed primary and secondary arrivals from OBS-4 record sections and those calculated from the 2-D model 43 3.7 Observed primary and secondary arrivals from OBS-2 airgun source record sec-tions and those calculated from the 2-D model 44 3.8 Observed data for OBS-1 shown with corresponding synthetic data calculated from the 2-D model . 47 3.9 Observed data for OBS-4 shown with corresponding synthetic data calculated from the 2-D model 48 3.10 Velocity model produced from ray tracing and synthetic seismogram analysis. 49 3.11 S-wave velocity structure established from P- to S-wave converted phase ob-servations 52 3.12 Alternative 2-D velocity model: normal oceanic crustal thickness. . . . . . . 54 3.13 Alternative 2-D velocity model: uniform upper mantle velocities 55 3.14 Observed data for OBS-4 shown with corresponding synthetic data calculated from 3 alternative 2-D models . 59 4.1 Free-air/Bouguer gravity field over coastal and offshore British Columbia. . \ 61 4.2 3-D perspective of offshore British Columbia: free-air gravity draped over ba-thymetry 63 4.3 Gravity station locations and the contoured free-air gravity anomaly within the survey area 65 4.4 Gravity response for density model calculated from the best fit velocity model. 69 4.5 Velocity versus density for serpentinized peridotites 72 4.6 Gravity response for density model modified to include a body of serpentinized peridotite 73 4.7 Density changes resulting from adjustment of layer boundaries 75 ix 4.8 Gravity response for final density model. 76 5.1 Proximal seismic experiments 78 5.2 2-D velocity model of Au and Clowes (1982) 79 5.3 Seismic reflection section of line 85-07 from Hasselgren and Clowes (1995) . 81 5.4 Free-air gravity map overprinted with observed low velocity regions 83 5.5 Synthetic and observed free-air gravity anomaly for profile EX3 84 5.6 Serpentine stability regions for conductively cooling oceanic lithosphere. . . . 86 5.7 Crustal age map for the northern Juan de Fuca plate . . 89 5.8 Reorientation of the Juan de Fuca ridge system: 5 Ma to 3 Ma . . . . . . . . 90 5.9 Schematic of a propagating rift tip 92 5.10 Basement topography beneath the VISPII profile. . . 93 .5.11 Model for serpentinization of the northern Juan de Fuca plate 95 A. l Geometry of a two-dimensional polygon 107 B. l Schematic for model of conductively cooling lithosphere 110 B.2 Conductive cooling profiles for oceanic lithosphere 112 x Acknowledgements I wish to thank Ron Clowes for his guidance and support over the duration of this study. I am also indebted to Phil Hammer and Kim Welford for many helpful discussions and their endless reserves of patience while providing impromptu technical support. To my friends/flatmates Pat Hayman and Chad Hewson, thanks for helping me adapt to Ca-nadian life and for doing the dishes on a regular basis. Finally, I wish to thank Taryn Fay whose companionship and encouragement has enriched my studies and my life here in Canada. xi To my family xii Chapter 1 INTRODUCTION The Juan de Fuca plate has undergone extensive investigation since the inception of plate tec-tonic theory in the 1960s. In this region, the striped character of marine magnetic anomalies was originally discovered [ Raff and Mason, 1961], providing the impetus for the sea-floor spread-ing model. Subsequent wide-angle marine seismic surveys undertaken on the Juan de Fuca plate and elsewhere have enhanced our understanding of how oceanic plates form because they have shown that oceanic crust exhibits the same basic velocity structure everywhere. Variations of the normal model do occur in several anomalous geological settings in ocean basins. These include fracture zones, seamounts, elevated regions surrounding hot spots and the extremely young crust under active spreading centres. From an extensive compilation of seismic refraction surveys within the Pacific and Atlantic oceans, White et al. [1992] have reported that the igneous section of oceanic crust averages 7.1 ± 0.8 km in thickness and exhibits a characteristic one-dimensional velocity structure away from anomalous geological settings. Geophysical and petrological investigations of crustal rocks dredged from ocean basins and ophiolite sequences on land also have shown that normal ig-neous crust has a characteristic composition and exhibits a predictable velocity-density rela-tionship [ Salisbury and Christensen, 1978; Carlson and Raskin, 1984]. The normal velocity model of oceanic crust is used as a template to analyse a wide-angle seismic data set collected over the northern Juan de Fuca plate. If Juan de Fuca crust from this region was formed at a normal spreading ridge setting, then a velocity model constructed from the data should match the model of normal oceanic crust. In addition, a density model converted 1 Chapter!. INTRODUCTION 2 from the velocity model using standard velocity-density relationships should agree with the ob-served gravity field. If either the seismic or gravity data sets show no correlation with the model of normal oceanic crust, an anomalous geological setting will be implicated in the formation of northern Juan de Fuca oceanic crust. 1.1 The 1980 Vancouver Island Seismic Project The 1980 Vancouver Island Seismic Project (VISP 80) involved a series of refraction and re-flection seismic experiments collectively designed to elucidate structure to upper mantle depths within the Cascadia subduction zone (figure 1.1). Complete details of the survey are available in Ellis et al. [1983]. The experiment was the first major seismic study across the continental margin of western British Columbia and Washington, from the oceanic Juan de Fuca plate to the inland volcanic arc. Most of the results have been published and have led to the refinement of Riddihough's [1979] original structural model for the Cascadia subduction zone. Spence et al. [1985] presented a structural model for the complete 350 km onshore-offshore line I. Line II (VISP II hereafter) was positioned to provide control on the crustal structure of the Juan de Fuca plate seaward of the subduction trench. Further details of the survey design will be presented in Chapter 2. An interpretation of VISP II has yet to be published and this study will present an independent analysis of the refraction data. 1.2 Tectonic Setting The Juan de Fuca microplate system lies between 40 ° N and 52 ° N latitude along the Pacific coast of the North American continent. This system is a remnant of the former Farallon plate which began to fragment at about 55 Ma [ Atwater, 1989]. Since that time, a complex history of plate fragmentation has continued in response to the approach of the spreading ridge system on the Cascadia subduction zone. The present configuration is divided into three distinct plates: Chapter 1. INTRODUCTION 3 Figure 1.1: Location map for VISP 80 experiment showing the refraction lines (I-IV), reflection line (RL), and the air gun and CSP lines. Site OBS 6 was occupied for line III and OBS 6-2 for line I. Bathymetry is in metres. Reproduced from Ellis et al. [1983]. the larger Juan de Fuca plate bounded to the north and south by the smaller Explorer and Gorda plates, respectively (figure 1.2). The Juan de Fuca ridge system generates new crust to the Pacific and Juan de Fuca plates and spreads at rates of between 4 and 6 cm/yr . The ridge system has remained at roughly the same distance offshore from the continental margin since at least 6 Ma [ Riddihough and Hyndman, 1991]. The Juan de Fuca plate subducts beneath North America at the Cascadia subduction zone with a convergence rate of between 4 and 5 cm/yr (figure 1.2). The oldest magnetic anomalies Chapter 1. INTRODUCTION 4 on the Juan de Fuca plate have been dated at about 8 to 9 Ma [ Davis and Currie, 1993]. Figure 1.2: Juan de Fuca microplate system. Current plate motions for the Juan de Fuca, Pacific and North American plates relative to an assumed fixed "hot spot" framework are displayed with velocities in cm/yr [ DeMets et al, 1994]. The dotted line defines the study area displayed in figures 1.3 and 1.4. 1.2.1 The Explorer triple junction The tectonics within the Explorer plate to the north are complex (figures 1.3, 1.4). Within this region the Cascadia subduction trench, Juan de Fuca Ridge system and Queen Charlotte Fault meet at a triple junction [ Riddihough and Hyndman, 1989]. The northern terminus of Juan de Fuca ridge is offset from right-lateral transform motion of Queen Charlotte Fault by a series of interconnected spreading centres, the Sovanco fracture zone and minor transform faults bound-ing the western edge of Explorer plate. Chapter 1. INTRODUCTION 5 Queen Charlotte ; Transform zone — A — Deformation front • OBS on refraction line Figure 1.3: Oblique perspective of topography of the northern Juan de Fuca plate for the region outlined in figure 1.2. Overprinted are major tectonic boundaries for the Juan de Fuca (JDF), Explorer (EXP), North American (NAM) and Pacific (PAC) plates. Also shown is the location of VISP II. Whether there is any component of subduction between the Explorer and North American plates is not clear. As crustal age decreases toward the north of Explorer Plate, where spreading centres approach the continental margin, the increased buoyancy greatly diminishes subduc-tion driving forces at the Cascadia subduction zone and increases frictional contact between the subducting plate and the North American continent [ Davis and Currie, 1993]. Since 4 Ma, inhibited subduction within the northern Explorer region has forced the Explorer plate to fragment from, and rotate clockwise relative to, the Juan de Fuca plate [ Riddihough, 1984; Wilson, 1993]. Independent motion of the Explorer plate since 4 Ma is marked by successive clockwise rotation of the Explorer ridge system. Successive rotation brings the Explorer spread-ing centres into an orientation perpendicular to the North American plate boundary and reorients the transform faults parallel to the southern end of the Queen Charlotte Fault [ Braunmiller and Nahelek, 2002]. Kreemer et al. [1998] use strain rates derived from earthquakes to identify a single Explorer Chapter 1. INTRODUCTION 6 transform zone within Explorer plate. Their model substantiates Rohr and Furlong's [1995] pseudo-plate model, in which the Pacific-North American plate boundary is migrating south-ward through Explorer plate such that the new triple junction locus is establishing itself at north-ernmost Juan de Fuca Ridge. They argue that a diffuse northwest band of seismicity linking the Queen Charlotte Fault and northern Juan de Fuca Ridge represents a new transform zone. East of this boundary, Explorer plate has ceased subducting beneath North America; west of this boundary, Explorer plate and Explorer Ridge are being transferred to the Pacific plate (figure 1.4). Braunmiller and Nahelek [2002] relocated several magnitude 5 + epicentres using a joint epicentre determination technique. Their updated locations define a narrow band of seismicity extending from the northern end of Explorer Ridge to the southeast end of the Sovanco fault zone. The relocated earthquakes are not mapped in figure 1.4 but their interpreted transform zone, which represents an actively forming plate boundary between the Explorer and Pacific plates is displayed. Source parameters from select earthquakes are used to calculate instantan-eous rotation poles for the Explorer, Juan de Fuca and Pacific plates. Their plate models predict some component of convergence between the Explorer and North American plates increasing from ~ 0.5 cm/yr near the Tuzo Wilson Knolls to ~ 2 cm/yr at the Nootka fault zone. This con-vergent motion may indicate continued subduction of Explorer plate beneath North America, from the Nootka fault zone to as far north as the Brooks Peninsula (figure 1.4). S-wave velocity structure beneath Vancouver Island established from teleseismic data sug-gests that the northern limit of subducted oceanic plate lies beneath the Brooks Peninsula which is ~ 100 km north of the extension of the Nootka fault zone [ Cassidy et al., 1998]. The fragment of subducted plate observed north of the Nootka fault zone extension is interpreted to represent Explorer plate that was subducted prior to the initiation of the fault zone at 4 Ma. Chapter 1. INTRODUCTION 1 132"W 131"W 130"W 129T/V 128'W 127TW 126'W 125TW 124TW 132-W 131'W 130*W 1291/V 128-W 127TW 126"W 125"W 124TW Figure 1.4: The Explorer triple junction accommodates motion between the Pacific (PAC), North American (NAM) and Juan de Fuca (JDF) plates. Recent (1980 to 2002) epicentres for magnitude 4 + earthquakes are shown from a catalogue maintained by the Geological Survey of Canada. The red band represents Braunmiller and Nabelek's [2002] interpretation of the new Pacific-North America plate boundary initiating within the Explorer (EXP) plate. The new transform boundary implies that the locus of a revised transform-ridge-transform triple junction is establishing itself at the northern tip of the Juan de Fuca Ridge. Yellow bands highlight the Sovanco and Nootka fault zones. The blue and green bands define the area beneath North Amer-ica of the subducted Explorer and Juan de Fuca plates respectively [ Cassidy et al, 1998]. Chapter 1. INTRODUCTION 8 Differential motion between the Explorer and Juan de Fuca plates is largely accommodated by left-lateral motion across the Nootka fault zone, first identified by Hyndman et al. [1979]. The zone is shown from recent seismicity to extend from the northern tip of the Juan de Fuca ridge system to the Cascadia deformation front and beneath the continental shelf to Vancouver Island (figure 1.4). From seismic reflection profiling across the zone defined from seismicity, Hyndman et al. [1979] report basement uplift of several hundred metres and sediment deform-ation extending to the surface. 1.3 Formation of Oceanic Lithosphere The basic structure of oceanic lithosphere is relatively simple compared to that of continen-tal lithosphere. Unlike the complicated histories for accretion of the continents, oceanic crust forms by the same basic process at spreading centres across the globe. New oceanic crust is formed by decompression melting of asthenospheric mantle welling up beneath diverging plates at spreading centres (figure 1.5). That is, the mantle is assumed to rise adiabatically beneath the ridge and undergoes pressure-release melting when it rises above the intersection of the adia-batic temperature with the solidus. Petrological studies show that magmas at spreading ridges have remarkably consistent basaltic compositions. Rare earth element concentrations from mid-ocean ridge basalt (MORB) indicate that the source melts originate at depths of 70 km or less beneath a spreading ridge [ White et al, 1992]. The lithosphere cools and loses buoyancy as it moves away from the spreading ridge and the shallow basaltic crust may experience changes in physical properties with age. These include a predicted decrease in porosity of the shallow basaltic crust in response to lithostatic loading and the precipitation of hydrothermal minerals within void spaces. The underlying lithospheric mantle also cools but the depth at which man-tle rocks readily deform increases with age away from the ridge. Consequently the thickness of the rigid plate moving over deeper mantle increases with distance from the spreading axis. Chapter 1. INTRODUCTION 9 Figure 1.5: Schematic of a mid-ocean spreading ridge. The inset shows the generalized struc-ture and composition of normal oceanic crust derived from marine seismic refraction studies [ Ewing and Houtz, 1979] and petrological studies of ophiolite complexes [ Salisbury and Chris-tensen, 1978]. The seismic layers represent changes in velocity gradient and do not imply a change in rock composition. 1.4 Composition of the Oceanic Lithosphere A number of recent studies provide insight into the regional structure of the Juan de Fuca plate and these observations are incorporated into the velocity and density models. These include physical properties determined through geophysical investigation and lithological information gained by direct sampling methods. The gross constituents of oceanic lithosphere are a blanket of marine sediments, igneous basement crust and upper mantle. 1.4.1 Marine sediments High erosion rates from mountain belts on the North American continent generate high rates of sedimentation over the Cascadia subduction zone and adjacent Juan de Fuca Plate. This regional sedimentary trap is termed the Cascadia Basin. The dominantly terrigenous turbidite sediments trapped within the Cascadia Basin acquire thicknesses as great as 1-2 km [ Au and Clowes, 1982; Chapter 1. INTRODUCTION 10 Ellis et al, 1983]. Drill core samples, recovered 100 km east of Juan de Fuca Ridge crest, show two distinct periods of sedimentation. Initial deposition of hemipelagic muds and fine turbidites is followed by a Pleistocene sequence of variably thick silt- and sand-turbidites interbedded with hemipelagic muds [ Underwood and Hoke, 2000]. P-wave velocities of deep marine sediments change rapidly with compaction and vary from velocities close to that of sea water (~ 1.5 km/s) at the seabed to about 3.5 km/s above igneous basement [ Jones, 1999]. 1.4.2 Igneous basement Pioneering marine seismic refraction studies (e.g. Hill [1957]) first demonstrated that, world-wide, oceanic crust is divisible into 2 layers: an upper layer with P-wave velocities of the range 4.0-6.0 km/s and a lower layer with velocities between 6-7 km/s. These were termed Seismic Layers 2 and 3 respectively. Advances in survey design led to higher definition of oceanic struc-ture. Layer 2 was further subdivided into three layers (2A, 2B and 2C) while two distinct ve-locity zones were,identified in Layer 3 (3A and 3B) [ Ewing and Houtz, 1979]. The general structure and composition of the igneous basement is summarized in figure 1.5. Much of what we know of the composition of the upper igneous basement comes from drilling. The Layer 2A/2B boundary is interpreted as the contact between basaltic lavas and a sheeted dyke complex [ Jones, 1999]. The steep velocity gradients observed in the upper oceanic crust are understood to result from closure of cracks and pore spaces with depth. High porosities in the shallow crust, created by rapid cooling and contraction, are reduced with depth by lithostatic loading, crystallization of hydrothermal minerals within void spaces and intrusion of igneous bodies. Land based evidence from ophiolite sequences indicates that the lower crust (Seismic Layer 3) is composed of plutonic diorites and gabbros and is underlain by mantle peri-dotites [ Salisbury and Christensen, 1978]. Layer thicknesses and velocities also vary laterally with age from spreading centres. At very Chapter 1. INTRODUCTION 11 close distances to the ridge, seismic velocities and densities within the shallow basaltic crust in-crease with age in response to lithostatic loading and a decrease in porosity. However, Johnson and Semyan [1993] show evidence that for longer time scales (> 10 million years), the shal-low basaltic crust may exhibit a systematic increase in porosity and a corresponding decrease in seismic velocity and rock density. Regardless of this interpretation, the variability of their measurements shows that the shallow crust (approximately equivalent to seismic layer 2A) ex-hibits significant variation in both seismic velocity and density. Results from a variety of seismic refraction surveys on oceanic crust away from anomalous geological settings show that the mean thickness of igneous crust is 7.1 (± 0.8) km [ White et al, 1992]. 1.4.3 Upper mantle Remnants of oceanic upper mantle are found in outcrop around the world, including ophiolite sequences (e.g. Salisbury and Christensen [1978]). The exposed rocks have an olivine-rich peridotite composition. Mantle and cummulate peridotites are often exposed in slow-spreading ridges on fault escarpments and transform faults. These peridotites are often altered by hydra-tion to form serpentinites. Upper mantle P-wave velocities exceed 7.6 km/s. Modelling results from a variety of ma-rine seismic refraction studies show that velocity gradients beneath the Moho are small with velocities reaching 8.4 km/s at depths of about 55 km [ Ewing and Houtz, 1979]. Velocities of the range 8.0-8.2 km/s are typically used to model the mantle near the base of the crust. Seismic anisotropy can cause large variation in P-wave velocity with changing azimuth of wave propagation. Keen and Barrett [1971] showed that upper mantle P-wave velocities in the northeast Pacific are fastest parallel to sea floor spreading directions (~ 8.3 km/s) and slowest perpendicular to them (~ 7.8 km/s). Seismic anisotropy within the upper mantle is attributed Chapter 1. INTRODUCTION 12 to the alignment of olivine crystals in response to convective mantle flow. The VISP II profile is oriented approximately 45 ° from the Juan de Fuca Ridge spreading direction(figure 1.3). Therefore, a sub-Moho P-wave velocity intermediate of the two extremes reported in the north-east Pacific is expected along strike of VISP II. If the direction of mantle flow in the study area varied little at the time of formation, the sub-Moho P-wave velocities sampled beneath VISP II from forward and reversed arrays should be relatively uniform. 1.5 Propagating Rift Pseudofaults Recent changes in relative plate motion between the Juan de Fuca and North American plates are recorded by the disrupted pattern of marine magnetic anomalies observed astride the Juan de Fuca ridge (figure 1.6). The current trend of the Juan de Fuca ridge and that of the older magnetic lineations on the Pacific plate indicate that the ridge has undergone 20 ° of clockwise rotation since 9 Ma [ Davis and Currie, 1993]. Wilson et al. [1984] proposed that the rotation of segments of the ridge system took place by a series of propagating rift events. According to the propagating rift model [ Hey, 1977], new spreading centres develop or-thogonal to the new spreading direction, growing and replacing the rift associated with the old spreading pattern. As the rift propagates through older lithosphere, the juxtaposition of new crust against existing crust results in the development of V-shaped discontinuities within the magnetic anomaly pattern which are termed pseudofaults (figure 1.6). The implications of rift propagation for the formation of Juan de Fuca crust will be discussed further in section 5.3.3. 1.6 Previous Work Numerous marine seismic surveys have been undertaken over the Juan de Fuca microplate sys-tem but most have focussed on determining crustal structure across tectonic boundaries at the margins of the plates. Chapter 1. INTRODUCTION 13 Of the refraction surveys that traverse segments of the spreading ridge axes, White and Clowes [1990] and Cudrak and Clowes [1993] present results over Endeavour segment of the Juan de Fuca Ridge. Both studies interpret crustal sections broadly similar to the model of normal oceanic crust and find no evidence for a low velocity anomaly consistent with an axial magma chamber. However, Menke et al. [2002] image a low velocity zone consistent with a shallow-crustal magma chamber beneath Coaxial segment on a more southerly part of the Juan de Fuca ridge system. Similar evidence for a shallow magma chamber beneath Explorer Ridge is inter-preted by Malecek and Clowes [1978]. A seismic survey over Axial volcano at the intersection of the Cobb hotspot chain with the Juan de Fuca Ridge by West et al. [2003] modelled 1-2 km of excess crust which they attribute to excess melt supply from the coalascence of two separate melting features. The trench axis (characteristic of most subduction zones) at the eastern margin of the Juan de Fuca plate is masked by thick sediments of the Cascadia Basin. However, several seismic studies have probed beneath the sediments and interpreted models [ Ellis et al, 1983; Clowes et al, 1987; Hyndman et al, 1990] show the classic arrangement of accretionary wedge sediments over subducting lithosphere. Waldron et al. [1990] reviewed the seismic data from Line 1 of the VISP 80 experiment (figure 1.1) and incorporated nearby drill hole data and results from an intersecting seismic reflection profile to produce a refined 2-D velocity model for Line 1. Their interpreted model exhibits a normal oceanic crust with igneous thickness of ~ 7 km for Juan de Fuca plate seaward of the trench. Chapter 1. INTRODUCTION 14 Figure 1.6: Isochron and fracture zone data derived from the magnetic anomaly reversal se-quence for the Juan de Fuca plate (reproduced from Wilson [1993]). Dots show individual cross-ings of isochrons measured from ship tracks. Solid lines are the interpreted locations of pseu-dofaults. Age scale at lower left shows the age in Ma of the selected isochrons, identified by grey level, with the full reversal sequence in black. Inset diagram shows the propagating rift model of Hey et al. [1988]. Grey bands: lithosphere formed by rift propagation; black bands; lithosphere produced at the doomed rift; diagonal stipple: transferred lithosphere. Chapter 1. INTRODUCTION 15 Au and Clowes [1982] produced 1-D velocity models from a seismic refraction survey ac-ross the Nootka fault zone. Their models show a significant amount of lateral variation of the igneous crustal thickness, from 6.4 to 11 km. The authors suggest that the abnormally thick crust may have resulted from adjustments at the spreading ridge in response to the initiation of the new plate boundary. The structure of Juan de Fuca crust immediately west of line VISP II is revealed via two seismic reflection lines and interpretations made by Hasselgren and Clowes [1995]. Continu-ous reflections attributable to oceanic Moho suggest that the majority of crust beneath the two profiles formed with a uniform thickness that is representative of normal oceanic crust. The results of this survey and those from the study of Aa and Clowes [1982] will be discussed fur-ther in section 5.1. VISP II extends immediately south of the Nootka fault zone and does not traverse a plate boundary. In the absence of preliminary evidence suggesting an anomalous geological setting,. it is expected that crust forming this part of the Juan de Fuca plate should fit White et al. 's [1992] model of normal oceanic crust. 1.7 Thesis Objectives Initial analysis of the seismic refraction data was made by H. White andR. Clowes [unpublished study, 1990]. Their preliminary interpretation indicated significant crustal thickness changes for the northern Juan de Fuca plate. However, these changes are not reflected in the gravity field. An independent analysis of the seismic data and reconciliation with the gravity field provides the impetus for this study. More specifically, the objectives of this thesis are to: 1) Identify characteristic phase arrivals from the marine refraction data that can be used to construct a velocity model. Chapter 1. INTRODUCTION 16 2) Apply ray tracing techniques to construct a velocity model for the northern Juan de Fuca plate that conforms to the traveltime picks and amplitude data. 3) Convert the velocity model to a density model using velocity-density relationships that satisfy the observed gravity data. 4) Use the velocity and density models to interpret the crustal structure and composition of the northern Juan de Fuca plate. 5) Relate interpretations from the models to the past and present plate tectonic regime. Chapter 2 ANALYSIS OF REFRACTION DATA The two-dimensional velocity structure of the Juan de Fuca plate, parallel to the Cascadia sub-duction zone, is interpreted from the marine seismic refraction data set. The timing and pat-tern of recorded phases reveals information on the crustal velocity structure with depth. The identification of characteristic raygroups and careful observation of arrival times is essential to producing an accurate velocity model. 2.1 Data Acquisition and Processing The marine seismic refraction data were collected as part of the 1980 Vancouver Island Seis-mic Project (VISP 80). Profile VISP II (figure 2.1) was positioned to provide insight into the structure of the oceanic crust immediately south of the Nootka fault zone on the northernmost part of Juan de Fuca plate. The refraction data were recorded on three ocean bottom seismometers (OBSs) deployed at the ends and middle of a 110 km line oriented parallel to the deformation front of the Casca-dia Subduction Zone. Explosive charges of varying size were fired at ~2.5 km intervals into all OBSs while 32 L air gun sources at 0.5 km spacings were recorded principally on OBS-2. Large 500 kg explosive charge sizes were used at the ends of the profile to ensure sufficient transmission of energy at up to 110 km offset. A 100 kg charge was fired at the centre of the profile and smaller (25-50 kg) charges were detonated at ~2.5 km intervals between the large charges. The arrivals at each OBS were recorded on vertical and horizontal geophone compo-nents and a supplemental hydrophone. Additional single channel, continuous seismic profiling 17 Chapter 2. ANALYSIS OF REFRACTION DATA 18 was carried out using a 5 L air gun to determine the structure of the sediments and to provide two-way travel times to basement [ Ellis et ai, 1983]. 128°W 127°W 49°N 48°N \ OBS-4 (v, hyd) V v. v « k OBS-2 (v,h) \ \ Vv \ —— Explosive source coverage Airgun source coverage V OBS-1 (v, h) \ im-49°N km 50 100 48°N 128°W I27°W Figure 2.1: VISP II refraction survey profile. The individual component record sections uti-lized in this study are listed beside each OBS: (v) vertical, (h) horizontal, (hyd) hydrophone. Bathymetry contours are in metres. At the time of acquisition, timing and positional systems were less developed than they are at present. For example, the explosive source shots were off-loaded from the shooting ship and fuse detonated. The origin time for the shot was calculated from the arrival time of the direct water wave back to the ship. Origin time errors arising from the uncertainty of the ship's speed may be as large as ± 15 ms for the smaller charges and ± 35 ms for the large charges. Additional Chapter 2. ANALYSIS OF REFRACTION DATA 19 timing error was introduced due to the time drift of the internal OBS clocks. Ellis et al. [1983] calculate a maximum time error of ± 60 ms for synchronization of the shot and receiver for all ocean bottom seismometers used in the VISP 80 experiment. The original unprocessed data were recorded on various magnetic tapes, all of which are now missing or, after 20 years of deterioration, unreadable. The present study is reliant on data extracted from hard copies of the plotted seismogram sections. Unfortunately, the hard copies available do not represent a complete set of the recorded data. Figure 2.1 illustrates the data, coverage used in this study. The seismogram sections were originally compiled by digitizing the analogue OBS signals with sampling rates in the range 163 to 230 samples/s. The signals were bandpass filtered (1-10 Hz) and seismograms were amplified by a factor of r2 to compensate for geometric spreading with distance from the source. All the sections were plotted as receiver gathers with a reducing velocity of 6.0 km/s to ac-centuate first-arrival phases propagated through the igneous section of the crust. These were the hard copies available for this study. However, a digital format was required such that the sections could be manipulated to enhance individual phases and assist the picking process. The hard copy data were redigitized into SEG-Y format with sampling rates in the range 400 to 500 samples/s. Once digitized, the SEG-Y file headers were updated with survey geometry and tim-ing information for input into the various PLOTSEC seismic processing routines [ Amor, 1996]. 2.2 First Arrivals: Interpretation and Picking The data quality is variable, ranging from high signal-to-noise ratios (S/N hereafter) for most receiver gathers at near offset to poor at far offsets. Figure 2.2 illustrates the data quality for OBS-1 and OBS-4 receiver gathers which together represent reversed sections of the VISP II profile. The seismograms are scaled to each trace's maximum root-mean-square (RMS) amp-litude. This choice of amplitude scaling best enhances the first refracted P-wave arrivals. Chapter 2. ANALYSIS OF REFRACTION DATA 20 o 0 I E H •ml i f ? Offset Distance (km) o VD i Offset Distance (km) Figure 2.2: Vertical component receiver gathers for OBS-1 and OBS-4. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). First P-wave arrivals are best observed from the vertical component data. As the wave prop-agates through the basement/sediment interface, an abrupt change from relatively high velocity basement to low velocity sediment is encountered. Consequently, the upcoming ray is refracted to near-vertical incidence beneath each OBS. First arrivals are easy to pick where the preceding trace is relatively noise free. Where it becomes difficult to distinguish the first recorded signal from background noise, picks are made on the basis of waveform correlation between traces. The observed first P-wave arrivals are termed PPP. Two distinct phases are interpreted from PPP: rays turning within the crust (classified in this study as Pc) and those that turn within the Upper Mantle (Pn). Pc is well observed out to shot offsets of around 60 km on each receiver Chapter 2. ANALYSIS OF REFRACTION DATA 21 gather. The Pn phase arrivals appear to coincide with the degradation in S/N beyond this dis-tance (figure 2.2). There are few near offset (< 10 km) traces available from the explosive source data. First arrivals from OBSs 1,2 and 4 receiver gathers exhibit velocities of ~6.0 km/s. Near offset observations are available from the OBS-2 airgun sections from which a phase with velocities between 2-4 km/s, typical of marine sediments, is visible. The sedimentary phase ob-servations provide a limited indication of the velocity structure of the marine sediments beneath OBS-2. Results from continuous seismic profiling show significant variation in the thickness of the sedimentary blanket across the profile and suggest a laterally variable velocity structure. Short wavelength variations of Pc and Pn phases probably result from lateral velocity varia-tions within the sedimentary layer. The travel time path through the sedimentary layer can vary considerably as a result of basement controlled topography, creating variable sediment thick-ness. A difference of 200 m in sediment thickness at different distances along the profile, as-suming an average sediment velocity of 2 km/s, will produce at least 0.1 s difference in travel time through the sedimentary layer. 2.2.1 Characteristics of Pc phase Oceanic mantle velocities are typically greater than 7.5 km/s [ Jones, 1999] and I have used this criteria to discrminate Pc from Pn. Disregarding short wavelength variations, general Pc phase trends are observed from cursory inspection of the data. At 25 km shot offset, OBS-1 receiver gather first arrivals exhibit an abrupt change in velocity from ~ 6.0 km/s to ~7.0 km/s (figure 2.2). This velocity discontinuity is not evident on OBS-4 record section but this may be a consequence of pick quality deteriorating with shot offset. Velocities greater than 7.0 km/s do not appear on OBS-2 and OBS-4 vertical component record sections until at least 35 km and 55 km shot offset, respectively. OBS-2 receiver gather displays traces for shots fired both north and south along the profile Chapter 2. ANALYSIS OF REFRACTION DATA 22 (figure 2.3). First arrival times show little difference out to offsets of about 35 km on either side of the OBS. However, between 35 and 40 km offset north of OBS-2 the apparent first ar-rivals indicate a rapid velocity increase from ~ 7.0 km/s to ~ 8.0 km/s, as would be observed for a change from a crustal phase (Pc) to a mantle phase (Pn). However, this evidence for the appearance of Pn at such close shot offset is tenuous as it is based on picks from just 2 traces. Offset Distance (km) Figure 2.3: Vertical component receiver gather for OBS-2. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). Pc on OBS-4 is clearly observed to about 40 km shot offset. Between 40 and 60 km offset, first arrivals are characterized by one cycle of low amplitudes preceding 2-3 cycles of higher amplitudes (figure 2.4). Beyond 60 km offset, where S/N is greatly diminished, the transition from Pc to Pn is difficult to pick. 2.2.2 Characteristics of Pn phase Pn appears on all of the explosive source vertical component receiver gathers (OBS-1,2 and 4) beyond 50 km shot offset. It is characterized by low S/N and a velocity in the range of 7.6 to 8.1 km/s. OBS-1 and OBS-4 provide the best record of Pn as they represent forward and reversed sections of the entire profile. The most striking observation of Pn from OBS-1 and OBS-4 is Chapter 2. ANALYSIS OF REFRACTION DATA 23 OBS-4 South (Vertical) LL > 4J Low amplitude j first arrivals • Distance (km) Figure 2.4: Pc phase picks (circles) between 30 km and 70 km shot offset on OBS-4 vertical component receiver gather. a clear discrepancy in Pn arrival times (figure 2.5). Although difficult to resolve, at 110 km offset, OBS-1 first arrivals appear at around -0.5 s (reduced travel time), while first arrivals on OBS-4 appear to arrive at around 0.5 s. These differences imply significantly different velocity structure toward either end of the profile. 2.3 Secondary Arrivals: Interpretation and Picking A number of secondary phases are observed from both the vertical and horizontal component receiver gathers. These phases result from rays converting to a different phase or reflecting at discontinuous boundaries before arriving at the OBS. Secondary phases identified from the re-ceiver gathers are classified using terminology from a previous interpretation of the data [H. White and R. Clowes, unpublished study, 1990]. The phases of interest are best observed from the OBS-2 airgun sections (figure 2.6). Raypaths for these phases are illustrated in figure 2.7. Chapter 2. ANALYSIS OF REFRACTION DATA 24 Offset Distance (km) Offset Distance (km) Figure 2.5: Pn phase observations on OBS-1 (left) and OBS-4 (right). Seismogram traces for shots fired at > 50 km offset are displayed. OBS-4 record section is plotted with a reversed orientation for direct comparison with OBS-1. The onset of Pn from Pc is interpreted to coincide with a decrease in S/N. First break picks are plotted as a continuous grey line. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). 2.3.1 Characteristics of PPS phase The PPS phase (figures 2.6, 2.7) results from the P- to S- phase conversion at the base of the sediment layer beneath each OBS. The ascending ray encounters large velocity and density con-trasts at this boundary and the resulting S-wave phase is clearly observed on all the horizontal component receiver gathers. The low velocity sediments refract the upcoming ray to near ver-tical incidence such that the S-wave branch of the ray path differs little with shot offset. Conse-quently, the PPS phase is observed with a constant time delay of 1.5 to 2 seconds after the first arrival. The uniformity of the delay across each section provides an additional check on first arrival picks. 2.3.2 Characteristics of PSS phase The PSS phase conversion occurs where the descending ray encounters the sediment/basement boundary and then travels as a shear phase for the remainder of its path. This phase is not de-tected from the explosive source data because the sections become noisy with reverberations Chapter 2. ANALYSIS OF REFRACTION DATA 25 Offset Distance (km) Figure 2.6: Multiple phases recorded on OBS-2 receiver gather from the airgun sources. Shots are recorded on both vertical (top) and horizontal (bottom) components. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). after the PPS arrival; however, it is visible on both components of the airgun sections and ap-pears between 12 km and 27 km shot offset (figure 2.6). The PSS arrivals exhibit a velocity of ~ 3.5 km/s corresponding to the S-wave velocity of the shallow igneous basement. 2.3.3 Characteristics of PPP1 phase PPP1 is the reflection from the free water surface back to the OBS. This water-bottom mul-tiple is most apparent from the hydrophone records of the explosions, of which only OBS-4 hydrophone section remains (figure 2.8). PPP1 can also be observed from the airgun data and it is most prominent on the vertical component between 20 km and 30 km shot offset (figure 2.6). Chapter 2. ANALYSIS OF REFRACTION DATA 26 Figure 2.7: Simplified geometry and terminology of converted and reflected raypath phases ob-served on the OBS-2 (airgun) sections. Bold solid lines represent P-wave raypaths and dashed lines represent S-wave raypaths. The PPP phase classification encompasses both Pc and Pn phases. PPP1 arrives with a constant time delay (~ 3.5 seconds) after PPP, varying with the different water depths at OBS-2 and OBS-4. 2.4 Picking and Uncertainties All picks were made using the PLOTSEC PICK utility [ Amor, 1996]. Picking difficulties aris-ing from low S/N on some of the seismograms require the correlation of phase waveforms from trace to trace. Noise in this instance includes wave train reverbarations subsequent to the first arrival. The 'trace flattening' utility is employed to aid identification of the phase from each seismogram to the next. Traces are levelled to a picked datum and the waveforms immediately following the picks are inspected to ensure correlation. Some error may still be introduced where the first few cycles of a phase are indiscernible Chapter 2. ANALYSIS OF REFRACTION DATA 27 — ao — T O — eo — B O — 4 0 — 3 0 —ao — 1 0 Offset Distance (km) Figure 2.8: OBS-4 hydrophone receiver gather. Shear waves are absent because the hydrophone measures pressure differences within water. The PPP1 phase appears with a constant time delay of ~ 3.5 seconds corresponding to the two way travel time in the water column above OBS-4. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). from background noise. Where such a pick would provide redundant information, the trace is ignored. The far offset (> 60 km) traces of the OBS-1 and OBS-4 record sections provide useful information on the deep crust and mantle and first arrival picks are included, albeit with large uncertainties. If the trace is relatively noise free and the phase break is clearly observed, an uncertainty corresponding to 1/4 of a waveform cycle is assigned. This criteria accounts for the inability to pick exactly where the first break occurs on the plotted seismogram and is used for most first arrivals on traces within 40 km shot offset. If a trace exhibits significant noise, the signal wave-form may become smeared and distinguishing 1/4 of the waveform cycle is difficult. In this case, uncertainties are calculated on the basis of signal-to-noise ratio. The PLOTSEC AMPPK routine f Amor, 1996] calculates pick uncertainty by comparing the energy within a specified time window before and after the pick. The S/N values are divided into bins to which uncer-tainties are assigned by the user. A S/N value considered low is allocated an uncertainty of 1/4 Chapter 2. ANALYSIS OF REFRACTION DATA 28 of its perceptible waveform cycle. The largest uncertainty assigned using this process is +/-0.150 seconds. This is a subjective technique but ensures that uncertainties are allocated in a consistent fashion. The PPP phase is most difficult to pick on the OBS-4 record section at > 60 km shot off-set. Beyond 60 km offset it is possible to pick several phases potentially attributable to PPP and the differences in pick times for each trace are as great as 0.6 seconds (figure 2.9). Hence there are two types of pick error that need to be accounted for: the uncertainty of the phase and the uncertainty of the pick for that phase. First break picks for these traces are calculated as the midpoint between the earliest and latest picked phases. The phase uncertainty is the time between the earliest and latest phase picks of each trace. An additional pick uncertainty is al-located on the basis of S/N over this time window (+/- 0.150 seconds for most of these traces) and this is added to the phase uncertainty to collectively produce a total first break uncertainty for each seismogram. _ no — 1 0 0 -oo -ao - 7 o -oo -oo Distance (km) Figure 2.9: Far offset PPP picks and uncertainties for OBS-4 vertical component record section plotted with a reducing velocity of 7.8 km/s. Where more than one coherent phase has been picked, circles indicate the earliest and latest pick on each trace. Horizontal lines represent the midpoint between these picks and vertical, thick grey bars define the combined uncertainty of the phase and the pick for each phase. Chapter 2. ANALYSIS OF REFRACTION DATA 29 2.5 Summary First arrival phases are picked from receiver gathers covering the entire length of the VISP II profile for all three OBSs. Secondary phase information is more limited but two converted phases and a water-bottom multiple phase are identified. Of the first arrival phases, the crustal (Pc) phase exhibits large amplitudes and velocities < 7.6 km/s while the mantle (Pn) phase shows low amplitudes and velocities > 7.6 km/s. Decreasing S/N with offset away from the receivers makes the Pn phase difficult to pick. A l l of the picked phases are used for comparison with predicted travel time data which is calculated by ray tracing through a velocity model. The allocation of an uncertainty of the ap-propriate size for each pick is crucial for it determines how well parts of the predicted model should fit the corresponding data. The picks are assigned variable uncertainties which generally increase with shot offset. The modelling process can highlight picks that appear out of place. Some picks were reviewed after modelling had commenced but care was taken not to introduce bias for the sake of a favourable model. Chapter 3 MODELLING OF REFRACTION DATA Crude seismic velocity models are easily calculated from refraction data if simple ray path geome-tries are assumed. That is, the subsurface is composed of layers separated by plane and possibly dipping interfaces, seismic velocities are uniform within each layer and layer velocities increase with depth. Complex models that are based on more realistic assumptions require the imple-mentation of sophisticated computer algorithms such as RAYINVR [ Zelt and Smith, 1992]. 3.1 RAYINVR Model Setup The forward modelling component of the RAYINVR program is based on asymptotic ray theory in two-dimensional media [ Zelt and Ellis, 1988]. In the forward modelling process, rays are traced from a designated source through a parameterized velocity model to an array of receivers. Travel times calculated for the rays are compared to those observed at the receivers. Model pa-rameters are updated by the user and the algorithm repeated until a sufficient match is obtained. A full description of the forward modelling theory used in RAYINVR is available in Zelt and Ellis [1988]. 3.1.1 Velocity model parameterization RAYINVR uses trapezoids, of variable size, to model two-dimensional velocity structure. A model constructed of a few simple trapezoids, as opposed to one built from a large matrix of cells of uniform velocity, has a greatly reduced number of variable parameters. This choice of parameterization makes ray tracing more efficient and reduces the time required to update the 30 Chapter 3. MODELLING OF REFRACTION DATA 31 model. Greater complexity is introduced by adding more trapezoids. Node values are input by the user as points on layer boundaries and RAYINVR automatically divides the model into the required number of trapezoidal blocks (figure 3.1). Undefined node values on each layer bound-ary are calculated by linear interpolation between those values specified by the user. RAYINVR determines undefined node parameters by interpolating between defined node parameters such that each trapezoid is defined by 4 nodes each with an x and z coordinate and a P-wave velocity. The velocity along upper and lower trapezoidal boundaries is determined by linear interpolation between nodes. The velocity at any point within a trapezoid is determined by linear interpola-tion between the upper and lower boundary along a vertical path. The method of interpolation allows for horizontal and vertical velocity gradients within each trapezoid. Velocity disconti-nuities are incorporated into the model by specifying different velocities at the bounding nodes on either side of a layer boundary. x (distance) Figure 3.1: Example of the RAYINVR velocity model parameterization. This three-layer model is defined by 20 independent model parameters: 10 boundary nodes (squares) and 10 velocity nodes (solid circles). RAYINVR automatically divides the model into 8 trapezoidal blocks. The velocity at any point within a block is calculated by linear interpolation. A layer boundary has a constant velocity on that side where only 1 velocity node is specified. Modified from Zelt and Smith [1992]. Chapter 3. MODELLING OF REFRACTION DATA 32 3.1.2 Ray tracing through the model Two-dimensional ray tracing equations, which describe the path of a propagating ray through the velocity model, are a pair of first-order ordinary differential equations. These can be written in two forms: — = — , 8 6 = (v'~v*^e\ (3.1) 8x tan#' 8x v ' 8x „ 89 (v, tan 9 — vx) — = tan0, j - ^ ^ 1 (3.2) oz dz V with initial conditions x = x0, Z = Z0, 9 — 90 6 is the angle between the tangent to the ray and the z axis [ Zelt and Ellis, 1988]. The variable v is the wave velocity and vx and vz are partial derivatives of velocity with respect to the x and z coordinates, respectively (z is positive downward). These equations are solved by Runge-Kutta numerical integration [ Sheriff and Geldart, 1983]. System 3.1 is solved with x as the integration variable when the ray path is near-horizontal, and system 3.2 is solved with z as the integration variable when the ray path is near-vertical. Snell's law is applied where a ray intersects a layer boundary. The ray step length, A, used to solve 3.1 or 3.2 varies depending on the local derivative of the velocity field. It is adjusted at each point along the ray path according to the relationship av A = — — (3.3) \Vx \ + \Vz\ The constant a is specified by the user. The appropriate choice of a is a compromise be-tween accuracy and computational time. Equation 3.3 yields a small ray step where the local velocity gradient is high so that the proceeding ray path is calculated accurately where the ray is bending sharply. Chapter 3. MODELLING OF REFRACTION DATA 33 The total travel time t between a source and receiver along a ray path for a velocity model is the sum of the travel times for each ray segment: where n is the total number of ray segments, which will vary depending on both the choice of a and the complexity of the velocity model, and k and Vi are the path length and velocity of the ith ray segment, respectively. Refracted, reflected and head waves may be traced through the velocity model. Through se-lections made by the user, ray paths may also include multiple reflections and converted phases. The user may specify that rays with a wide range of take-off angles are traced through the model or that the algorithm search for the fastest ray path linking source and receiver. The source or sources may be positioned anywhere within the model but receivers are assumed to be located at the top of the model. Consequently the desired survey configuration, where multiple shots at the sea surface are fired into one receiver on the seabed, is excluded. To overcome this prob-lem, the shots are specified as receivers, the OBSs are specified as shots and rays are traced in reverse. 3.2 Developing the 1-D Starting Model The initial velocity model is built from one-dimensional structure and two-dimensional com-plexity was subsequently added to satisfy observations from each OBS. The 1-D model is di-vided into 9 layers (including 1 layer for sea water), representing the typical layer-cake struc-ture of oceanic lithosphere established from worldwide marine refraction surveys and studies of ophiolite complexes (e.g. Ewing and Houtz [1979], Salisbury and Christensen [1978]). A simple model is constructed that approximately conforms to OBS-2 first arrival picks (figure 3.2). The velocity model is essentially 1-D but incorporates the 2-D geometry of the seafloor (3.4) Chapter 3. MODELLING OF REFRACTION DATA 34 and basement topography defined from the continuous seismic profiling (section 2.1). Table 3.1 lists layer thicknesses and P-wave velocities used in the 1-D starting model. Rays are traced through 2-D model space that is 130 km in length and extends to 25 km in depth. Velocity (km/s) 1 3 <i 3 9 7 9 1, •3 1-D igneous crust and upper mantle OBS-1 Offset Distance: 0 O B S - 2 Offset Distance: -80 OBS-4 Offset Distance: -120 1 40 60 Model Distance (km) 20 40 60 80 -60 -40 -20 0 20 -100 -80 -60 -40 -20 100 120 40 0 Figure 3.2: Initial 1-D subsediment velocity model is shown on the left. On the right are the ob-served first arrival picks, including their estimated error, from the explosive source shots to each OBS and those calculated by RAYINVR for the velocity model on the left. Picks and calculated arrival times are colour coded for each OBS: OBS-1 (red), OBS-2 (green) and OBS-4 (blue). Inverted triangles represent the position of each OBS relative to profile distance. A reducing velocity of 6.0 km/s is applied. Previous compositional and structural studies of the Juan de Fuca Plate, in addition to ubi-quitous global observations of 1-D velocity structure, contribute a priori information that is included in the starting model (outlined in section 1.4). Two layers are used to represent a bimodal sedimentary sequence in the starting velocity model: an upper layer with a velocity of 1.8 km/s overlying a layer with a velocity of 2.3 km/s. These velocities are consistent with those used by Ellis et al. [1983] in their interpretation of continuous seismic profiling (CSP) data along the VISP I profile and those used by Au and Chapter 3. MODELLING OF REFRACTION DATA 35 Layer a (km/s) AZ(km) Sea water 1.49 2.53 - 2.60 Sediments (1A) 1.80 0.41 - 0.63 Sediments (IB) 2.30 0.36-0.91' 2A 3.50 - 4.80 0.03 - 0.87 2B 4.80 - 5.90 0.80 2C 5.90 - 6.20 0.50 3A 6.50 - 6.90 2.50 3B 7.10-7.50 2.50 Upper mantle 7.90 - 8.20 14.00 Table 3.1: Layer P-wave velocities (a) and vertical thicknesses ( A 27) for the 1 -D starting model shown in figure 3.2. The upper mantle layer extends to a maximum model depth of 25 km and does not represent the true thickness of oceanic upper mantle. A P-wave velocity of 8.20 km/s is fixed at 25 km depth. Clowes [1982] to model refraction data across the nearby Nootka fault zone. The total thick-ness of these layers varies with the topography of the underlying basement. This topography is constrained using two-way travel times from CSP data (see section 2.1). The geometry of the two sedimentary layers remained unchanged throughout the modelling process; however, velocities were changed to incorporate lateral variation. 3.3 Constructing the 2-D Model Travel times generated from the initial model satisfy most observations from OBS-2 and most OBS-1 picks from the northwestern end of the profile. The model yields a poor fit to most OBS-4 observations and OBS-1 picks from the southeastern end of the profile. Overall, calculated travel times for OBS-1 picks are uniformly late and for OBS-4 picks they are uniformly early. Before proceeding to improve the model, it is necessary to establish a measure of how well any improved model fits the data. Observed travel time picks are of variable quality and the degree to which parts of the model are updated to fit each pick should be weighted appropriately. Chapter 3. MODELLING OF REFRACTION DATA 36 The Reduced Chi Squared (xl) statistical test provides the best measure of the goodness of fit for the model derived data because it takes into account measurement uncertainty of the observed data. It is defined as: (3.5) where T{ is the difference between the time observed and the time calculated for each pick (the travel time residual), <J is the uncertainty of the individual pick and n is the total number of cal-culated travel times. A xl value of 1 represents the optimal fit to observations. A xl value well in excess of 1 (e.g. 10) implies either 1) that the model is a poor fit to the travel time picks, 2) that the picks themselves are poor or, 3) that the assigned uncertainties are overly optimistic. A value much less than 1 (e.g. 0.1) indicates either that the model contains unnecessary hetero-geneity or that pick uncertainties are too pessimistic. The xl and the root-mean-square (TRMs) of travel time residuals for the Pc and Pn phases of the 1-D subsediment model from the three OBSs are displayed in table 3.2. Phase TRMS (S) n Pc 0.121 6.271 125 Pn 0.415 2.614 35 Pc and Pn 0.223 5.432 161 Table 3.2: T R M S and xl residual travel times of explosive source first arrivals to the three OBSs for the 1-D subsediment velocity model displayed in figure 3.2. Note that the large xl value for the Pc phase indicates that the crustal velocity structure is a poor fit to observations. Chapter 3. MODELLING OF REFRACTION DATA 37 3.3.1 2-D forward modelling RAYINVR allows rays to be traced in groups, layer by layer. The modelling process, of incor-porating 2-D structure into the igneous crust, involved first adjusting the geometry of the shal-lowest layer (seismic layer 2A) by adding nodes at the lower boundary, and then calculating travel times for rays traced through that layer. Once a satisfactory fit to observed travel times was achieved, the process was repeated for the next layer. Observed picks from the airgun data, fired within 30 km offset of OBS-2, provided the best constraint on the thickness of seismic layers 2A and 2B. Velocity node values remained unchanged for seismic layers 2A, 2B and 2C because slight variations in the thickness of these layers produced a satisfactory fit to observa-tions. The fit to observations within 60 km shot offset is readily improved by introducing small lateral heterogeneity to the shallow seismic layers 2A, 2B, 2C and 3A. Observations beyond 60 km offset, specifically from OBS-1 and OBS-4, present a more difficult problem as illustrated by figure 3.2. Beyond 60 km offset on OBS-1 (70-120 km model distance in figure 3.2), calculated arrivals are slightly late and require either a decrease in the thickness of layers 3A and 3B or, increased velocities for layers 3A, 3B and the upper mantle. Modelled arrivals for the same shot offset on OBS-4 (10-60 km model distance in figure 3.2) require either a significant increase in the thickness of layers 3 A and 3B or, reduced velocities for layers 3 A, 3B and the upper mantle. A compromise is required that goes some way to satisfying observations from both OBSs. A combination of increased layer thickness and reduced seismic velocity was applied to lay-ers 3B and the upper mantle at the northwestern end of the model to minimize deviation from the normal thickness of oceanic crust. Delaying arrival times through this part of the model significantly improved the degree of fit to OBS-4 picks but increased the misfit to OBS-1 picks between 100 and 110 km shot offset. Some improvement was made to the OBS-1 picks, with-out significantly degrading the fit to OBS-4 picks, by increasing the velocities within layer 3A Chapter 3. MODELLING OF REFRACTION DATA 38 at the northwestern end of the profile. Despite improving the fit to lower crustal and upper mantle phase picks from OBS-1 and OBS-4, modelled arrival times fall outside pick uncertainties. The best fit to far offset picks (> 60 km) is one which yields calculated arrivals delayed relative to observed arrivals from OBS-1 and too early for those from OBS-4. It is clear from ray tracing through the model constructed with seismic layers that some degree of crustal thickening, between 60 km and 130 km along the model profile, is necessary to satisfy Pn picks from OBS-2 and OBS-4 record sections (figure 3.3). SE Model Distance (km) NW 0 . 20 40 60 80 100 120 1 1 1 1 i J . 0 20 40 60 80 100 120 Model Distance (km) Figure 3.3: Pn rays traced through the 2-D velocity model from OBS-2 and OBS-4 are shown above and the calculated travel times are displayed below. The dashed white line represents the modelled position of the crust-mantle boundary. A reducing velocity of 6 km/s is applied. 33.2 Inversion subsequent to forward modelling RAYINVR incorporates an inversion routine that allows velocity or boundary nodes to be up-dated such that travel time residuals are minimized using a damped least-squares solution [ Zelt Chapter 3. MODELLING OF REFRACTION DATA 39 and Smith, 1992]. This routine was utilized to fine tune the fit of calculated travel times to obser-vations without significantly altering the forward model. No structural or velocity nodes were added to the model. Travel time inversion is a non-linear problem because the travel times are a function of both velocity structure and the raypath geometries, which themselves are velocity dependent. The problem is overcome by using a starting model and an iterative approach. Initial model para-meters are perturbed by a prescribed extent and rays are traced through the perturbed model such that travel time perturbations corresponding to each model parameter can be calculated analyt-ically (see Zelt and Smith [1992]). The results yield a partial derivative matrix containing the elements 8ti/6m.j which describe how the ith travel time changes with the jth model parameter. The linear problem is thus formulated A A m = At , (3.6) where A is the partial derivative matrix, A m is the model parameter adjustment vector , and At is the travel time residual vector. The travel time residual vector is determined by raytracing prior to perturbing the velocity model. After tracing rays through the perturbed model, the vec-tor of model parameter adjustments is found from the combination of the travel time residual vector and the inverse of the updated partial derivative matrix A m = A 1 At . (3.7) An optimal solution to the overdetermined problem is found using damped least-squares and the new model parameter adjustment vector is applied to the previous starting model. Succes-sive iterations of the inversion step may be required until a desired misfit threshold is satisfied. Chapter 3. MODELLING OF REFRACTION DATA 40 The sparse ray coverage of the lower crust and mantle does not lend to a robust model in-version. To avoid introducing unnecessary complexity, boundary nodes remained fixed and ve-locity nodes only were allowed to vary. Velocity nodes were set to vary a maximum of 0.1 km/s for each inversion. An inversion was applied more than once to specified velocity nodes where a significant shift in velocity was required. This was deemed appropriate in shallow regions of the model where there is a high degree of raypath coverage and in deeper regions where a velocity shift greater than 0.1 km/s did not introduce anomalous features. Figure 3.4 illustrates the perturbations made to the velocity model as a result of updating the velocity nodes. The changes are subtle but illustrate a general increase in velocities in the mid to upper crust at the northwestern end of the profile. The slight polarization of the lower mantle, into higher velocities beneath OBS-1 and lower velocities beneath OBS-4, is consistent with the forward modelling approach. The modelling statistics before and after the application of the inversion routine are also displayed in table 3.3. The xl values for the Pc phase subsequent to inversion are significantly reduced and represent an improvement to the velocity model in the mid to upper crust. Conversely, xl values for the Pn phase are marginally reduced and suggest that introducing further complexity into the structure of the lower crust and upper mantle will not significantly improve the model. Before After Phase TRMS (S) n TRMS (S) Xr n Pc 0.074 1.877 153 0.070 1.611 153 Pn 0.284 1.646 36 0.279 1.620 36 Pc and Pn 0.145 1.804 189 0.137 1.604 189 Table 3.3: TRMS and xl residual travel times of first arrivals for the 2-D velocity model con-structed from forward modelling both before and after application of the inversion routine. Chapter 3. MODELLING OF REFRACTION DATA 41 SE NW ° 1 OBS-1 OBS-2 OBS-4 f -0.15 0.00 0.15 Figure 3.4: Velocity model adjustments made by the application of the inversion routine. Posi-tive changes are shown in red and negative changes in blue. The bold dashed line indicates the modelled position of the boundary between lower crust and upper mantle. 3.3.3 Incorporating secondary phase observations Converted and reflected phase picks are included as a further test of the validity of the 2-D model. These include PPS, PSS and PPPI as outlined in section 2.3. PPS is observed on OBS-1 and OBS-2 record sections and is best observed from the hori-zontal component data. It arrives with a uniform delay after PPP providing a good constraint on the S-wave velocity structure of the two sedimentary layers from which Poisson's ratio is determined (outlined in section 3.4.1). In addition, the pattern of arrivals mimics the velocity gradients observed from PPP validating the picks used to model the deep crustal structure. Fig-ure 3.5 displays primary and secondary arrival picks and the travel times calculated from the 2-D model for OBS-1 and OBS-2. The explosive source PPS phase observations for OBS-1 and OBS-2 are well fit by those calculated from the model and yield a combined xl of 2.519 s from 62 calculated travel times. This value represents a fit of slightly lesser quality than that calculated from first P-wave arrivals Chapter 3. MODELLING OF REFRACTION DATA 42 .0 SE NW P P S I'll P P P . . . I 1 , Observed Calculated 1 I M , , OBS-1 Model Distance (km) SE i Observed Calculated P P S 1 i • , , T I t : » ' t ' P P P i , 1 . * * ! ! 1 Pn " I M M O B S - 2 T GO SO Model Distance (km) NW Figure 3.5: Observed primary and secondary arrivals from OBS-1 (top) and OBS-2 (bottom) record sections and those calculated from the 2-D model. The transitions from Pc to Pn within the PPP phase are indicated with black arrows. Line coverage of the OBS-1 horizontal compo-nent record is limited and PPS observations are only picked between 10 km and 90 km profile distance. A reducing velocity of 6 km/s is applied. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). (table 3.3). No horizontal component data were available from OBS-4 from which PPS could be picked. However, the PPP1 multiple was clearly picked from the OBS-4 hydrophone record section (fig-ure 2.8). The observed uniform time delay provides information on water depth only at OBS-4. Much like PPS, velocity gradients observed from the pattern of arrivals are consistent with those observed from PPP. Observed and calculated travel times for the PPP and PPP1 phases from OBS-4 are shown in figure 3.6. PPP1 is also picked from OBS-2 airgun source data between Chapter 3. MODELLING OF REFRACTION DATA 43 15 km and 30 km shot offset to the north (figure 3.7). SE Observed Calculated Pn 40 GO • : 80 Distance (km) NW IT ! ' I PPP1 PPP OBS-4 100 120 Figure 3.6: Observed primary and secondary arrivals from OBS-4 record sections and those calculated from the 2-D model. The transition from Pc to Pn within the PPP phase is indicated with a black arrow. The OBS-4 hydrophone record is limited to within 80 km shot offset and as such PPS observations are only picked between 40 km and 120 km profile distance. A reducing velocity of 6 km/s is applied. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). PPS, PPP1 and PSS were all picked from the OBS-2 airgun data (figure 3.7). The gradient of the PSS phase provides a good measure of the S-wave velocity structure of the upper igneous crust. All of the phases observed from the OBS-2 airgun data are suitably fit by calculations from the model and demonstrate that the P- and S-wave velocity structure of the sediments and upper crust (to ~ 7 km depth) beneath OBS-2 is replicated well by the model velocity layers. The xl values (table 3.4) illustrate a good fit to observations. The low value for PPP1 suggests that uncertainties assigned to picks for this water surface multiple are probably larger than they should be. Chapter 3. MODELLING OF REFRACTION DATA 44 Q I 4 PSS |NW| / 1 SE 1 pss '% \ \ \ f % f t f » % PPPI pps P ^ V ^ ' t l M ^ PPS PPP "* 71}*^ PPP 1 Observed • Calculated O B S - 2 • 1 1 40 60 1 80 1 1 100 120 Model Distance (km) Figure 3.7: Observed primary and secondary arrivals from OBS-2 airgun source record sections and those calculated from the 2-D model. The airgun source PPPI was observed only between 10 km and 30 km shot offset north of OBS-2. A reducing velocity of 6 km/s is applied. The inset shows the featured OBS position and line coverage (shaded black) relative to the rest of the line and other OBSs (shaded grey). The different converted phase arrivals all showed a satisfactory fit to observations and con-sequently no additional alterations were made to the velocity model. 3.3.4 Synthetic seismograms As a secondary check of the validity of the 2-D model, synthetic receiver gather seismograms are produced from the model for comparison with the recorded data. Amplitude corrections were previously applied to each of the recorded traces to account for variation in the size of the explosive source and to compensate for the effect of geometric spreading [R. Clowes, personal communication, 2003]. Corrections made to the amplitudes of the original recorded seismo-grams prevent the comparison of absolute amplitudes between the recorded and synthetic data. Chapter 3. MODELLING OF REFRACTION DATA 45 Phase TRMS (S) Xr n PPP 0.058 1.087 86 PPS 0.084 1.586 74 PSS 0.099 2.034 56 PPP1 0.060 0.281 22 Table 3.4: TRMS and xl residual travel times of primary and secondary phase observations for the 2-D velocity model from OBS-2 airgun data (figure 3.7). However, the relative pattern of amplitude variation with shot offset is not significantly altered by these corrections. The pattern of amplitude variation can be used to distinguish between first-order discontinuities in velocity structure and velocity gradients within the crust. Large velocity gradients cause seismic rays to turn more rapidly than low velocity gradients. As a result, rays propagating through a velocity field with a high gradient will not disperse as rapidly as they would travelling though a field of low velocity gradient and more energy is focussed back to the surface receiver. The amplitude of turning and reflected waves is calculated by zero-order asymptotic ray theory using the TRAMP algorithm [ Zelt and Ellis, 1988]. The amplitude A of the ray at its endpoint is given by the expression: A = (3.8) L where A0 is the initial ray amplitude and L is the geometrical spreading. The factor q accounts for energy partitioning at velocity discontinuities and is calculated using Zoeppritz displace-ment amplitude coefficients [ Zelt and Ellis, 1988]. The factor e is the integer number of times the direction of positive displacement changes for the SV segment of the raypath with respect to the ray coordinate system. This expression assumes a point source with uniform directional characteristics. The time and complex amplitude for all rays traced to each shot location are su-perpositioned. The result is a series of impulses with the appropriate amplitude and phase shift. Chapter 3. MODELLING OF REFRACTION DATA 46 The final seismdgram is produced by convolution with a source wavelet function. The source wavelet applied here is produced by stacking first arrivals from the observed data. A quality factor (C?) is assigned to each layer of the model to account for energy loss by anelastic attenuation. Table 3.5 lists the compressional wave Qp values assigned to each layer of the model. Qp values assumed for the sea water and sediment layers are based on acoustic measurements made in marine sediments by Best etal. [2001]. Attenuation parameters defining layers 2A and 2B are derived from the results of Jacobson and Lewis [1990] while values for the lower crust and mantle are assumed from Spudich and Orcutt [1980]. Layer QP Seawater 1000 Sediment 1 10 Sediment 2 10 2A 25 2B 150 2C 450 3A 450 3B 450 Upper Mantle 450 Table 3.5: P-wave attenuation parameters assigned to the velocity model Synthetic seismograms produced for OBS-1 and OBS-4 (both vertical component, explo-sive source) are shown in figures 3.8 and 3.9: The relative amplitudes of the first arrival syn-thetic data appear to mimic those for the observed sections. The general trend observed from OBS-1 synthetic section (figure 3.8) is a pattern of large amplitudes out to 60 km shot offset fol-lowed by an abrupt change to relatively low amplitudes. The low amplitudes exhibited by the first 3 traces coincide with a shadow zone between rays turning at the base of Seismic Layer 2C and rays turning within uppermost Seismic Layer 3A. This feature is not exhibited in the recorded data and indicates the model velocity gradient in Seismic Layer 2C may be slightly Chapter 3. MODELLING OF REFRACTION DATA 47 too large beneath OBS-1. More importantly, the low amplitudes observed beyond 60 km offset are coincident with the recorded data and indicate a low velocity gradient as would be expected within the upper mantle. 20 40 60 80 100 120 0 20 40 60 80 100 120 Shot Offset (km) 0 20 40 60 80 100 120 M o d e l Dis tance (km) Figure 3.8: Observed data (top) for OBS-1 shown with corresponding synthetic data (middle) calculated from the 2-D model. Both the observed and synthetic sections are plotted with an r2 distance correction factor. The source wavelet is plotted to the same scale as the travel time axis. Pc (blue) and Pn (green) ray coverage for the receiver gathers is shown below. The same change, from high amplitudes to low amplitudes, is apparent from OBS-4 syn-thetic data (figure 3.9) at approximately 80 km shot offset. It is difficult to locate this feature Chapter 3. MODELLING OF REFRACTION DATA 48 with the same accuracy in the recorded data as it is masked by significant noise. Large ampli-tudes are coherent to at least 70 km offset suggesting that the upper mantle, with a low velocity gradient, resides at a greater depth beneath OBS-4 than it does for OBS-1. M o d e l Dis tance ( k m ) Figure 3.9: Observed data (top) for OBS-4 shown with corresponding synthetic data (middle) calculated from the 2-D model. Both the observed and synthetic sections are plotted with an r 2 distance correction factor. The source wavelet is plotted to the same scale as the travel time axis. Pc (blue) and Pn (green) ray coverage for the receiver gathers is shown below. Chapter 3. MODELLING OF REFRACTION DATA 49 3.4 Features of the Velocity Model The final velocity model was updated to include results from synthetic seismogram modelling. The model (figure 3.10) conforms to first arrival, converted phase and amplitude observations. It is constructed to satisfy these observations and to show the minimum possible deviation from the structure of normal oceanic crust. SE NW OBS-1 OBS-2 OBS-4 0 10 20 30 40 50 60 70 80 90 100 110 120 130 Model Distance (km) | | ~| ^ ^ ^ ^ ^ M ~ T " ^ ^ (km/s) 1.48 1.50 4.00 5.00 6.00 6.70 7.20 7.60 8.00 8.40 Figure 3.10: The final velocity model produced from ray tracing and synthetic seismogram _ lysis. Velocity contours are in km/s. The dashed white line represents the inferred position of the crust-mantle boundary (CMB). The region of the model with zero raypath coverage is masked out. Vertical exaggeration is 2.5:1. ana-The model displays little lateral variation in velocity down to 7 km depth. Below 7 km depth, velocity contours indicate a significant horizontal velocity gradient decreasing to the northwest. Thickened crust in the 2-D velocity model (figure 3.10) is manifest by velocities representative of crustal material (3.5-7.6 km/s) extending to a greater depth at the northwestern end of the profile than at the southeastern end. The change in subsediment crustal thickness, from ~ 7 km at the southeastern end to ~ 10 km at the northwestern end, is modelled by collectively enlarging each seismic layer across the profile. This was considered the most appropriate way to model the travel time data because trends identified from the Pc phase could not be attributed Chapter 3. MODELLING OF REFRACTION DATA 50 to individual layers. In the absence of an identifiable PmP phase the boundary between crust and upper mantle is resolved using turning rays only. As outlined in Chapter 2, first arrival picks are classified as Pc or Pn. Boundary nodes defining the crust-mantle boundary (CMB) are manipulated such that the different phases are traced to their correct pick location. At the northwestern end of the profile the CMB is modelled as a first-order discontinuity between lower crustal velocities of < 7.3 km/s and upper mantle velocities of > 7.7 km/s. This discontinuity is less pronounced at the southeastern end of the profile where lower crustal velocities as high as 7.8 km/s were required to fit calculated travel times from far offset rays to arrivals observed at OBS-1. However, the vertical velocity gradient in the upper mantle is significantly greater at the northwestern end of the profile than at the southeastern end. Consequently, the transition from crustal to mantle phases exhibited from the synthetic amplitude data is less pronounced from OBS-4 than it is on OBS-1 (figures 3.8, 3.9). Significant features of the favoured velocity model are an increase in crustal thickness from ~ 7 km to ~ 10 km coupled with decreased upper mantle and lower crustal velocities from southeast to northwest along profile. This interpretation of crustal thickening by collectively, and progressively, enlarging individual seismic layers along profile implies that thick crust is a consequence of an increase in melt generation at the Juan de Fuca ridge. However, the model is not so well constrained that greater structural complexity can be ruled out. 3.4.1 S-wave velocities and Poisson's ratios Poisson's ratio a is an elastic constant which provides an additional physical property meas-urement that can be related to rock composition. It is defined from the ratio of P- to S-wave velocities: Chapter 3. MODELLING OF REFRACTION DATA 51 Hey - 2 2(a/8Y - 2 (3.9) where a is P-wave velocity and 3 is S-wave velocity. The ratio is dimensionless and for crustal rocks can range between 0 and 0.5 (fluid). Most crustal rocks yield Poisson's ratios between 0.25 and 0.30 [ Shearer, 1999]. RAYINVR converts P-wave velocities to S-wave velocities from a Poisson's ratio defined for each layer. The elastic property of the material in each layer is thus assumed to be laterally homogeneous. S-wave velocities are established from arrival times of the P- to S- converted phases. Rather than adjust S-wave velocity values at each individual velocity node, a values for each layer are updated until a satisfactory fit to converted phase observations is achieved. PPS arrivals are observed on all horizontal component record sections covering the entire profile length. They provide extensive observations from which the S-wave velocity structure of the sedimentary layers across the complete profile is defined. However, S-wave velocities for the igneous crust are derived only from PSS observations observed on OBS-2 airgun sections. Ray coverage from these sections is limited to within 30 km shot offset and the maximum model depth for these turning rays is ~ 7 km, which corresponds to Seismic Layer 3A (figure 3.11). As such, the S-wave velocity structure at greater depth and beyond 30 km offset of OBS-2 is unknown. Table 3.6 lists the range of P- and S-wave velocities for each layer together with the assigned values for Poisson's ratio. The ratios determined for the sedimentary layers are very close to 0.5 indicating a high degree of water saturation. Poisson's ratio's for the upper igneous crust match those typical of oceanic basalts [ Jones, 1999]. Chapter 3. MODELLING OF REFRACTION DATA 52 Figure 3.11: S-wave velocity structure established from P- to S-wave converted phase observa-tions. Also shown is the limited ray coverage for the PSS phase from which the S-wave velocity of the igneous basement is determined. The maximum depth of penetration of this ray group is into the shallow part of seismic layer 3A. 3.5 M o d e l U n c e r t a i n t i e s The xl and TRMS values provide a measure of how well the final model fits the observed data relative to other model choices. The favoured model does not satisfy all observations, especially picks beyond 60 km shot offset on the explosive source record sections. The optimal model solution is a compromise yielding travel times too early for OBS-2 and OBS-4 picks and too late for OBS-1 picks. 3.5.1 A l t e r n a t i v e m o d e l s Travel times calculated from the 1-D model (figure 3.2) demonstrate that first arrivals to OBS-2 and OBS-4 picks greater than 60 km shot offset need to be delayed by anomalous velocity struc-ture in the lower crust or upper mantle. The favoured 2-D model achieves this using a combi-nation of thickened crust and low velocities in the upper mantle. It should be noted, however, that the travel time data can be satisfied by variations of this model. Figures 3.12 and 3.13 rep-resent alternative end-member models that produce TRMS and xl values comparable to those Chapter 3. MODELLING OF REFRACTION DATA 53 Layer AZ (km) a (km/s) 0 (km/s) cr Sea water 2.53 - 2.60 1.49 - 0.5 Sediments (1A) 0.41 - 0.63 1.73 - 1.90 0.47 - 0.52 0.46 Sediments (IB) 0.36-0.91 2.28 - 2.39 0.69 - 0.72 0.45 2A 0.27 - 0.99 3.46 - 4.54 1.85-2.43 0.30 2B 0.30-1.30 4.96 - 5.93 2.86 - 3.42 0.25 2C 0.25 - 1.10 5.94 - 6.33 3.43 - 3.65 0.25 3A 2.30 - 3.20 6.19-7.02 3.57 - 4.05 0.25 3B 3.30-4.30 7.04 - 7.80 - -Upper mantle - 7.72 - 8.20 - -Table 3.6: Vertical layer thickness, velocity and Poisson's ratio parameters for the final velocity model (figures 3.10, 3.11). S-wave velocities (J3) and Poisson's ratios ( c r ) are indeterminate for Seismic Layer 3B and the upper mantle. calculated from the favoured 2-D model (Table 3.7). Figure 3.12 displays alternative model 1 (ALT1). Here, the crust is modelled with a relat-ively uniform thickness of ~ 7 km. The upper mantle is assigned an anomalously low velocity of 7.6 km/s at the northwestern end of the profile to produce the required Pn travel time delays. An oceanic upper mantle velocity of 7.6 km/s represents the low extreme of mantle velocities observed worldwide and is normally associated with elevated mantle temperatures typically ob-served beneath spreading centres and plume affected crust [ White et al., 1992]. The localization of low mantle velocities at the north-western end of the profile requires a large lateral velocity gradient of 8.1 km/s to 7.6 km/s from southeast to northwest. The model presented in figure 3.13 (ALT2) uses upper mantle velocities > 8.0 km/s. The upper mantle exhibits a low lateral velocity gradient and Pn travel times are delayed by thicken-ing Seismic Layer 3B. This interpretation requires the total thickness of igneous oceanic crust to increase from ~ 7 km at the southeastern end of the profile to 12.5 km at the northwestern end. White et a/.'s [1992] global compilation of marine seismic survey results indicates that thick crust of this magnitude is exclusively a consequence of mantle plume interaction with the Chapter 3. MODELLING OF REFRACTION DATA 54 SE NW OBS-1 OBS-2 OBS-4 T 1 T 90 100 110 120 130 Model Distance (km) I I I (km/s) 1.48 1.50 4.00 5.00 6.00 6.70 7.20 7.60 8.00 8.40 Figure 3.12: Alternative 2-D velocity model (ALT1). The total igneous crustal thickness is a near uniform 7 k m across the profile. Observed Pn travel time delays are replicated using an-omalously low mantle velocities at the north-western end of the profile. Velocity contours are in km/s. The dashed white line represents the inferred position of the crust-mantle boundary ( C M B ) . The region of the model with zero raypath coverage is masked out. Vertical exaggera-tion is 2.5:1. spreading ridge. 3.5.2 Validating with amplitude data In addition to travel time data, the most appropriate model choice may be determined from amp-litude information. Synthetic seismogram sections were generated from calculated first arrivals for all 3 models. The sections plotted for comparison with the OBS-4 explosive source, vertical component receiver gather exhibit the most variation between models and are displayed in fig-ure 3.14. Chapter 3. MODELLING OF REFRACTION DATA 55 OBS- OBS -2 OBS-4 20 -I 1 1 1 1 1 1 1 1 1 1 1 1 1-0 10 20 30 40 50 60 70 80 90 100 110 120 130 Model Distance (km) 1.48 1.50 4.00 5.00 6.00 6.70 7.20 7.60 8.00 8.40 Figure 3.13: Alternative 2-D velocity model (ALT2). Upper mantle velocities are > 8.0 km/s and exhibit a low lateral velocity gradient. Observed Pn travel time delays are replicated by thickening the igneous crust. The total igneous crustal thickness increases from ~ 7 km to 12.5 km south-east to north-west. Velocity contours are in km/s. The dashed white line represents the inferred position of the crust-mantle boundary (CMB). The region of the model with zero raypath coverage is masked out. Vertical exaggeration is 2.5:1. Model TRMS (S) xl n Favoured 0.138 1.627 189 ALT1 0.144 1.692 188 ALT2 0.133 1.674 187 Table 3.7: TRMs and xl travel time residuals of first Pn arrivals for the favoured 2-D velocity model, displayed in figure 3.10, and the alternative velocity models shown in figures 3.12 and 3.13. Chapter 3. MODELLING OF REFRACTION DATA 56 The synthetic seismogram section produced for model ALT2 displays a first arrival ampli-tude pattern that is a poor match to the observed data. The ~ 0.5 s travel time delay of the largest amplitudes on seismograms between -30 and -45 km shot offset is not evident in the observed data. In addition, prominent waveform amplitudes are calculated to occur out to almost 90 km shot offset and are followed by low waveform amplitudes resulting from mantle phase arrivals. This feature is not visible from the observed section but may be masked by significant noise. The two discrepancies combine to make model ALT2 an unsuitable candidate and suggest that the favoured model may represent the limit to which thick crust can occur beneath OBS-4. Unlike the favoured model and model ALT2, the synthetic data for model ALT1 displays no clear transition from large to small amplitudes attributable to the onset of the Pn phase from the Pc phase with offset away from OBS-4. It should be noted, however, that high S/N makes this characteristic difficult to identify from the real data. Otherwise, the pattern of synthetic data presented for model ALT1 shows reasonable consistency with the waveforms of the recorded section and should not be ruled out on the basis of amplitude information alone. 3.5.3 Consequences of thick crust and low mantle velocities All three models presented adequately fit the travel time data but their anomalous lower crustal and upper mantle velocity structures stray considerably from that of normal oceanic crust. The favoured model and ALT2 implicate the influence of a mantle plume beneath the Juan de Fuca ridge to generate the increased melt required to form thicker than average crust. The thicker crust in model ALT2 necessitates interaction with a more voluminous piume head than that re-quired by the favoured model. Similarly, given that profile VISP II is a significant distance from an active spreading centre, a mantle plume appears to be the only explanation for lower than av-erage upper mantle velocities interpreted from model ALT1 and the favoured model. The low mantle velocities used to model ALT1 are more extreme than that of the favoured model and Chapter 3. MODELLING OF REFRACTION DATA 57 indicate that the plume size would be more of a factor to increase upper mantle temperatures to produce lower than expected velocities. Alternatively, the low mantle velocities could be at-tributed to mantle anisotropy, but the large lateral velocity gradient exhibited by model ALT1 requires a significant change in the direction of mantle flow beneath the crust. Measurements of upper mantle velocities at different shooting azimuths (see section 1.4.3) suggest that a reduc-tion in upper mantle velocity from 8.1 km/s to 7.6 km/s would require a change in the direction of mantle flow from parallel with the orientation of the VISP II profile beneath OBS-1 to per-pendicular to the profile orientation beneath OBS-4. Magnetic anomaly lineations on the Juan de Fuca plate that intersect VISP II do not indicate a dramatic change in orientation as would be expected for such a large change in the direction of sub-crustal mantle flow (figure 1.6). There-fore, mantle anisotropy does not provide a suitable explanation for the low mantle velocities exhibited by model ALT1. Elevated basement topography, formed in response to mantle plume interaction with oceanic crust, is observed at numerous locations worldwide on both large scales (e.g. the Hawaiian Is-lands and the Kerguelen Plateau) and small scales (e.g. the Heck Seamount Chain west of the Juan de Fuca Ridge). Given that the effects of a mantle plume are not apparent in the crustal basement topography observed from CSP data recorded along the VISP II profile, the favoured model is the preferred choice because it minimizes the extent to which a mantle plume would interact with the crust. 3.6 Summary Picked first arrivals are used to construct an initial velocity model using the RAYINVR ray trac-ing algorithm and supplemental inversion routine. Additional constraints are enforced using secondary phase picks and amplitude information. The integration of travel time and ampli-tude data disqualify model ALT2 and imply that the total vertical thickness of igneous crust Chapter 3. MODELLING OF REFRACTION DATA 58 beneath the northwestern end of profile VISP II should be limited to no more than 10 km. The lack of evidence for a mantle plume at the same location suggests that the favoured model is a better candidate than model ALT1 because it requires a smaller plume to generate the low upper mantle velocities interpreted from the data. For these reasons the favoured model is selected to test compatibility with gravity data in the following chapter. Chapter 3. MODELLING OF REFRACTION DATA 59 Offset Distance (km) -120 -100 -80 -60 -40 -20 0 -120 -100 -80 -60 -40 -20 0 Offset Distance (km) Figure 3.14: Observed data for OBS-4 (top) shown with corresponding synthetic data calculated from the 3 alternative 2-D models displayed in figures 3.10, 3.12 and 3.13 respectively. Both the observed and synthetic sections are plotted with an r 2 distance correction factor. The source wavelet is plotted to the same scale as the travel time axis. Chapter 4 GRAVITY MODELLING Compressional P-wave velocity can be expressed as: where A and p are elastic moduli that relate P-wave velocity (a) to rock density (p). The elas-tic moduli are themselves a function of pressure, temperature and mineralogy and these para-meters may be calculated using experimental results. A more practical method of translating velocity to density utilizes empirically derived velocity-density relationships established from direct measurements of seismic velocities and density for particular rock types. In this chapter the velocity model derived from forward modelling of the VISP II seismic re-fraction data and various velocity-density relationships are used to predict 2-D density models. A gravity field is calculated for these models that should match observations from the free-air gravity field. If no match is obtained, an explanation for the mismatch is required. 4.1 Data acquisition and processing Reduced gravity data for western Canada are available from the National Gravity Data Base maintained by the Geophysical Data Centre of the Geological Survey of Canada. The data are stored in grid format with a 2 km grid spacing. Free-air (offshore) and Bouguer (onland) grid data for coastal and offshore British Columbia are shown in figure 4.1. All measurements are reduced to the International Gravity Standardization Network 1971 (IGSN71) datum. The Ge-odetic Reference System 1967 (GRS67) gravity formula is used to calculate theoretical gravity 60 Chapter 4. GRAVITY MODELLING 61 values. Terrain corrections were applied to the Bouguer data. The Geological Survey of Canada quotes the offshore data as accurate to ± 2 mgal. 4.1.1 Characteristics of the gravity field 20 5 0 -5 -10 -15 -20 -30 -35 -45 -60 -70 -80 -90 -100 -105 -110 -115 -120 -125 -130 -135 -145 Gravity (mgal) -130° -128° -126° Figure 4.1: Gravity field over coastal and offshore British Columbia. The free-air anomaly is shown over marine areas and the Bouguer anomaly, computed with a reduction density of 2670 kg/m3, over land masses. Also shown are the coastline, major tectonic features of the northern Juan de Fuca plate system and the location of the VISP II profile. The long wavelength features of the combined free-air/Bouguer gravity field exhibit a general northwest-southeast trend. From southwest to northeast a simple alternating anomaly pattern is evident. Oceanic crust seaward of the Cascadia deformation front is characterized by a slightly positive anomaly. The gravity field is more negative over the deformation front and continental Chapter 4. GRAVITY MODELLING 62 shelf, positive over Vancouver Island, and large negative values represent the area of the Coast Mountain Belt. The VISP II refraction profile (figure 4.1) is aligned roughly parallel to the strike of the gravity field over the northern part of Juan de Fuca plate. The cause of some of the small scale features in the free-air anomaly is revealed from an ana-lysis of bathymetry. Figure 4.2 shows a sub-region of the gravity grid in the vicinity of VISP II draped over a 3-D perspective of the corresponding bathymetry. A number of short wavelength free-air anomalies over oceanic crust correlate to prominent bathymetric features. These in-clude positive anomalies over the Heck seamount chain and the Paul Revere and Winona ridges. Bathymetric relief across the Nootka fault zone is small and does not produce a distinguishable gravity anomaly. In contrast to the Pacific plate, seamount chains are conspicuously absent from bathymetry over the Juan de Fuca plate. Thick sediments mask the basement topography but Calvert et al. [1990] imaged a ~ 10 km wide seamount at 47.3 ° N and 128.1 ° W from a seismic reflection line on the Juan.de Fuca plate. A gravity anomaly is not detectable from this feature because the density contrast between the basaltic basement and the enveloping sediments is small. 4.1.2 Previous gravity studies Numerous models have been constructed to explain the gravity field over the Cascadia subduc-tion zone (e.g. Riddihough [1979], Clowes et al. [1997]). Of those that extend over the relat-ively quiet part of the gravity field west of the Cascadia deformation front, very simple structure is used to model the Juan de Fuca plate. Clowes et al. [1997] present 2.5-D models from a series of lines that transect the north-ern Juan de Fuca plate at an orientation normal to the strike of the Cascadia subduction zone. Oceanic crust is modelled using two layers to approximate a 2-4 km thick sequence of sediments and one 6-7 km thick layer to replicate the igneous crust. The gravity low over and to the east Chapter 4. GRAVITY MODELLING 63 Free-Air Gravity (mgal) -59 -43 -35 -30 -27 -23 -20 -17 -14 -12 -10 -9 -8 -6 -5 -4 -2 0 5 Figure 4.2: 3-D perspective of offshore British Columbia: free-air gravity draped over bathy-metry. Major tectonic features are shown including the approximate location of the Cascadia deformation front. Prominent bathymetric highs west of this boundary coincide with positive gravity anomalies. of the deformation front is modelled by a wedge of accretionary sediments and the increasing depth of the oceanic lithosphere as it subducts beneath North America. 4.1.3 Resolution of the grid data The gravity measurements used to produce the offshore grid data are contributed from various surveys carried out aboard marine research vessels. Ship-board gravity measurements are typic-ally made along a traverse at uniform intervals. The ship tracks themselves may have a variable orientation and spatial distribution such that some regions are better sampled than others (e.g. figure 4.3). Each gravity observation made along a ship track has a location (X,Y) and a gravity value (Z). To enable easy manipulation of the gravity data, XYZ data is best represented in the form of a regularly spaced grid. The process of gridding takes XYZ data and interpolates the value of Z at the nodes of a grid. Grid data from the National Gravity Database are interpolated using a minimum curvature Chapter 4. GRAVITY MODELLING 64 method (e.g. Briggs [1974]) such that the smoothest possible surface is fit to the given data values. The resolution of the grid data is dependent on both the spacing between grid nodes and the density of the original gravity observations. The grid data are subject to spatial aliasing where the node interval is large relative to the observations. In this instance significant short wavelength anomalies are smoothed over. Conversely, if gravity observations are sparse, short wavelength features will not be sampled at all and a finely spaced grid may appear to have a higher resolution than in actuality. The accuracy of the contoured grid data shown in figure 4.3 is higher along ship track lines than between them. The majority of the ship-board gravimetry traverses over the northern Juan de Fuca Plate cross the VISP II profile at an almost perpendicular orientation and are spaced ~ 10 km apart. The fine structure of the gravity field along the northern part of the profile is controlled by a sub-parallel ship track, on which gravity observations are spaced a maximum 2 km apart, that extends no more than 10 km to the west between OBS-2 and OBS-4. Thus, it appears that the region of anomalous velocity structure, between OBS-2 and OBS-4, is covered by a sufficient density of observations and that a 2 km grid spacing is sufficiently small to sample the gravity expression expected from a density anomaly at the same location. 4.2 Development of the gravity model A density model constructed in two dimensions is desired to complement the existing 2-D ve-locity model. The broad regional gradient of the free-air gravity anomaly conveniently strikes subparallel to the VISP II refraction profile (figure 4.3). Thus, the regional negative gravity gra-dient attributed to the subducting lithosphere and accretionary wedge sediments to the east will not affect the relative pattern of the gravity field over the VISP II profile. Contributions to the profile free-air anomaly from out-of-plane structure that can not be approximated from 2.5-D modelling are considered negligible. Chapters GRAVITY MODELLING 65 -''29° -128° -127" -126° Figure 4.3: Gravity station locations and the contoured free-air gravity anomaly within the sur-vey area. Ship tracks that traverse the VISP II profile are shown as dotted blue lines. Gravity measurements were made a maximum of 2 km apart along each ship track. The contour interval for the free-air gravity anomaly is 5 mgal. Contour depths to the crust-mantle boundary estab-lished from the 2-D VISP II velocity model and a 1-D model of refraction line EX2R [ Au and Clowes, 1982] are shown as dashed red lines. Chapters GRAVITY MODELLING 66 4.2.1 Data procedures and modelling method The gravity data are manipulated using GEOSOFT commercial processing software. Profile gravity values, extending 100 km beyond each end of the VISP II profile, were interpolated from the grid data using GEOSOFT. Coordinates for the spatial data are reprojected from geographic latitude and longitude coordinates to x and y coordinates using the WGS84 reference ellipsoid. The observed gravity values are sampled at ~ 2 km intervals to be used as input for the SAKI 2.5-D forward modelling program [ Webring, 1985]. SAKI requires a density model parameterization similar to that of RAYINVR (described in section 3.1.1). 2-D polygons with a variable number of nodes are used to model homogeneous mass blocks of finite or infinite length that extend into and out of the 2-D profile plane. The theoretical basis is discussed more thoroughly in Appendix A. The gravity effect of the given density distribution is calculated relative to a specified reduction density. The model is updated until a satisfactory fit to observed gravity values is achieved. A reduction density equal to the density of the upper mantle (3300 kg/m3) is subtracted from the layer densities to calculate the synthetic gravity anomaly. This value was chosen so that the gravity effect of assigning fixed dimensions to layer 3B perpendicular to the profile plane could be calculated relative to the upper mantle density. A floating datum is used in the forward modelling program sq that relative changes in the calculated field can be directly compared to the observed anomaly. Edge effects are suppressed by extending the model length an additional 200 km at each end of the gravity profile. 4.2.2 Rock densities The Nafe-Drake curve summarizes the general relationship between compressional wave veloc-ity and density of commonly occurring sedimentary and igneous rocks [ Ludwig et al., 1970]. The curve is compiled from laboratory measurements from a variety of rock types and provides Chapter 4. GRAVITY MODELLING 67 an approximation of velocity/density relationships for typical marine geology. Barton [1986] cautions against directly converting from velocity to density using the relationship because there is a considerable range of densities for any given velocity. Velocity/density relationships are most accurate where the rock composition can be inferred such as occurs for surveys of oceanic crust. Densities used to describe the 2 sedimentary layers are extracted from Hamilton's [1978] ve-locity/density relationships for turbidites and mudstone seafloor sediments. Carlson and Raskin [1984] determined a velocity/density relationship from measurements of igneous rocks that com-pose typical oceanic crust. Density values were determined from this relationship using the average P-wave velocity of seismic layers 2A, 2B, 2C, 3A and 3B established from the veloc-ity model. A value of 3300 kg/m3 is assigned to the upper mantle consistent with Clowes et al. 's [1997] results from gravity modelling of crustal structure across the Cascadia subduction zone. 4.3 Initial model Much like the velocity model, the relative simplicity of oceanic crustal structure allows mass blocks to be approximated as 9 subhorizontal layers of infinite extent perpendicular to the pro-file plane. The geometry of the favoured VISP II velocity model was replicated to produce an initial density distribution model (figure 4.4). Layer thicknesses at 145 km model distance are defined from Au and Clowes' [1982] interpretation of refraction line EX2R across the Nootka fault zone. Otherwise, layer thicknesses are manipulated to fit the gravity values measured be-yond the extent of the VISP II profile. The structure of the crust beyond each end of the VISP II profile exhibits little variation from that modelled beneath OBS-1, which represents the typ-ical structure of normal oceanic crust. The required transition at the northwestern end of the profile, from abnormal crust that is 10 km thick back to crust with a normal thickness (7 km), Chapter 4. GRAVITY MODELLING 68 takes place northewest of the VISP II profile between 110 km and 145 km model distance. The change coincides with the approximate position of the Nootka fault zone and implies that the fault zone juxtaposes normal oceanic crust of the Explorer Plate against anomalously thick crust of the northern Juan de Fuca Plate. Bathymetric relief was sampled every 5 km from grid data to define the geometry of the seafloor. The depth of the model is limited to 15 km. Table 4.1 lists the parameters used in the construction of the initial density distribution model (figure 4.4). Layer AZ (km) p (kg/m3) cr Sea water 2.53 - 2.60 1030 0.5 Sediments (1A) 0.41-0.63 1900 0.46 Sediments (IB) 0.36-0.91 2200 0.45 2A 0.27 - 0.99 2500 0.30 2B 0.30-1.30 2750 0.25 2C 0.25 -1.10 2850 0.25 3A 2.30 - 3.20 2900 0.25 3B 3.30-4.30 2970 -Upper mantle - 3300 -Table 4.1: Vertical layer thickness, density and Poisson's ratio parameters for the initial gravity model (figure 4.4) established from the favoured velocity model (see section 3.4.1). Poisson's ratios (cr) are indeterminate for Seismic Layer 3B and the upper mantle. The model is displayed in figure 4.4 and shows the prominent ~ 25 mgal gravity low not evident in the observed data. The observed and modelled data are clearly discordant and appear to disqualify the velocity/density relationships used to derive the density model. Some of the discrepancy can be attributed to velocity heterogeneity in each layer that is not accounted for when a uniform density is applied. Adding this complexity to the density model may change the short wavelength amplitudes but will not significantly reduce the broad low that extends for ~ 100 km. Chapter 4. GRAVITY MODELLING 69 — Calculated kilometres -200 -100 0 100 200 300 1 ' ' 1 ' 1 ' 1 1 1 1 ' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1030 1900 2200 2500 2750 2850 2900 2970 3300 Density (kg/m 1) Figure 4.4: The initial gravity model calculated from the parameterization of the best fit RAY-INVR velocity model (figure 3.10). The dashed white line indicates where the crust-mantle boundary is inferred from the VISP II velocity model and the 1-D velocity model of refrac-tion line EX2R [ Au and Clowes, 1982]. Also shown are OBS positions (inverted triangles) and the location of the Nootka fault zone (cross hatching), with the direction of motion indicated by sense symbols. The dotted white line bounds the region of the density model constrained by ray tracing through the favoured velocity model. The gravity responses relative to observed values are shown above. The dashed line represents the calculated response when lower crust at greater than 11 km depth is given fixed dimensions of 20 km into and 20 km out of the pro-file plane. The dotted line shows the response when the lower crustal protrusion is completely removed. Vertical exaggeration is 12:1. 4.3.1 Source of the calculated anomaly The velocity-density relationships used to derive the density model have done a poor job recon-ciling the seismic data with the gravity data. The density model requires a significant increase in mass beneath the northwestern end of the VISP II profile for the calculated response to match the observed data. Such an increase in mass would require that unusual material exists within the crust or upper mantle that has unexpected physical properties. The limited S-wave veloc-ity data used to derive Poisson's ratios for each of the seismic layers (section 3.4.1, table 4.1) Chapters GRAVITY MODELLING 70 indicates that the mid to upper crust exhibits elastic properties typical of basalts and cumulate gabbros that exist in normal oceanic crust (cr = 0.25-0.3). Elastic properties for Layer 3B and the upper mantle could not be determined and without evidence to the contrary suggest that this part of the model is the most likely place for anomalous material to exist. Parameters defining layer 3B are manipulated to test how the fit of the calculated response to the observed data can be improved. Figure 4.4 shows the gravity response when the crustal root is completely removed and replaced with a density equal to that used to model the upper mantle (3300 kg/m3). The good agreement between calculated and observed values indicates that the mass deficient crustal root has the greatest contribution to the calculated gravity low. The misfit of the calculated anomaly is marginally reduced by limiting the areal extent of the over-thickened crust that protrudes into the upper mantle. The dashed line exhibited in figure 4.4 represents the calculated response when the lower crustal protrusion is limited to 20 km into and 20 km out of the profile plane. Smaller dimensions would reduce the feature to an unrealistic. plate-like shape with a strike coincident to that of the profile. Based on these results, anomalous material that has a density more akin to that of the upper mantle resides in the lower crust. 4.4 Serpentinization of peridotite mantle Results from velocity and gravity modelling require that in the interpreted region of the lower crust material characterized by velocities in the range of 7.0 r 7.2 km/s and densities of > 3000 kg/m3 exists. The alteration of peridotite by serpentinization can produce the desired phys-ical properties and has been used to satisfy apparently discordant gravity and seismic datasets elsewhere. For example, Canales et al. [2000], in a seismic refraction and gravity study across a nontransform offset of the Mid-Atlantic Ridge, invoke the presence of serpentinized peridotites, between 3.8 km and 6.2 km depth below the seafloor, to explain modelled velocities of 7.0-7.4 km/s coincident with densities of 3000-3250 kg/m3. .' Chapter 4. GRAVITY MODELLING 71 Serpentinites are rocks composed of the serpentine-group minerals; principally lizardite, chrysotile and antigorite which have the approximate composition Mg3[Si205](OH)4. Serpent-ine minerals readily replace olivine and orthopyroxene in ultramafic rocks. Serpentinized ultra-mafic rocks are frequently recovered from ocean basins on fault escarpments within spreading ridges and along transform offsets between spreading segments [ Fryer, 2002]. The most sig-nificant change in composition induced by the serpentinization of peridotite is the addition of Up to 15 weight % H 2 0 [ O'Hartley, 1996]. The addition of water can greatly reduce a peridotite's seismic velocity and density. Christensen [1972; 1966] analyzed a suite of peridotites serpentinized to varying degrees to determine relationships between compressional- and shear-wave velocities arid density at vary-ing pressures. Figure 4.5 illustrates the results obtained at 4 kbar (equivalent to ~ 13 km depth) and shows how P-wave velocities decrease in a linear fashion with the increasing extent of ser-pentinization. Note that peridotites consisting of ~ 25% serpentine exhibit a P-wave veloc-ity of ~ 7.2 km/s and a density of ~ 3100 kg/m3. This velocity is equivalent to that used to model the base of seismic layer 3B in the region of overthickened crust. This density is sub-stantially higher than that predicted by velocity/density relationships for cumulate gabbros and from Christensen's [1996] measurements of gabbroic rocks at pressures equivalent to 13 km depth. Thus, it is possible that a partially serpentinized body of peridotite mantle occupies the region of over-thickened crust interpreted from the velocity model. 4.4.1 Incorporating a serpentinized body For the purposes of minimizing structural change, layer 3B is first subdivided into 2 bodies. A serpentinized body within seismic layer 3B is modelled where velocities of 7.2 ± 0.05 km/s were used to model the region of interpreted thick crust (figure 3.10). This region is reassigned a density of 3150 kg/m3 and extends between 40 and 145 km model distance (figure 4.6). The Chapter 4. GRAVITY MODELLING 72 Volume % Serpentine 100 75 50 25 0 8.51 1 1 1 1 1 1 1 1 2500 2600 2700 2800 2900 3000 3100 3200 3300 Density (kg/m3) Figure 4.5: P-wave velocities at 4 kbar versus density and vol % serpentine for peridotites and serpentinized peridotites. The solid line represents the least squares fit and the dotted lines rep-resent the standard deviation of the residuals from the measurements of Christensen [1966]. Also shown is the mean value (solid circle) and standard deviation (dotted ellipse) for a suite of gabbro samples from Christensen [1996]. Ap indicates the contrast in density at 7.2 km/s for gabbros and peridotites containing 25 % serpentine. updated density lies within 1 standard deviation of the value calculated for peridotites contain-ing 25% serpentine by volume (figure 4.5). The fit between the calculated gravity response from the model with just 1 serpentinized body and the observed data is greatly improved from dif-ferences as great as 25 mgal in figure 4.4 to less than 10 mgal. The root-mean-square (RMS) of the residuals is reduced from 9.16 mgal to 3.03 mgal. Further improvement is made by ex-tending the serpentinization to include velocities of 7.1 ± 0.05 km/s within seismic layer 3B. Following the relationship of figure 4.5, a velocity of 7.1 km/s matches a perdidotite containing ~ 30 % serpentine by volume which corresponds to a density of ~ 3050 kg/m3. This body with a new density of 3050 kg/m3 is incorporated into the model. Figure 4.6 shows how the Chapter 4. GRAVITY MODELLING 73 long wavelength component (> 50 km) of the calculated anomaly from both serpentinized bod-ies is close to an optimal fit to the observed data. The residual RMS is reduced from 3.03 mgal to 2.47 mgal. The remaining misfit results from the shorter wavelength (< 50 km) variations of the calculated data and may be reduced through modifications to the density structure of the shallower part of the crust. 10 i 0 o ca 00 ,§-10 1 -20 O -30 1030 1900 2200 2500 2750 2850 2900 2970 3050 3150 3300 Density (kg/m3) Figure 4.6: The density model and calculated gravity response modified to include a region of serpentinized peridotite modelled as 2 bodies (horizontal and vertical stripe pattern). The dashed black line is the response calculated from the modification of the model in figure 4.4 by the inclusion of body S1 only. The solid black line is the response calculated from the inclusion of both S1 and S2 bodies. Further explanations are contained in the caption of figure 4.4. The model suggests that the two bodies of partially serpentinized peridotite replace seismic layer 3B of the Juan de Fuca plate only, starting at 40 km before tapering in the region of the Nootka fault zone ( 140 km). With the removal of part of seismic layer 3B, oceanic crust at the northwestern tip of the Juan de Fuca plate now appears to thin approaching the Nootka fault Chapter 4. GRAVITY MODELLING 74 zone. 4.4.2 Reducing the misfit of short wavelength variations Some of the misfit illustrated in figure 4.6 may be attributed to model parameter uncertainties. These include the uncertainty of the position of the structural nodes which define a layer's geom-etry and the uncertainty of velocity/density relationships used to assign density values. In addi-tion, as mentioned in section 4.3, lateral heterogeneity within each seismic layer is essentially smoothed out when approximated with an average density. The nodes defining the boundaries of 1 A, IB, 2A, 2B and 2C were adjusted until the calcu-lated response produced a reasonable match to the observations. Nodes defining the boundary between the sediments and the top of the igneous basement were not moved as their position is defined from the CSP profiling (section 2.1). Figure 4.7 illustrates the changes made to the model by this process. The RMS of the residuals for the model is reduced from 2.47 mgal to 1.87 mgal which, compared to the defined measurement uncertainty of ± 2 mgal, represents a good fit to the observations. The final model is presented in figure 4.8. The same changes to the density model layer boundaries are applied to the velocity model to evaluate the new misfit to the picked first arrival phases. The statistics of the first arrival misfit from ray tracing through the velocity model before and after the adaptation are displayed in table 4.2. Before After TRMS (S) n 0.137 1.604 189 0.143 2.229 189 Table 4.2: TRMS and xl residual traveltimes of first arrivals for the 2-D velocity model before and after reconciliation with the density model illustrated in figure 4.8. Chapter 4. GRAVITY MODELLING 75 o 20 Model Distance (km) 40 60 80 100 120 1A 2A IB -300 -100 -50 50 100 300 Density Change (kg/m3) Figure 4.7: Density changes resulting from adjustment of boundaries defining layers 1A, IB, 2A, 2B, 2C and 3A. The region of the model not constrained by the velocity model ray coverage is masked out. The modifications have increased the misfit but not to an unsatisfactory extent. Of all the measurements used to integrate the 2 models, those used to define the velocity-density relation-ships carry the greatest uncertainty. For example, Christensen's [1996] measurements of 252 basalt specimens at 2 kbar (equivalent to 6-7 km depth) exhibit standard deviations for P-wave velocity and density of 0.546 km/s and 139 kg/m3 respectively. Perturbing just layer 2A of the best fit velocity model by 0.5 km/s yields a calculated xl misfit of 2.218. Therefore, the differ-ence in misfit between the best fit velocity model and the velocity model modified to reconcile with the density model is within acceptable error. 4.5 Summary Gravity data over the VISP II profile do not reveal a gravity low over the region of thick crust postulated from the velocity model. The absence of a correlative anomaly suggests that there exists material in the lower crust with anomalous physical properties: P-wave velocities of the range 7.0 - 7.2 km/s but densities of at least 3000 kg/m3. The seismic refraction and gravity data Chapter 4. GRAVITY MODELLING 76 10 r 1030 1900 2200 2500 2750 2850 2900 2970 3050 3150 3300 Density (kg/m3) Figure 4.8: The final density model in which layers in the mid to upper crust are modified such that the calculated gravity response best matches the observed data. Further explanations are contained in the caption of figure 4.4. are satisfied by incorporating a broad zone of partially serpentinized peridotite into the region previously interpreted as thick crust. The model confines the region of serpentinization to the Juan de Fuca plate and now indicates that northern Juan de Fuca crust thins approaching the Nootka fault zone in the vicinity of VISP II. Chapter 5 DISCUSSION AND CONCLUSIONS The low velocities observed from the favoured velocity model would suggest anomalously thick lower crust predominantly composed of gabbroic material. However, this interpretation is not consistent with the gravity field. The densities for the same model geometry required to match the gravity field indicate the presence of serpentinized peridotites beneath anomalously thin crust. Other seismic experiments within the vicinity of VISP II may add credence to, and help define the regional extent of, the postulated zone of serpentinization. In addition, the geophys-ical constraints aside, how any inferred zone of serpentinized peridotite developed requires fur-ther consideration. 5.1 Supporting Evidence Two independent seismic experiments in the immediate vicinity of the VISP II profile provide complementary evidence for unusual velocity structure and crustal characteristics within the northern Juan de Fuca plate (figure 5.1). 5.1.1 The 1977 Nootka fault zone wide-angle experiment results The unusual travel time delays observed from the VISP II receiver gathers, which yield the model of upper mantle serpentinization (figure 4.8), are replicated by a refraction line from an earlier experiment over the Nootka fault zone [ Au and Clowes, 1982]. Refraction line EX3 is oriented parallel and to the south of the Nootka fault zone (figure 5.1). This line intersects two lines, EX1 and EX2, that traverse the fault zone to the north. 77 Chapter 5. DISCUSSION AND CONCLUSIONS 7 8 130°W 129°W 128'W 127'W 126*W 125'W 124°W 130"W 129"W 128°W 127*W 126°W 125°W 124"W Figure 5:1: Proximal seismic experiments over the northern Juan de Fuca plate. The reflection lines are indicated by bold lines and.refraction lines as dotted lines with solid circles defining each OBS position. Arrows indicate the shooting direction for profiles EX1, EX2 and EX3. The shaded zone represents the inferred regional extent of serpentinization established from the seismic experiments. Variations in interpreted crustal thickness from 1-D models for forward and reversed sec-tions required the implementation of a 2-D technique to model the laterally varying structure beneath line EX3. A simple 2-D model was constructed from a combination of three 1-D mod-els with smoothing applied to account for the lateral variation between each. By using an early ray tracing technique, synthetic seismograms were produced that generally approximated the observed data. Figure 5.2 is a reconstruction of this model and shows the dramatic change in the interpreted depth to the crust-mantle boundary, from 9 km near the Juan de Fuca Ridge to 13.6 km at the northeastern end of the profile and at the crossing with the VISP II profile. These Chapter 5. DISCUSSION AND CONCLUSIONS 79 depths equate to igneous crustal thicknesses of 6.4 and 11 km at the respective ends of the pro-file. From their model the onset of crustal thickening is inferred to occur at a position along profile equivalent to 48.65 ° N and 128.25 ° W. This position corresponds to approximately 25 km offset in a southwest direction normal to the V I S P II profile. SW CMB EX2 VISP n Crossing I EX3R 20 30 40 Model Distance (km) NE EX3 CMB 50 60 (km/s) 1.48 1.50 4.00 5.00 6.00 6.70 7.20 7.60 8.00 8.50 Figure5.2: 2-D velocity model of line EX3 (figure5. l)fromAw and Clowes [1982]. The model is a piecewise construction of 3 1-D models: EX2, EX3 and EX3R (the reversal of EX3). Dotted lines define the proportions used from each 1-D model and the zones of smoothing between to account for lateral variation. Dashed white line indicates the crust-mantle boundary. Also shown is the crossing with the V I S P II profile. The point of intersection corresponds to 120 km model distance on V I S P II. The lateral velocity gradient in the upper mantle is attributed to mantle anisotropy as lines EX3 and EX3R are oriented closer to the ridge spreading direction than EX2. 2-D models were not constructed for lines EX1 and EX2. Forward and reversed 1-D models for EX2 suggest a crust of relatively uniform, normal thickness. A model produced for line EX 1 was constructed solely from a forward section and indicates a crust with a thickness of ~ 8 km. Chapter 5. DISCUSSION AND CONCLUSIONS 80 This is an intermediate value, between the thick crust of the Juan de Fuca and the relatively nor-mal crustal thickness exhibited on the Explorer plate side of the Nootka fault zone, and results from averaging the laterally varying crustal thickness to produce a 1-D model. Interpretation of the models confines thick crust to the Juan de Fuca plate only. 5.1.2 The 1985 Juan de Fuca plate reflection experiment results In 1985, several multichannel seismic reflection (MCS) lines were recorded off the west coast of British Columbia for the purpose of investigating the lithospheric structure over the continental margin and adjacent Juan de Fuca plate. Two of these lines (85-07 and 85-09), perpendicular to and southwest of the VISP II profile (figure 5.1), are presented by Hasselgren and Clowes [1995]. Results from an additional two lines (85-01 and 85-02), perpendicular to and northeast of VISP II are described by Hyndman et al. [1990]. Hasselgren and Clowes [1995] observe continuous, coherent reflection Moho at 6-6.5 s two-way travel time (TWT) on both lines 85-07 and 85-09. By assuming an average velocity through the crust of 6000 m/s, the depth to the top of the crust-mantle transition (CMT) reflector varies as little as 600 m for most of the length of each profile. However, at the eastern end of line 85-07 (shot points 1-600 in figure 5.3).CMT reflectors disappear and are replaced by an anomalous zone extending from 5 s to at least 7 s TWT containing dipping, discontinuous reflectors. In-terval velocities derived from stacking velocities for the deep reflectors revealed velocities of the range 6.6-7.5 km/s, typical of lower crustal gabbroic material. The authors suggest that the complex reflectivity results from the interlayering of materials with contrasting mafic content and ascribe this feature to an abrupt change in crustal thickness coincident with the trace of a propagating rift pseudofault. The same pattern of reflectivity was also observed at the crossing with another pseudofault on line 85-09. The onset of anomalous reflectivity (shot points 1-600 in figure 5.3) coincides with a change Chapter 5. DISCUSSION AND CONCLUSIONS 81 VISP II Crossing V.E. -.2:1 Figure 5.3: Coherency filtered section (top) and interpreted line drawing (bottom) of reflec-tion line 85-07. Note the onset of dipping reflectors at times below the interpreted crust-mantle transition (CMT) over a zone extending east of shot point 600. PF7 refers to the crossing of a pseudofault as interpreted by Wilson et al. [ 1984]. Also shown is the crossing of line 85-07 with the VISP II profile. The point of intersection corresponds to 65 km model distance on VISP II. S is sediments; D is disrupted sediments, B is top of igneous crust; W is water bottom multiple. Hasselgren and Clowes [1995] interpreted IC as intracrustal reflectors and U as an underplated zone. Reproduced from Hasselgren and Clowes [1995]. from relatively smooth to rough topography of the basement surface. The roughness of base-. ment topography can be used as a proxy for the rate of sea-floor spreading. This is exempli-fied by Pacific-type and Atlantic-type topography: basement straddling the fast-spreading East Pacific Rise is considered relatively smooth when compared to that of the slow-spreading Mid-Atlantic Ridge [ Jones, 1999]. Slower spreading rates are associated with decreased melt produc-tion to the accreting plate and consequently thinner crust [ White et al., 1992]. Therefore, the relatively rough topography exhibited east of shot point 600 on line 85-07 should imply thin-ning of the crust rather than an increase in thickness. The deep lenticular reflectivity may sug-gest interlayering of ultramafic material exhibiting contrasting degrees of serpentinization. The absence of similar reflectivity at the eastern end of line 85-09 limits the southern extent of the low velocity zone to somewhere between lines 85-07 and 85-09. Chapter 5. DISCUSSION AND CONCLUSIONS 82 Seismic sections for reflection lines 85-01 and 85-02 east of VISP II (figure 5.1) from Hyn-dman et al. [1990] show little evidence for anomalous reflectivity which may be interpreted as serpentinized upper mantle. However, reflectivity associated with the CMT is also difficult to resolve on line 85-02 as it traverses just 30 km of oceanic crust seaward of the Cascadia de-formation front. Therefore, the eastern-most boundary of the serpentinized zone is tentatively positioned just a few kilometres east of the VISP II profile (figure 5.1). 5.1.3 Comparison with free-air gravity map Unlike seismic line VISP II, profiles EX3 and 85-07 are oriented perpendicular to the strucutural grain of the Cascadia subduction zone. Consequently, at these orientations, variations in the gravity field resulting from local structure are overprinted with the regional variation of larger scale structure attributable to the adjacent subduction zone (figure 5.4). For example, accre-tionary wedge sediments east of the deformation front and the deepening of the lithosphere as it subducts beneath North America contribute to the large negative gravity gradient in a north-easterly direction approaching the Cascadia subduction zone. The gravitational effects of the Cascadia subduction zone make it difficult to discriminate a gravity low expected from the onset of thick crust along seismic lines EX3 and 85-07. In or-der to make this distinction it is necessary to predict the local gravitational effect of thick crust at these locations to compare with the regional gravity field. Figure 5.5 shows free-air gravity measured along a profile coincident with seismic refraction line EX3. By using the geometry of the interpreted velocity model (figure 5.2), a simple 2-D synthetic gravity anomaly is gener-ated for the mass deficiency expected where lower crust replaces upper mantle. The synthetic anomaly has a large -40 mgal amplitude that does not correlate with the observed data (figure 5.5). Furthermore, any correlation is further reduced when long wavelength components (> 250 Chapter 5. DISCUSSION AND CONCLUSIONS 83 -128° -127° -126° Figure 5.4: Free-air gravity map overprinted with observed low velocity regions. Seismic lines are displayed in red where arrows indicate forward shooting directions. Portions of each seismic line inferred to coincide with regions of serpentinized upper mantle are shown with thick green bars. The Cascadia deformation front is shown in blue. Gravity contour interval is 5 mgal. km) of the free-air anomaly, attributed to larger scale deep structure within the Cascadia sub-duction zone, are filtered out. The same argument can be used to refute the presence of thick crust beneath the eastern end of seismic line 85-07. From analysis of the gravity data and their proximity to the VISP II profile, the low velocity zones observed from seismic lines EX3 and 85-07 are interpreted to form a broad zone of partially serpentinized peridotite mantle (figure 5.1). Chapter 5. DISCUSSION AND CONCLUSIONS 84 sw 5 \ -Q 10 15 E X 3 profile distance (km) E X 3 velocity jnoddlengh^ Figure 5.5: Synthetic gravity anomaly (dotted line) and observed free air gravity anomaly for profile EX3 (solid line). The synthetic anomaly is generated using a simple approximation of a 2-D mass-deficient crustal root with density contrast of -300 kg/m3. The geometry of the root is taken from the 2-D velocity model of EX3 (figure 5.2) and is assumed to be limited to the extent of the EX3 profile. The dashed line represents the free-air gravity anomaly following the application of a high pass filter with a 250 km cutoff wavelength (wavelength components greater than 250 km are removed). 5.2 Serpentinization of Peridotite in Ocean Basins Prior to Christensen's [1966; 1972] measurements of both compressional and shear velocities for partially serpentinized peridotites, some authors suggested that P-wave velocities for seis-mic layer 3, observed from marine refraction studies, indicated a composition of partially ser-pentinized peridotite for typical oceanic lower crust (e.g., Hess [1962]). Subsequent global measurements of P- and S-wave velocities discounted this hypothesis because the observed Poisson's ratios (0.25-0.28) were incompatible with those measured from serpentinite samples Chapter 5. DISCUSSION AND CONCLUSIONS 85 in the laboratory (> 0.3). These results suggested that serpentinization does not always accom-pany the formation of new oceanic crust and, as an alteration event, requires particular envir-onmental conditions for the reaction to arise. 5.2.1 Thermal regime of oceanic lithosphere Serpentinization is an exothermic reaction that only proceeds at temperatures below 500 ° C within the uppermost mantle [ O'Hanley, 1996]. As such, ultramafic rocks crystallizing at a spreading centre require a sufficient period of cooling for serpentinization to procede. Beneath a spreading ridge, the temperature profile in the asthenospheric wedge is adiabatic and equi-valent to ~ 13000 [ Turcotte and Schubert, 1982]. As the newly crystallized part of the plate moves away from the ridge it undergoes conductive cooling. This cooling can be approximated by a simple one-dimensional model of a uniform half-space with a fixed surface temperature. Figure 5.6 illustrates geothermal profiles generated by this method for lithosphere of different ages. Derivation of the model and its limitations are discussed in Appendix B. Figure 5.6 provides an estimate for the depth to which the serpentinization reaction can pro-ceed. The zone of serpentinized peridotite observed from the final VISP II model (figure 4.8) ex-tends from approximate depths of 4 km to 9 km beneath the surface of the igneous oceanic crust. For serpentine minerals to remain stable at these depths it is required, from the simple model of conductively cooling oceanic plates, that the hydrated lithosphere have undergone at least 5 million years of cooling before the onset of serpentinization. However, the model of conduct-ive cooling neglects the influence of hydrothermal currents which significantly increase the rate of lithospheric cooling at distances close to the spreading ridge (see Appendix B). Therefore, the upper mantle may obtain sub-500 ° C temperatures within only a short distance of formation beneath the ridge crest. Chapter 5. DISCUSSION AND CONCLUSIONS 86 1400 Temperature (°C) Figure 5.6: Geotherms for conductively cooling oceanic lithosphere with ages 1, 2, 5 and 10 million years. The 500 ° C isotherm represents the upper limit for the stability of serpentine minerals. The solid grey bar defines the depths over which serpentinization is inferred to occur from the VISP II model (figure 4.8). 5.2.2 Serpentinization models Hess [1962] originally proposed that serpentinization of peridotite took place as upwelling up-per mantle cooled below 5000 C and was hydrated by its own juvenile water. Chemical studies of serpentinite indicate that serpentinization of oceanic mantle requires large volumes of water much greater than that of the original rock [ Francis, 1981]. Francis [1981] argued that a more likely mechanism of hydration takes place by hydrothermal circulation of sea water into the up-per mantle beneath recently created oceanic crust. A necessary condition for serpentinization of upper mantle rocks to occur is a high flow rate of sea water into the mantle. The maintain-ance of a high flow rate requires permeable pathways, like faults or fractures, for sea water to penetrate into the mantle and vigorous hydrothermal circulation such as occurs in the vicinity of Chapter 5. DISCUSSION AND CONCLUSIONS 87 a spreading axis. Vigorous hydrothermal circulation also allows upper mantle temperatures to be lowered to within the stability field of serpentine minerals at distances close to the spreading ridge axis. Macdonald and Fyfe [1985] show that at temperatures of ~ 300 ° C the complete serpentinization of a 1 km thick layer of peridotite can take place in as little as 1 million years. Fryer [2002] and O 'Hanley [ 1996] describe a process of oceanic crustal formation and modi-fication that may explain why serpentinized peridotites are frequently recovered from fault es-carpments within slow-spreading (1 to 5 cm/yr) ridge axes such as the Mid-Atlantic Ridge and not from fast-spreading (greater than 9 cm) ridge axes like the East Pacific Rise. Unlike fast-spreading ridges, the supply of magma to slow-spreading ridges is episodic, suggesting that the ridge will endure amagmatic periods. The crust will thin without the addition of new magma to the spreading plates and does so by listric normal faulting. Salisbury and Keen [1993] show clear images of listric normal faults from deep seismic reflection profiles off the coast of Nova Scotia and suggest that they are a common feature of crust formed at slow-spreading ridges. In addition, the absent heat source beneath a ridge during a period of amagmatism allows iso-therms to deepen such that temperatures within the mantle are conducive to serpentinization. The normal faults provide the necessary pathways for sea water penetration into the upper man-tle. Francis [1981] suggests that serpentinization at a normal fault is more likely to take place beneath the upthrown block as this is where isotherms are displaced downward relative to the footwall. The volume expansion that accompanies the serpentinization process increases the bouyancy of the upthrown block, providing additional uplift and lengthens the expression of the fault scarp. Segments of the Juan de Fuca ridge system spread at rates of anywhere between 4 and 6 cm/yr [ Riddihough and Hyndman, 1991]. The average rate of spreading at the Juan de Fuca Ridge has been relatively consistent over the last 6 Ma [ Riddihough and Hyndman, 1991]. Therefore, the spreading rate at the time of the formation of crust on the northern Juan de Fuca Plate falls into the high end of the slow-spreading ridge category. Why then should this region Chapter 5. DISCUSSION AND CONCLUSIONS 88 of serpentinization, illustrated in figure 5.1, be exclusively confined to the northernmost part of the Juan de Fuca plate? An understanding may be achieved by reconstructing the history of crustal formation in this region. 5.3 Timing of Serpentinization The hydration of the upper mantle beneath the northern Juan de Fuca region may be related to the timing of crustal formation at the Juan de Fuca Ridge. Hasselgren et al. [ 1992] proposed that anomalous reflectivity observed on seismic reflection lines 85-07 and 85-09 coincided with the initiation of 2 propagating rift pseudofaults. Further insight into the crustal formation process is gained by the construction of a crustal age map for the northern Juan de Fuca plate. 5.3.1 Crustal Age Data Magnetics data, gridded at 2 km, were obtained from the Geophysical Data Centre of the Geo-logical Survey of Canada. The magnetic anomaly pattern is displayed by polarity in figure 5.7. Age contours and the position of pseudofault traces PF5 and PF7 are determined following the interpretation of Wilson [1993]. Ages for the Explorer plate are difficult to determine, but it is known that the Nootka fault zone juxtaposes younger crust of the Juan de Fuca plate against plder crust on the Explorer side because the fault zone has been accommodating left-lateral mo-tion since at least 4 Ma [ Riddihough, 1984]. Chapter 5. DISCUSSION AND CONCLUSIONS 89 129°W 128'W 127'W 129"W 128'W 127'W Figure 5.7: Magnetic anomaly map for the northern Juan de Fuca plate. Positive anomalies are shaded grey; negative anomalies are left white. Overlain are age contours in Ma for the crust (dashed red), pseudofaults PF5 and PF7 (solid red) and the location of fossil ridge axis FRA5 (dotted red). Included are the locations of the reflection (solid black) and refraction (dotted) lines with solid circles defining the OBS positions (refer to figure 5.1 for further identification). Green boxes illustrate the serpentinized regions inferred from the seismic lines. 5.3.2 Propagating rift history on the northern Juan de Fuca plate. Wilson et al. [1984] present a history of tectonic evolution of the Juan de Fuca plate by prop-agating rifting. Their two-plate spreading model successfully identified seven propagation se-quences, each successively rotating the Juan de Fuca ridge system in a clockwise direction (fig-ure 5.8). Propagating rift 5 first initiated at 8.56 Ma and ceased at ~ 5 Ma, generating pseu-dofaults on each side of the spreading axis. Pseudofault 5 (PF5) results from the southward propagation of the new spreading direction and is recorded by the disruption of the magnetic Chapter 5. DISCUSSION AND CONCLUSIONS 90 anomaly pattern on the northern Juan de Fuca plate (figure 5.7). Davis and Currie [1993] note a significant pole rotation shift corresponding to Chron 3a (about 5 Ma [ Wilson, 1993]) that involved at least 6 0 of clockwise rotation over a period of less than 1 million years. Wilson et al. [1984] suggest the change in spreading direction required a reorganization of the ridge that led to propagating rift events 6 and 7 at 4.5 Ma (figure 5.8). Pseudofault 7 (PF7) (figure 5.7) is a remnant of the northward propagating rift 7 event and appears to terminate at about 1.5 Ma where it intersects the Nootka fault zone. 5.0+ 5.0" 4.5 4.0 3.5 3.0 Figure 5.8: Reorientation of the Juan de Fuca ridge system: 5 Ma to 3 Ma. Propagating rift sequences are numbered. Reproduced from Wilson et al. [1984]. 5.3.3 Effects of ridge propagation on the crustal formation process From the close proximity and sub-parallel alignment of PF7 to the VISP II profile (figure 5.7), the development of PF7 may have been influential in the formation of the inferred zone of ser-pentinization. From about 4 Ma to 1.5 Ma new crust was being generated against existing crust east of PF7. The age contrast across most of PF7 is an approximately uniform 2 million years. Chapter 5. DISCUSSION AND CONCLUSIONS 91 This age difference has implications for the structure of the crust formed on the younger side of PF7. The relatively cold edge of lithosphere that faces the spreading ridge axis will perturb temperatures in the wedge of asthenosphere rising beneath the propagating ridge tip as it passes through older, cooler lithosphere. This is a well documented process at ridge-transform-ridge plate boundaries [ Fox and Gallo, 1984]. Fox and Gallo's [1984] generalized model for the structure of oceanic lithosphere proximal to a ridge-transform offset predicts that this effect will be more pronounced for slowly-slipping ridge-transform intersections (~ 2 cm/yr slip) than for fast-slipping ones (~ 18 cm/yr slip). The effectively thinned crust of the new plate may be ob-served at distances up to 20 km away from the edge of the colder lithosphere. In a refinement to the original propagating rift model, Hey et al. [1980] predict that there should be rough topography in the transition zones along pseudofaults because of a period of slow spreading at the tip of a propagating rift (figure 5.9). This is recorded at the Galapogas 95.5 ° W propagator where ~ 21 km of irregular rift tip spreading precedes normal spreading at the trailing ridge axis [ Kleinrock and Hey, 1989]. On the Juan de Fuca plate, this is observed on reflection line 85-07 over a region that coincides with the crossing of PF7 (figure 5.3). Similarly, a zone of rough basement topography should be apparent where VISP II intersects PF5, an inner pseudofault. Basement relief for the velocity model is coarsely sampled but does exhibit a zone of more variable topography between the fossil ridge axis intersection (FRA5) at 40 km and the crossing of PF5 at 60 km model distance (figure 5.10). This region should correspond to a predicted fossil shear zone that formed as the propagating rift moved away from the doomed rift (figure 5.9). 5.3.4 Implications of the initiation of the Nootka fault zone Riddihough [1984] concluded that the Explorer plate fragmented from the Juan de Fuca plate at about 4 Ma. From figure 5.7, the Juan de Fuca lithosphere east of PF7 appears to have formed Chapter 5. DISCUSSION AND CONCLUSIONS 92 Figure 5.9: Schematic of transitional crust created at a propagating rift tip. Active spreading . axes with full spreading rate are shown as heavy lines; active axes with transitional spreading rates are shown as dashed lines; fossil axis is shown as a dotted line. A transitional spreading zone exists between full spreading on the offset ridge axes: from left to right, the spreading rate decreases from the full spreading rate to zero on the dying spreading centre, while it increases from zero to the full rate on the propagating rift. New crust formed at the propagating rift tip should be characterized by rough basement topography associated with a decreased spreading rate. A shear zone should exist where the propagating rift tip moves away from the dying rift tip (shaded grey). Adapted from Hey et al. [1980]. prior to, or contemporaneously with, the initiation of the Nootka fault zone. Juan de Fuca plate west of PF7 formed post 2 Ma, well after left-lateral movement began within the fault zone. Juan de Fuca crust immediately south of the Nootka fault zone and west of PF7 may therefore be characterized by decreased melt supply as the new plate was forming at the spreading ridge, adjacent to the truncating Explorer plate. In a reconstruction of the spreading history of the Juan de Fuca ridge system, Wilson [1993] noted that magnetic anomalies 3 and 3A (4-6.5 Ma) south of the Nootka fault zone and north of Chapter 5. DISCUSSION AND CONCL USIONS 93 SE Model Distance (km) NW S ° ( J 20 40 1 1 60 80 1 1 100 120 I FRA5 PF5 Seawater Sediments Basement Figure 5.10: Basement topography beneath the VISP II profile determined from continuous seismic profiling. The most variable topography occurs over a zone between pseudofault 5 (PF5) and the fossil ridge axis (FRA5). PF5 (figure 5.7) did not fit a model of rigid plate motion. Wilson [1993] postulated that slight deformation of this area must have occurred either as the Explorer plate started to fragment from the Juan de Fuca plate at 4 Ma or, after the initial separation, as the diffuse plate boundary mi-grated across this area to its present position at the Nootka fault zone. Faulting associated with this distributed deformation may represent potential pathways for sea water penetration down to upper mantle depths and could explain why the inferred serpentinization is largely confined to this region of the northern Juan de Fuca plate. 5.4 Serpentinization model for the northern Juan de Fuca plate Figure 5.11 illustrates the proposed serpentinization model for the northern Juan de Fuca plate. The model shows two distinct regions of serpentinization within lithosphere of different ages. A narrow zone of hydrated lithosphere between 2 Ma and 4 Ma in age is inferred to coincide with a margin on the inboard side of outer pseudofault 7. This correlates with the observations of travel time delays on the forward and reverse sections of refraction line EX3 (section 5.1.1) and the zone of anomalous reflectivity on reflection line 85-07 (section 5.1.2). Crust formed in this region should be symptomatic of a slower accretion process associated with slower spread-ing as the propagating rift tip penetrated existing older crust formed at the dying rift (figure Chapter 5. DISCUSSION AND CONCLUSIONS 94 5.9). Rough topography predicted to form at this zone of slow-spreading is observed where the margin crosses line 85-07 (figure 5.3). Following O'Hartley's [1996] model for serpentinzation of peridotite at a divergent margin, outlined in section 5.2.2, the intermittent supply of magma to a slow-spreading ridge implies that the ridge will experience amagmatic episodes. Without magma to crystallize onto the spreading plates, the crust will thin by normal faulting and with-out a heat source, isotherms will collapse downward. Sub-500 ° C temperatures in the upper mantle allow serpentinization to proceed and the normal faults provide conduits for penetration of the mantle by sea water. An interlayered zone of ultramafic rocks and serpentinites beneath thinner than average crust would explain anomalous sub-crustal reflectivity observed beneath PF7 on line 85-07 (figure 5.3). The hydration of adjacent lithosphere aged 4 Ma to 6 Ma may have occurred over the same time interval (2-4 Ma). Upper mantle juxtaposed against newly created lithosphere will have undergone 2 million years of cooling by the time the propagating ridge passed along the trace of PF7. Temperatures within the upper mantle of the older lithosphere would have been conducive for serpentine formation and susceptible to percolating hydrothermal waters from the proximal spreading ridge. Initial separation of the Explorer plate from the Juan de Fuca plate occurred at ~ 4 Ma, prior to the time the rift was propagating through this area. A diffuse plate bound-ary or the migration of the Nootka fault system northward to its present location deformed the plate in the vicinity of the 4-6 Ma lithosphere [ Wilson, 1993]. Fracturing associated with this deformation may have provided additional pathways for circulating sea water to penetrate the upper mantle at distances up to 20 km from the passing spreading ridge. Propagating rift 7 had reached the present-day position of the Nootka fault zone by 1.5 Ma. Active deformation within a wide fault zone would have provided numerous pathways for sea water to penetrate the upper mantle at this time This may explain why the greatest thickness of serpentinized peridotite, in-ferred from the model of VISP II, occurs close to the present-day position of the Nootka fault zone (figure 4.8). Chapter 5: DISCUSSION AND CONCLUSIONS 95 129°W 128°W 127°W 49"N 48°N 129°W 128"W 127°W 2-4 Ma serpentinized lithosphere J ^ E E E 4-6 Ma serpentinized lithosphere Figure 5.11: Model for serpentinization of the northern Juan de Fuca plate. Included are pseudo-fault traces PF5 and PF7 (solid grey) and fossil ridge axis FRA5 (dotted grey). Also shown are the locations of seismic refraction lines (dotted black), OBS positions (solid black circles) and seismic reflection line 85-07 (thick black line). The zone of 2-4 Ma serpentinized lithosphere defines the margin of slow spreading crust on the inboard side of outer pseudofault PF7 which is predicted to accompany rift tip propagation. Crust either side of this margin is predicted to have formed at spreading rates comparable to the present-day Juan de Fuca ridge system (4-6 cm/yr) [ Riddihough and Hyndman, 1991]. Chapter 5. DISCUSSION AND CONCLUSIONS 96 5.4.1 Are propagating rift events elsewhere associated with serpentinization of the upper mantle? The formation of anomalous lithosphere at a pseudofault is not unique to the northern Juan de Fuca plate [ Kruse et al, 2000]. Likewise, the scenario of lithosphere hydration and ser-pentinization should be observed at analogous rift propagation events elsewhere. Kruse et al. [2000] show positive mantle Bouguer anomalies coinciding with an outer pseudofault trace on the eastern margin of the Juan Fernandez microplate. From this observation they postulate that outer pseudofaults may be associated with either thin and/or dense crust. In addition, the grav-ity signatures show similarities with those observed over comparable age offset nontransform discontinuities on the slow-spreading mid-Atlantic ridge. From a geophysical study at the mid-Atlantic ridge, Canales et al. [2000] suggest that nontransform offsets between slow-spreading ridge segments provide suitable environments for lithosphere hydration and serpentinization. If the slow-spreading crust formed near outer pseudofaults exhibits similar crustal accretion char-acteristics as slow-spreading crust formed at nontransform offsets, it follows that upper mantle near a pseudofault may be equally susceptible to serpentine formation. The Cobb propagator initiated at the Juan de Fuca ridge system at 4 Ma (Propagator 4 on figure 5.8). This event propagated northward through seafloor with a ~ 1.5 million year age contrast and is recorded by outer pseudofault 4 (PF4) on the Juan de Fuca plate. Where reflec-tion line 85-09 crosses the trace of PF4, Hasselgren and Clowes [1995] note a ~ 20 km wide zone of interlayered dipping reflectors immediately east of PF4 that is. remarkably similar to that observed beneath PF7 on line 85-07. However, in contrast to the zone observed on line 85-07, the reflectivity is apparent only beneath crust on the older side of the pseudofault. Future studies of oceanic lithosphere formed at pseudofaults that integrate seismic and grav-ity data may confirm whether serpentinization is a common process associated with propagating rift events. Chapter 5. DISCUSSION AND CONCLUSIONS 97 5.5 Conclusions Two-dimensional forward and inverse modelling of seismic refraction data reveals the presence of anomalously low velocities within the northern Juan de Fuca plate. The zone of low veloc-ities is modelled beneath the VISP II profile to extend for approximately 100 km and resides at depths typically associated with seismic layer 3B and the upper mantle of normal oceanic lithosphere. Synthetic seismograms generated for first arrival phases from the model are in good agreement with the recorded sections and confirm the existence of a lateral velocity gradient at sub-sediment depths of 3 km to 10 km. Observations of P- to S-wave phase conversions at the sediment/basement interface show that the overlying sediments are heavily saturated with Poisson's ratios of 0.45 to 0.46. The observation of a limited S-wave phase that passes through the upper crust indicates Poisson's ratios of 0.25 to 0.30, equivalent to measurements of typ-ical oceanic basalts. The completed velocity model suggests that sub-sediment oceanic crust increases in thickness from 7 km to 10 km from southeast to northwest (a thickness increase of 43%). A two-dimensional density model, constructed from the conversion of P-wave velocities us-ing standard velocity/density relationships for oceanic crust, produces a large (25 mgal) gravity low with a wavelength of ~ 140 km that is not evident in the free-air gravity anomaly. The in-corporation of gravity data discounts the thick crust hypothesis and indicates the existence of a zone of upper mantle containing 25-30% serpentine by volume. The modelled zone of partially serpentinized peridotite has P-wave velocities of 7-7.2 km/s and densities of 3050-3150 kg/m3. Exclusive of this region, seismic refraction and gravity modelling reveal seismic velocities and densities typical of normal oceanic crust. Observations from proximal seismic refraction lines EX 1,EX2, EX3 [Au and Clowes, 1982] and reflection line 85-07 [ Hasselgren and Clowes, 1995] confirm the existence of unusual ve-locity structure and imply that the inferred region of serpentinized upper mantle extends up to Chapter 5. DISCUSSION AND CONCLUSIONS 98 distances of 20 km away from the VISP II profile but is confined to the Juan de Fuca plate. The region of serpentinized upper mantle appears to coincide with 2-6 Ma aged crust on either side of a pseudofault trace that terminates at the Nootka fault zone. The process of serpentinization is interpreted to have resulted from a northward propagating rift that juxtaposed lithosphere with a 2 million year age contrast. Slow accretion rates at the propagating rift tip between. 2-4 Ma were conducive to lithosphere hydration and subsequent serpentinization which occurs over a zone on the inboard side of the trailing outer pseudofault. Upper mantle temperatures within older lithosphere on the outboard side of the propagating rift cooled to within the serpentine stability field at the time of rifting. The formation of the Nootka fault zone pervasively fractured 4-6 Ma aged crust on the outboard side of the propagating rift tip allowing hydrothermal currents that emanated at the propagating rift to efficiently percolate through the crust and serpentinize the underlying upper mantle. The inferred serpentinization regime is analogous to that at a slow-spreading nontransform offset [ Canales et al., 2000]. References Amor, J., 1996, Generalized software for seismic refraction data, LSPF Newsletter, 23: Lithop-robe Seismic Processing Facility, University of Calgary, Alberta. Atwater, T., 1989, Plate tectonic history of the northeast Pacific and western North America in Winterer, E. L., Hussong, D. M., and Decker, R. W., Eds., The Geology of North America, The Eastern Pacific Ocean and Hawaii: Geological Society of America, N, 21-72. Au, D., and Clowes, R. M., 1982, Crustal structure from an OBS survey of the Nootka fault zone off western Canada: Geophys. J.Roy. Astr. Soc, 68, 27-47. Barton, P. J., 1986, The relationship between seismic velocity and density in the continental crust - a useful constraint?: Geophys. J. R. astr. Soc, 87, 195-208. Best, A. I., Huggett, Q. J., and Harris, J. K., 2001, Comparison of in situ and laboratory acoustic measurements on Lough Hyne marine sediments: J. Acoust. Soc. Am., 110, 695-709. Braunmiller, J., and Nabelek, J., 2002, Seismotectonics of the Explorer region: / . Geophys. Res., 107, 2208, doi: 10.1029/2001JB000220. Briggs, I. C , 1974, Machine contouring using minimum curvature: Geophysics, 39, 39^ 48. Calvert, A. J., Hasselgren, E. A., and Clowes, R. M., 1990, Oceanic rift propagation - A cause of crustal underplating and seamount volcanism: Geology, 18, 886-889. Canales, J. P., Detrick, R. S., Lin, J., and Collins, J. A., 2000, Crustal and upper mantle seis-mic structure beneath the rift mountains and across a nontransform offset at the Mid-Atlantic Ridge (35 0 N): J. Geophys. Res., 105, 2699-2719. 99 References 100 Carlson, R. L., and Raskin, G. S., 1984, Density of the ocean crust: Nature, 311, 555-558. Cassidy, J. R, Ellis, R. M., Karavas, C , and Rogers, G. C , 1998, The northern limit of the subducted Juan de Fuca plate system: J. Geophys. Res., 103, 26,949-26,961. Christensen, N.I., 1966, Elasticity of ultrabasic rocks: J. Geophys. Res., 71, 5921-5931. Christensen, N. I., 1972, The abundance of serpentinites in the oceanic crust: J. Geology, 80, 709-719. Christensen, N. I., 1996, Poisson's ratio and crustal seismology: J. Geophys. Res., 101, 3139— 3156. Clowes, R. M., Yorath, C. J., and Hyndman, R. D., 1987, Reflection mapping across the con-vergent margin of western Canada: Geophys. J. R. astr. Soc, 89, 79-84. Clowes, R. M., Baird, D. J., and Dehler, S. A., 1997, Crustal structure of the Cascadia subduc-tion zone, southwestern British Columbia, from potential field and seismic studies: Can. J. Earth Sci., 34, 317-335. Cudrak, C. E , and Clowes, R. M., 1993, Crustal structure of Endeavour Ridge segment, Juan de Fuca Ridge, from a detailed seismic refraction survey: J. Geophys. Res., 98, 6329-6349. Davis, E. E., and Currie, R. G., 1993, Geophysical observations of the northern Juan de Fuca Ridge system: lessons in sea-floor spreading: Can. J. Earth Sci., 30, 278-300. Davis, E. E., Chapman, D. S., Forster, C. B., and Villinger, H., 1989, Heat-flow variations cor-related with buried basement topography on the Juan de Fuca Ridge flank: Nature, 342,533-537. References 101 DeMets, C , Gordon, R. G., Argus, D. R, and Stein, S., 1994, Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions: Geophys. Res. Lett., 21,2191-2194. Ellis, R. M., Spence, G. D., Clowes, R. M., Waldron, D. A., Jones, I. E , Green, A. G., Forsyth, D. A., Mair, J. A., Berry, M. J., Mereu, R. F., Kanasewich, E. R., Cumming, G. L., Hajnal, Z., Hyndman, R. D., McMechan, G. A., and Loncarevic, B. D., 1983, The Vancouver Is-land Seismic Project: a CO-CRUST onshore-offshore study of a convergent margin: Can. J. Earth Sci., 20,719-741. Ewing, J., and Houtz, R., 1979, Acoustic stratigraphy and structure of the oceanic crust in Tal-wani, M., Harrison, C. G. A., and Hayes, D. E., Eds., Deep drilling results in the Atlantic Ocean: American Geophysical Union, 1-14. Fox, P. J., and Gallo, D. G., 1984, A tectonic model for ridge-transform-ridge plate boundaries: implications for the structure of oceanic lithosphere: Tectonophysics, 104, 205-242. Francis, T. J. G., 1981, Serpentinization faults and their role in the tectonics of slow spreading ridges: J. Geophys. Res., 86, 11,616-11,622. Fryer, P., 2002, Recent studies of serpentinite occurences in the oceans: mantle-ocean interac-tions in the plate tectonic cycle: Chemie der Erde-Geochemistry, 62, 257-302. Hamilton, E. L., 1978, Sound velocity-density relations in seafloor sediments and rocks: Journal of the Acoustical Society of America, 63, 366-377. Hasselgren, E. O., and Clowes, R. M., 1995, Crustal structure of northern Juan de Fuca plate from multichannelreflection data: J. Geophys. Res., 100, 6469-6486. Hasselgren, E. O., Clowes, R. M., and Calvert, A. J., 1992, Propagating rift pseudofaults - zones References 102 of crustal underplating imaged by multichannel seismic reflection data: Geophys. Res. Lett., 19, 485-488. Hess, H. H., 1962, History of ocean basins in Engel, A. E. J., James, H. L., and Leonard, B. R, Eds., Petrologic studies: a volume in honor of A. F. Buddington: Geological Society of America, New York, 599-620. Hey, R., Duennebier, F. K., and Morgan, W. J., 1980, Propagating rifts on midocean ridges: J. Geophys. Res., 85, 3647-3658. Hey, R. N., Menard, H. W., Atwater, T. M., and Caress, D. W., 1988, Changes in direction of seafloor spreading revisited: J. Geophys. Res., 93, 2803-2811. Hey, R., 1977, A new class of 'pseudofaults' and their bearing on plate tectonics: a propagating rift model: Earth Planet. Sci. Lett., 37, 321-325. Hill, M. N., 1957, Recent geophysical exploration of the ocean floor: Physics and Chemistry of the Earth, 2, 129-163. Hyndman, R. D., Riddihough, R. P., and Herzer, R., 1979, The Nootka Fault Zone - a new plate boundary off western Canada: Geophys. J. Roy. Astr. Soc, 58, 667-683. Hyndman, R. D., Yorath, C. J., Clowes, R. M., and Davis, E. E., 1990, The northern Cascadia subduction zone at Vancouver Island: seismic structure and tectonic history: Can. J. Earth Sci., 27, 313-329. Jacobson, R. S., and Lewis, B. T. R., 1990, The first direct measurements of upper oceanic crustal compressional wave attenuation: J. Geophys. Res., 95, 17,417-17,429. Johnson, H. P., and Semyan, S. W., 1993, Age variation in the physical properties of oceanic basalts: Implications for crustal formation and evolution: J. Geophys. Res., 99, 3123-3134. References 103 Jones, E. J. W., 1999, Marine geophysics: John Wiley and Sons Ltd, England. Keen, C. E., and Barrett, D. L., 197.1, A measurement of seismic anisotropy in the northeast Pacific: Can. J. Earth Sci., 8, 1056-1064. Kleinrock, M. C , and Hey, R. N., 1989, Detailed tectonics near the Galapogas 95.5 ° W propag-ator: how the lithosphere tears and a spreading axis develops: J. Geophys. Res., 94, 13,801-13,838. Kreemer, C , Govers, R., Furlong, K. P., and Holt, W. E., 1998, Plate boundary deformation between the Pacific and North America in the Explorer region: Tectonophysics, 293, 225-238. Kruse, S. E., Tebbens, S. R, Naar, D. E , Lou, Q., and Bird, R. T., 2000, Comparisons of gravity anomalies at pseudofaults, fracture zones, and nontransform discontinuities from fast to slow spreading areas: J. Geophys. Res., 105, 28,399-28,410. Ludwig, J. W, Nafe, J. E., and Drake, C. L., 1970, Seismic refraction in Maxwell, A. E., Ed., The Sea: Wiley, 53-84. MacDonald, A. H., and Fyfe, W. S., 1985, Rate of serpentinization in seafloor environments: Tectonophysics, 116, 123-135. Malecek, S. J., and Clowes, R. M., 1978, Crustal structure near Explorer Ridge from a marine deep seismic sounding survey: J. Geophys. Res., 83, 5899-5912. Menke, W, West, M., and Tolstoy, M., 2002, Shallow-crustal magma chamber beneath the axial high of the Coaxial segment of Juan de Fuca Ridge at the source site of the 1993 eruption: Geology, 30, 359-362. References 104 O'Hanley, D. S., 1996, Serpentinites: records of tectonic and petrological history: Oxford Monographs on Geology and Geophysics, 34, 277p. Raff, A. D., and Mason, R. G., 1961, Magnetic survey off the west coast of North America, 400 N latitude to 50 ° N latitude: Geol. Soc. Am. Bull., 72, 1267-1270. Riddihough, R., and Hyndman, R. D., 1989, Queen Charlotte Islands margin in Winterer, E. L., Hussong, D. M., and Decker, R. W., Eds., The Geology of North America, The Eastern Pa-cific Ocean and Hawaii: Geological Society of America, N, 403-411. Riddihough, R. R, and Hyndman, R. D., 1991, Modern plate tectonic regime of the continen-tal margin of western Canada in Gabrielse, H., and Yorath, C. J., Eds., Geology of the Cor-dilleran Orogen in Canada: Geological Survey of Canada, Geology of Canada, 4, 435-455. Riddihough, R. P., 1979, Gravity arid structure of an active margin - British Columbia and Wash-ington: Can. J. Earth Sci., 16, 350-363. Riddihough, R., 1984, Recent movements of the Juan de Fuca plate system: J. Geophys. Res., 89,6980-6994. Rohr, K. M. M., and Furlong, K. P., 1995, Ephemeral plate tectonics at the Queen Charlotte triple junction: Geology, 23, 1035-1038. Salisbury, M. H., and Christensen, N. I., 1978, The seismic velocity structure of a traverse through the Bay of Islands ophiolite complex, Newfoundland, an exposure of oceanic crust and upper mantle: J. Geophys. Res., 83, 805-817. Salisbury, M. H., and Keen, C. E., 1993, Listric faults imaged in oceanic crust: Geology, 21, 117-120. Shearer, P. M., 1999, Introduction to seismology: Cambridge University Press, England. References 105 Sheriff, R. E., and Geldart, L. R, 1983, Exploration seismology, vol. 2: Data-processing and interpretation: Cambridge University Press, England. Spence, G. D., Clowes, R. M., and Ellis, R. M., 1985, Seismic structure across the active sub-duction zone of western Canada: J. Geophys. Res., 90, 6754-6772. Spudich, P., and Orcutt, J., 1980, Petrology and porosity of an oceanic crustal site: Results from waveform modeling of seismic refraction data: J. Geophys. Res., 85, 1409-1433. Talwani, M., Worzel, J. L., and Landisman, M., 1959, Rapid gravity computation for two-dimensional bodies with application to the Mendocino submarine fracture zone: J. Geophys. Res., 6 4 , 49-59. Turcotte, D. L . , and Schubert, G. S., 1982, Geodynamics: John Wiley and Sons, Inc., U.S.A. Underwood, M. B., and Hoke, K. D., 2000, Composition and provenance of turbidite sand and hemipelagic mud in northwestern Cascadia Basin in Fisher, A., Davis, E. E., and Escutia, C , Eds., Proceedings of the Ocean Drilling Program, Scientific Results: Ocean Drilling Pro-gram, 51-65. Waldron, D. A., Clowes, R. M., and White, D. J., 1990, Seismic structure of a subducting oceanic plate off western Canada in Green, A. G., Ed., Studies of Laterally Heterogeneous Structures Using Seismic Refraction and Reflection Data: Geological Survey of Canada, 105-113. Webring, M., 1985, SAKI: A fortran program for generalized linear inversion of gravity and magnetic profiles, Open file report: United States Department of the Interior, Geological Sur-vey, 85-122. West, M., Menke, W., and Tolstoy, M., 2003, Focused magma supply at the intersection References 106 of the Cobb hotspot and the Juan de Fuca ridge: Geophys. Res. Lett., 30, 1724, doi: 10.1029/2003GL017104. White, D. J., and Clowes, R. M., 1990, Shallow crustal structure beneath the Juan de Fuca Ridge from 2-D seismic refraction tomography: Geophys. J. Int., 100, 349-367. White, R. S., McKenzie, D., and O'Nions, R. K., 1992, Oceanic crustal thickness from seismic measurements and rare earth element inversions: J. Geophys. Res., 97, 19,683—19,715. Wilson, D. S., Hey, R. N., and Nishimura, C , 1984, Propagation as a means of reorientation of the Juan de Fuca ridge: J. Geophys. Res., 89, 9215-9225. Wilson, D. S., 1993, Confidence intervals for motion and deformation of the Juan de Fuca plate: J. Geophys. Res., 98, 16053-16071. Zelt, C. A., and Ellis, R. M., 1988, Practical and efficient ray tracing in two-dimensional media for rapid traveltime and amplitude forward modelling: Can. J. Expl. Geophys., 2 4 , 16-31. Zelt, C. A., and Smith, R. B., 1992, Seismic traveltime inversion for 2-D crustal velocity struc-ture: Geophys. J. Int., 108, 16-34. Appendix A Calculation of gravity response from 2-D modelling The SAKI gravity modelling algorithm is based on the method of Talwani et al. [1959]. A complex two-dimensional-like body may be approximated by an n-sided polygon of con-stant density that extends to infinity in the y-direction (figure A.l). The vertical component of gravitational attraction, Ag, for the polygonal body at a location on the surface is equal to a line integral around the perimeter of the body: Ag = 2Gpjz89 (A.l) where G is the universal constant of gravitation and p is the density of the body; • X : : -- a -— P S Figure A. 1: Geometry of a two-dimensional 5 sided polygon used to approximate an irregularly shaped body. Adapted from Talwani etal. [1959]. From figure A.l this is: 107 Appendix A. Calculation of gravity response from 2-D modelling 108 2Gp I z89 = 2G fB z89 + f z8d + f° z86 + ... [* z88 J JA JB Jc JE (A.2) From figure A. 1 we have the following geometric relations for an arbitrary point Q on BC: z = (x — ai) tan 7; = x tan 9, (A.3) and Zi+l tan 7; = x i+1 — zi+\ %i v Zi+l ~ %i (AA) From equation A.3 we can eliminate x to get z a, tan 7i tan 9 tan 7J — tan 9 and I z89= fc " ^ t a n f l ^ JBC JB tan 7; —tan 0 an 7^  — tan 9 This integral becomes (A.5) (A .6) / JBC z80 = cii sin ficosji I (9i — 9i+i) + tan 7; In cos #i(tan 0i — tan 7^ ) cos #;+i(tan 9i+\ — tan7j)J Zi (A.7) where 9i = tan 1 —, Appendix A. Calculation of gravity response from 2-D modelling 109 7i = tan and = tan" 1 ( — ) (A.8) can now be computed in terms of the x's and the z's. The total vertical component of gravitational attraction Ag at P is the sum of the contributions from each side of the polygon: 5 Ag = 2Gp^2Zi (A.9) i=i This calculation is modified to account for the limited extent of a body in the y direction by using approximate end-corrections. Appendix B Conductive cooling model of the oceanic lithosphere At a spreading ridge, the surface plates are created from upwelling asthenosphere and move ho-rizontally away from the ridge (figure B.l). Seawater efficiently cools the surface of the plates. The base of the lithosphere is a rheological boundary between the rigid plate and ductily flowing asthenosphere. This boundary may be approximated by an isotherm equivalent to the temper-ature at which mantle rocks deform in the laboratory (Tm = 1300 ° [ Turcotte and Schubert, 1982]). This boundary will deepen as the cooling lithosphere moves away from the ridge. Spreading Ridge T = TS T = T„ Asthenosphere z = 0 Figure B.l: Schematic for model of conductively cooling lithosphere. The cooling can be approximated by a one-dimensional model of a uniform half-space with a fixed surface temperature (figure B.l). The model, adapted from the derivation of Turcotte and Schubert [1982], assumes only vertical heat flow and vertical temperature variation (the lateral 110 Appendix B. Conductive cooling model of the oceanic lithosphere 111 flow of heat from one column to the next is assumed small relative to the vertical heat flow). We start with an initial half-space with a uniform temperature: Tm - 13000 C. The surface tem-perature, Ts = 0 0 C, is held fixed in time. Davis et al. [1989] estimated basement surface tem-peratures of 600 C for crust aged 3 Ma east of the Juan de Fuca ridge system but a difference of 60 ° C for Ts has little effect on the conductive cooling profiles (figure B.2). We now model the cooling of a vertical column as it moves away from either side of the spreading centre. Let T(z, t) = temperature as a function of depth, z, beneath the surface of the igneous base-ment and time, t. The temperature at t = 0 is Tm = 13000 C and is constant throughout the half-space. The solution for T(z, t) is: T(z,t) = Tt + (Tm^T.)erf(ri) (B.l) where (B.2) is the error function of rj which is a dimensionless ratio that relates depth to the thermal diffusion distance: " = 27S ( B ' 3 ) K is the coefficient of thermal diffusivity which is approximately 25 km2/MY for the earth [ Turcotte and Schubert, 1982]. This simple approximation neglects convective cooling by hydrothermal circulation within the crust. Efficient hydrothermal circulation will lower the upper boundary of conductive cool-ing to a greater depth within the crust comensurate with the depth of fluid circulation. In ad-dition, the model does not account for the thermal effects of the less permeable sedimentary Appendix B. Conductive cooling model of the oceanic lithosphere 112 0 200 400 600 800 1000 1200 1400 Temperature (°C) Figure B.2: Conductive cooling profiles for 1, 2, 5 and 10 Ma aged lithosphere. Profiles with a basement surface temperature of Ts = 0 ° C are shown in black. Profiles with Ts = 60 ° C are shown in grey. blanket that commonly overlies oceanic crust. The development of a less permeable sedimen-tary blanket may inhibit convective cooling of the basement surface but Davis et al. [1989] found, in a heat-flow study of sediment draped topography of crust aged 3 Ma on the flanks of the Juan de Fuca ridge, that pore-fluid convection was sufficiently vigorous to maintain uniform temperatures at the basement surface. 

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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052760/manifest

Comment

Related Items